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
main(int argc,char ** args)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