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"; 8 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 21 22 rhs 23 --- 24 0 3.0 25 1 12.0 26 2 6.0 27 28 solu 29 ---- 30 0 1.0 31 0 1.0 32 0 1.0 33 */ 34 35 #include <petscmat.h> 36 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; 50 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!"); 55 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)); 80 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 } 91 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 } 107 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 } 123 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)); 128 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 } 137 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 } 147 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 } 155 156 /*TEST 157 158 build: 159 requires: !defined(PETSC_USE_64BIT_INDICES) double !complex 160 161 test: 162 requires: datafilespath 163 args: -Ain ${DATAFILESPATH}/matrices/indefinite/afiro_A.dat 164 165 TEST*/ 166