Actual source code: ex78.c

  1: static char help[] = "Reads in a matrix in ASCII MATLAB format (I,J,A), read in vectors rhs and exact_solu in ASCII format.\n\
  2: Writes them using the PETSc sparse format.\n\
  3: Note: I and J start at 1, not 0, use -noshift if indices in file start with zero!\n\
  4: Input parameters are:\n\
  5:   -Ain  <filename> : input matrix in ascii format\n\
  6:   -rhs  <filename> : input rhs in ascii format\n\
  7:   -solu  <filename> : input true solution in ascii format\n\\n";

  9: /*
 10: Example: ./ex78 -Ain Ain -rhs rhs -solu solu -noshift -mat_view
 11:  with the datafiles in the following format:
 12: Ain (I and J start at 0):
 13: ------------------------
 14: 3 3 6
 15: 0 0 1.0
 16: 0 1 2.0
 17: 1 0 3.0
 18: 1 1 4.0
 19: 1 2 5.0
 20: 2 2 6.0

 22: rhs
 23: ---
 24: 0 3.0
 25: 1 12.0
 26: 2 6.0

 28: solu
 29: ----
 30: 0 1.0
 31: 0 1.0
 32: 0 1.0
 33: */

 35: #include <petscmat.h>

 37: int main(int argc, char **args)
 38: {
 39:   Mat         A = NULL;
 40:   Vec         b, u = NULL, u_tmp;
 41:   char        Ain[PETSC_MAX_PATH_LEN], rhs[PETSC_MAX_PATH_LEN], solu[PETSC_MAX_PATH_LEN];
 42:   int         m, n = 0, nz, dummy; /* these are fscaned so kept as int */
 43:   PetscInt    i, col, row, shift = 1, sizes[3], nsizes;
 44:   PetscScalar val;
 45:   PetscReal   res_norm;
 46:   FILE       *Afile, *bfile, *ufile;
 47:   PetscViewer view;
 48:   PetscBool   flg_A, flg_b, flg_u, flg;
 49:   PetscMPIInt size;

 51:   PetscFunctionBeginUser;
 52:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 53:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 54:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");

 56:   /* Read in matrix, rhs and exact solution from ascii files */
 57:   PetscCall(PetscOptionsGetString(NULL, NULL, "-Ain", Ain, sizeof(Ain), &flg_A));
 58:   PetscCall(PetscOptionsHasName(NULL, NULL, "-noshift", &flg));
 59:   if (flg) shift = 0;
 60:   if (flg_A) {
 61:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n Read matrix in ascii format ...\n"));
 62:     PetscCall(PetscFOpen(PETSC_COMM_SELF, Ain, "r", &Afile));
 63:     nsizes = 3;
 64:     PetscCall(PetscOptionsGetIntArray(NULL, NULL, "-nosizesinfile", sizes, &nsizes, &flg));
 65:     if (flg) {
 66:       PetscCheck(nsizes == 3, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must pass in three m,n,nz as arguments for -nosizesinfile");
 67:       m  = sizes[0];
 68:       n  = sizes[1];
 69:       nz = sizes[2];
 70:     } else {
 71:       PetscCheck(fscanf(Afile, "%d %d %d", &m, &n, &nz) == 3, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Badly formatted input file");
 72:     }
 73:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "m: %d, n: %d, nz: %d \n", m, n, nz));
 74:     PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of rows, cols must be same for this example");
 75:     PetscCall(MatCreate(PETSC_COMM_SELF, &A));
 76:     PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m, n));
 77:     PetscCall(MatSetFromOptions(A));
 78:     PetscCall(MatSeqAIJSetPreallocation(A, nz / m, NULL));
 79:     PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));

 81:     for (i = 0; i < nz; i++) {
 82:       PetscCheck(fscanf(Afile, "%d %d %le", &row, &col, (double *)&val) == 3, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Badly formatted input file");
 83:       row -= shift;
 84:       col -= shift; /* set index set starts at 0 */
 85:       PetscCall(MatSetValues(A, 1, &row, 1, &col, &val, INSERT_VALUES));
 86:     }
 87:     PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 88:     PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
 89:     fclose(Afile);
 90:   }

 92:   PetscCall(PetscOptionsGetString(NULL, NULL, "-rhs", rhs, sizeof(rhs), &flg_b));
 93:   if (flg_b) {
 94:     PetscCall(VecCreate(PETSC_COMM_SELF, &b));
 95:     PetscCall(VecSetSizes(b, PETSC_DECIDE, n));
 96:     PetscCall(VecSetFromOptions(b));
 97:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n Read rhs in ascii format ...\n"));
 98:     PetscCall(PetscFOpen(PETSC_COMM_SELF, rhs, "r", &bfile));
 99:     for (i = 0; i < n; i++) {
100:       PetscCheck(fscanf(bfile, "%d %le", &dummy, (double *)&val) == 2, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Badly formatted input file");
101:       PetscCall(VecSetValues(b, 1, &i, &val, INSERT_VALUES));
102:     }
103:     PetscCall(VecAssemblyBegin(b));
104:     PetscCall(VecAssemblyEnd(b));
105:     fclose(bfile);
106:   }

108:   PetscCall(PetscOptionsGetString(NULL, NULL, "-solu", solu, sizeof(solu), &flg_u));
109:   if (flg_u) {
110:     PetscCall(VecCreate(PETSC_COMM_SELF, &u));
111:     PetscCall(VecSetSizes(u, PETSC_DECIDE, n));
112:     PetscCall(VecSetFromOptions(u));
113:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n Read exact solution in ascii format ...\n"));
114:     PetscCall(PetscFOpen(PETSC_COMM_SELF, solu, "r", &ufile));
115:     for (i = 0; i < n; i++) {
116:       PetscCheck(fscanf(ufile, "%d  %le", &dummy, (double *)&val) == 2, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Badly formatted input file");
117:       PetscCall(VecSetValues(u, 1, &i, &val, INSERT_VALUES));
118:     }
119:     PetscCall(VecAssemblyBegin(u));
120:     PetscCall(VecAssemblyEnd(u));
121:     fclose(ufile);
122:   }

124:   /* Write matrix, rhs and exact solution in Petsc binary file */
125:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n Write matrix in binary to 'matrix.dat' ...\n"));
126:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_SELF, "matrix.dat", FILE_MODE_WRITE, &view));
127:   PetscCall(MatView(A, view));

129:   if (flg_b) { /* Write rhs in Petsc binary file */
130:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n Write rhs in binary to 'matrix.dat' ...\n"));
131:     PetscCall(VecView(b, view));
132:   }
133:   if (flg_u) {
134:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n Write exact solution in binary to 'matrix.dat' ...\n"));
135:     PetscCall(VecView(u, view));
136:   }

138:   /* Check accuracy of the data */
139:   if (flg_A & flg_b & flg_u) {
140:     PetscCall(VecDuplicate(u, &u_tmp));
141:     PetscCall(MatMult(A, u, u_tmp));
142:     PetscCall(VecAXPY(u_tmp, -1.0, b));
143:     PetscCall(VecNorm(u_tmp, NORM_2, &res_norm));
144:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n Accuracy of the reading data: | b - A*u |_2 : %g \n", res_norm));
145:     PetscCall(VecDestroy(&u_tmp));
146:   }

148:   PetscCall(MatDestroy(&A));
149:   if (flg_b) PetscCall(VecDestroy(&b));
150:   if (flg_u) PetscCall(VecDestroy(&u));
151:   PetscCall(PetscViewerDestroy(&view));
152:   PetscCall(PetscFinalize());
153:   return 0;
154: }

156: /*TEST

158:    build:
159:       requires: !defined(PETSC_USE_64BIT_INDICES) double !complex

161:    test:
162:       requires: datafilespath
163:       args: -Ain ${DATAFILESPATH}/matrices/indefinite/afiro_A.dat

165: TEST*/