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, (char *)0, 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\n", &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\n", &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\n", &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\n", &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*/