1 2 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\ 3 Writes them using the PETSc sparse format.\n\ 4 Note: I and J start at 1, not 0, use -noshift if indices in file start with zero!\n\ 5 Input parameters are:\n\ 6 -Ain <filename> : input matrix in ascii format\n\ 7 -rhs <filename> : input rhs in ascii format\n\ 8 -solu <filename> : input true solution in ascii format\n\\n"; 9 10 /* 11 Example: ./ex78 -Ain Ain -rhs rhs -solu solu -noshift -mat_view 12 with the datafiles in the following format: 13 Ain (I and J start at 0): 14 ------------------------ 15 3 3 6 16 0 0 1.0 17 0 1 2.0 18 1 0 3.0 19 1 1 4.0 20 1 2 5.0 21 2 2 6.0 22 23 rhs 24 --- 25 0 3.0 26 1 12.0 27 2 6.0 28 29 solu 30 ---- 31 0 1.0 32 0 1.0 33 0 1.0 34 */ 35 36 #include <petscmat.h> 37 38 int main(int argc, char **args) 39 { 40 Mat A = NULL; 41 Vec b, u = NULL, u_tmp; 42 char Ain[PETSC_MAX_PATH_LEN], rhs[PETSC_MAX_PATH_LEN], solu[PETSC_MAX_PATH_LEN]; 43 int m, n = 0, nz, dummy; /* these are fscaned so kept as int */ 44 PetscInt i, col, row, shift = 1, sizes[3], nsizes; 45 PetscScalar val; 46 PetscReal res_norm; 47 FILE *Afile, *bfile, *ufile; 48 PetscViewer view; 49 PetscBool flg_A, flg_b, flg_u, flg; 50 PetscMPIInt size; 51 52 PetscFunctionBeginUser; 53 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 54 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 55 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!"); 56 57 /* Read in matrix, rhs and exact solution from ascii files */ 58 PetscCall(PetscOptionsGetString(NULL, NULL, "-Ain", Ain, sizeof(Ain), &flg_A)); 59 PetscCall(PetscOptionsHasName(NULL, NULL, "-noshift", &flg)); 60 if (flg) shift = 0; 61 if (flg_A) { 62 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n Read matrix in ascii format ...\n")); 63 PetscCall(PetscFOpen(PETSC_COMM_SELF, Ain, "r", &Afile)); 64 nsizes = 3; 65 PetscCall(PetscOptionsGetIntArray(NULL, NULL, "-nosizesinfile", sizes, &nsizes, &flg)); 66 if (flg) { 67 PetscCheck(nsizes == 3, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must pass in three m,n,nz as arguments for -nosizesinfile"); 68 m = sizes[0]; 69 n = sizes[1]; 70 nz = sizes[2]; 71 } else { 72 PetscCheck(fscanf(Afile, "%d %d %d\n", &m, &n, &nz) == 3, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Badly formatted input file"); 73 } 74 PetscCall(PetscPrintf(PETSC_COMM_SELF, "m: %d, n: %d, nz: %d \n", m, n, nz)); 75 PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of rows, cols must be same for this example"); 76 PetscCall(MatCreate(PETSC_COMM_SELF, &A)); 77 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m, n)); 78 PetscCall(MatSetFromOptions(A)); 79 PetscCall(MatSeqAIJSetPreallocation(A, nz / m, NULL)); 80 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)); 81 82 for (i = 0; i < nz; i++) { 83 PetscCheck(fscanf(Afile, "%d %d %le\n", &row, &col, (double *)&val) == 3, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Badly formatted input file"); 84 row -= shift; 85 col -= shift; /* set index set starts at 0 */ 86 PetscCall(MatSetValues(A, 1, &row, 1, &col, &val, INSERT_VALUES)); 87 } 88 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 89 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 90 fclose(Afile); 91 } 92 93 PetscCall(PetscOptionsGetString(NULL, NULL, "-rhs", rhs, sizeof(rhs), &flg_b)); 94 if (flg_b) { 95 PetscCall(VecCreate(PETSC_COMM_SELF, &b)); 96 PetscCall(VecSetSizes(b, PETSC_DECIDE, n)); 97 PetscCall(VecSetFromOptions(b)); 98 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n Read rhs in ascii format ...\n")); 99 PetscCall(PetscFOpen(PETSC_COMM_SELF, rhs, "r", &bfile)); 100 for (i = 0; i < n; i++) { 101 PetscCheck(fscanf(bfile, "%d %le\n", &dummy, (double *)&val) == 2, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Badly formatted input file"); 102 PetscCall(VecSetValues(b, 1, &i, &val, INSERT_VALUES)); 103 } 104 PetscCall(VecAssemblyBegin(b)); 105 PetscCall(VecAssemblyEnd(b)); 106 fclose(bfile); 107 } 108 109 PetscCall(PetscOptionsGetString(NULL, NULL, "-solu", solu, sizeof(solu), &flg_u)); 110 if (flg_u) { 111 PetscCall(VecCreate(PETSC_COMM_SELF, &u)); 112 PetscCall(VecSetSizes(u, PETSC_DECIDE, n)); 113 PetscCall(VecSetFromOptions(u)); 114 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n Read exact solution in ascii format ...\n")); 115 PetscCall(PetscFOpen(PETSC_COMM_SELF, solu, "r", &ufile)); 116 for (i = 0; i < n; i++) { 117 PetscCheck(fscanf(ufile, "%d %le\n", &dummy, (double *)&val) == 2, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Badly formatted input file"); 118 PetscCall(VecSetValues(u, 1, &i, &val, INSERT_VALUES)); 119 } 120 PetscCall(VecAssemblyBegin(u)); 121 PetscCall(VecAssemblyEnd(u)); 122 fclose(ufile); 123 } 124 125 /* Write matrix, rhs and exact solution in Petsc binary file */ 126 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n Write matrix in binary to 'matrix.dat' ...\n")); 127 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_SELF, "matrix.dat", FILE_MODE_WRITE, &view)); 128 PetscCall(MatView(A, view)); 129 130 if (flg_b) { /* Write rhs in Petsc binary file */ 131 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n Write rhs in binary to 'matrix.dat' ...\n")); 132 PetscCall(VecView(b, view)); 133 } 134 if (flg_u) { 135 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n Write exact solution in binary to 'matrix.dat' ...\n")); 136 PetscCall(VecView(u, view)); 137 } 138 139 /* Check accuracy of the data */ 140 if (flg_A & flg_b & flg_u) { 141 PetscCall(VecDuplicate(u, &u_tmp)); 142 PetscCall(MatMult(A, u, u_tmp)); 143 PetscCall(VecAXPY(u_tmp, -1.0, b)); 144 PetscCall(VecNorm(u_tmp, NORM_2, &res_norm)); 145 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n Accuracy of the reading data: | b - A*u |_2 : %g \n", res_norm)); 146 PetscCall(VecDestroy(&u_tmp)); 147 } 148 149 PetscCall(MatDestroy(&A)); 150 if (flg_b) PetscCall(VecDestroy(&b)); 151 if (flg_u) PetscCall(VecDestroy(&u)); 152 PetscCall(PetscViewerDestroy(&view)); 153 PetscCall(PetscFinalize()); 154 return 0; 155 } 156 157 /*TEST 158 159 build: 160 requires: !defined(PETSC_USE_64BIT_INDICES) double !complex 161 162 test: 163 requires: datafilespath 164 args: -Ain ${DATAFILESPATH}/matrices/indefinite/afiro_A.dat 165 166 TEST*/ 167