Actual source code: ex78.c
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";
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
23: rhs
24: ---
25: 0 3.0
26: 1 12.0
27: 2 6.0
29: solu
30: ----
31: 0 1.0
32: 0 1.0
33: 0 1.0
34: */
36: #include <petscmat.h>
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;
53: PetscInitialize(&argc, &args, (char *)0, help);
54: MPI_Comm_size(PETSC_COMM_WORLD, &size);
57: /* Read in matrix, rhs and exact solution from ascii files */
58: PetscOptionsGetString(NULL, NULL, "-Ain", Ain, sizeof(Ain), &flg_A);
59: PetscOptionsHasName(NULL, NULL, "-noshift", &flg);
60: if (flg) shift = 0;
61: if (flg_A) {
62: PetscPrintf(PETSC_COMM_SELF, "\n Read matrix in ascii format ...\n");
63: PetscFOpen(PETSC_COMM_SELF, Ain, "r", &Afile);
64: nsizes = 3;
65: PetscOptionsGetIntArray(NULL, NULL, "-nosizesinfile", sizes, &nsizes, &flg);
66: if (flg) {
68: m = sizes[0];
69: n = sizes[1];
70: nz = sizes[2];
71: } else {
73: }
74: PetscPrintf(PETSC_COMM_SELF, "m: %d, n: %d, nz: %d \n", m, n, nz);
76: MatCreate(PETSC_COMM_SELF, &A);
77: MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m, n);
78: MatSetFromOptions(A);
79: MatSeqAIJSetPreallocation(A, nz / m, NULL);
80: MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
82: for (i = 0; i < nz; i++) {
84: row -= shift;
85: col -= shift; /* set index set starts at 0 */
86: MatSetValues(A, 1, &row, 1, &col, &val, INSERT_VALUES);
87: }
88: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
89: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
90: fclose(Afile);
91: }
93: PetscOptionsGetString(NULL, NULL, "-rhs", rhs, sizeof(rhs), &flg_b);
94: if (flg_b) {
95: VecCreate(PETSC_COMM_SELF, &b);
96: VecSetSizes(b, PETSC_DECIDE, n);
97: VecSetFromOptions(b);
98: PetscPrintf(PETSC_COMM_SELF, "\n Read rhs in ascii format ...\n");
99: PetscFOpen(PETSC_COMM_SELF, rhs, "r", &bfile);
100: for (i = 0; i < n; i++) {
102: VecSetValues(b, 1, &i, &val, INSERT_VALUES);
103: }
104: VecAssemblyBegin(b);
105: VecAssemblyEnd(b);
106: fclose(bfile);
107: }
109: PetscOptionsGetString(NULL, NULL, "-solu", solu, sizeof(solu), &flg_u);
110: if (flg_u) {
111: VecCreate(PETSC_COMM_SELF, &u);
112: VecSetSizes(u, PETSC_DECIDE, n);
113: VecSetFromOptions(u);
114: PetscPrintf(PETSC_COMM_SELF, "\n Read exact solution in ascii format ...\n");
115: PetscFOpen(PETSC_COMM_SELF, solu, "r", &ufile);
116: for (i = 0; i < n; i++) {
118: VecSetValues(u, 1, &i, &val, INSERT_VALUES);
119: }
120: VecAssemblyBegin(u);
121: VecAssemblyEnd(u);
122: fclose(ufile);
123: }
125: /* Write matrix, rhs and exact solution in Petsc binary file */
126: PetscPrintf(PETSC_COMM_SELF, "\n Write matrix in binary to 'matrix.dat' ...\n");
127: PetscViewerBinaryOpen(PETSC_COMM_SELF, "matrix.dat", FILE_MODE_WRITE, &view);
128: MatView(A, view);
130: if (flg_b) { /* Write rhs in Petsc binary file */
131: PetscPrintf(PETSC_COMM_SELF, "\n Write rhs in binary to 'matrix.dat' ...\n");
132: VecView(b, view);
133: }
134: if (flg_u) {
135: PetscPrintf(PETSC_COMM_SELF, "\n Write exact solution in binary to 'matrix.dat' ...\n");
136: VecView(u, view);
137: }
139: /* Check accuracy of the data */
140: if (flg_A & flg_b & flg_u) {
141: VecDuplicate(u, &u_tmp);
142: MatMult(A, u, u_tmp);
143: VecAXPY(u_tmp, -1.0, b);
144: VecNorm(u_tmp, NORM_2, &res_norm);
145: PetscPrintf(PETSC_COMM_SELF, "\n Accuracy of the reading data: | b - A*u |_2 : %g \n", res_norm);
146: VecDestroy(&u_tmp);
147: }
149: MatDestroy(&A);
150: if (flg_b) VecDestroy(&b);
151: if (flg_u) VecDestroy(&u);
152: PetscViewerDestroy(&view);
153: PetscFinalize();
154: return 0;
155: }
157: /*TEST
159: build:
160: requires: !defined(PETSC_USE_64BIT_INDICES) double !complex
162: test:
163: requires: datafilespath
164: args: -Ain ${DATAFILESPATH}/matrices/indefinite/afiro_A.dat
166: TEST*/