Actual source code: ex215.c
1: static char help[] = "Tests MatSolve(), MatSolveTranspose() and MatMatSolve() with SEQDENSE\n";
3: #include <petscmat.h>
5: int main(int argc, char **args)
6: {
7: Mat A, RHS, C, F, X;
8: Vec u, x, b;
9: PetscMPIInt size;
10: PetscInt m, n, nsolve, nrhs;
11: PetscReal norm, tol = PETSC_SQRT_MACHINE_EPSILON;
12: PetscRandom rand;
13: PetscBool data_provided, herm, symm, hpd;
14: MatFactorType ftyp;
15: PetscViewer fd;
16: char file[PETSC_MAX_PATH_LEN];
19: PetscInitialize(&argc, &args, (char *)0, help);
20: MPI_Comm_size(PETSC_COMM_WORLD, &size);
22: /* Determine which type of solver we want to test for */
23: herm = PETSC_FALSE;
24: symm = PETSC_FALSE;
25: hpd = PETSC_FALSE;
26: PetscOptionsGetBool(NULL, NULL, "-symmetric_solve", &symm, NULL);
27: PetscOptionsGetBool(NULL, NULL, "-hermitian_solve", &herm, NULL);
28: PetscOptionsGetBool(NULL, NULL, "-hpd_solve", &hpd, NULL);
30: /* Determine file from which we read the matrix A */
31: ftyp = MAT_FACTOR_LU;
32: PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), &data_provided);
33: if (!data_provided) { /* get matrices from PETSc distribution */
34: PetscStrcpy(file, "${PETSC_DIR}/share/petsc/datafiles/matrices/");
35: if (hpd) {
36: #if defined(PETSC_USE_COMPLEX)
37: PetscStrcat(file, "hpd-complex-");
38: #else
39: PetscStrcat(file, "spd-real-");
40: #endif
41: ftyp = MAT_FACTOR_CHOLESKY;
42: } else {
43: #if defined(PETSC_USE_COMPLEX)
44: PetscStrcat(file, "nh-complex-");
45: #else
46: PetscStrcat(file, "ns-real-");
47: #endif
48: }
49: #if defined(PETSC_USE_64BIT_INDICES)
50: PetscStrcat(file, "int64-");
51: #else
52: PetscStrcat(file, "int32-");
53: #endif
54: #if defined(PETSC_USE_REAL_SINGLE)
55: PetscStrcat(file, "float32");
56: #else
57: PetscStrcat(file, "float64");
58: #endif
59: }
61: /* Load matrix A */
62: #if defined(PETSC_USE_REAL___FLOAT128)
63: PetscOptionsInsertString(NULL, "-binary_read_double");
64: #endif
65: PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd);
66: MatCreate(PETSC_COMM_WORLD, &A);
67: MatLoad(A, fd);
68: PetscViewerDestroy(&fd);
69: MatConvert(A, MATSEQDENSE, MAT_INPLACE_MATRIX, &A);
70: MatGetSize(A, &m, &n);
73: /* Create dense matrix C and X; C holds true solution with identical columns */
74: nrhs = 2;
75: PetscOptionsGetInt(NULL, NULL, "-nrhs", &nrhs, NULL);
76: MatCreate(PETSC_COMM_WORLD, &C);
77: MatSetSizes(C, m, PETSC_DECIDE, PETSC_DECIDE, nrhs);
78: MatSetType(C, MATDENSE);
79: MatSetFromOptions(C);
80: MatSetUp(C);
82: PetscRandomCreate(PETSC_COMM_WORLD, &rand);
83: PetscRandomSetFromOptions(rand);
84: MatSetRandom(C, rand);
85: MatDuplicate(C, MAT_DO_NOT_COPY_VALUES, &X);
86: MatDuplicate(C, MAT_DO_NOT_COPY_VALUES, &RHS);
88: /* Create vectors */
89: VecCreate(PETSC_COMM_WORLD, &x);
90: VecSetSizes(x, n, PETSC_DECIDE);
91: VecSetFromOptions(x);
92: VecDuplicate(x, &b);
93: VecDuplicate(x, &u); /* save the true solution */
95: /* make a symmetric matrix */
96: if (symm) {
97: Mat AT;
99: MatTranspose(A, MAT_INITIAL_MATRIX, &AT);
100: MatAXPY(A, 1.0, AT, SAME_NONZERO_PATTERN);
101: MatDestroy(&AT);
102: ftyp = MAT_FACTOR_CHOLESKY;
103: }
104: /* make an hermitian matrix */
105: if (herm) {
106: Mat AH;
108: MatHermitianTranspose(A, MAT_INITIAL_MATRIX, &AH);
109: MatAXPY(A, 1.0, AH, SAME_NONZERO_PATTERN);
110: MatDestroy(&AH);
111: ftyp = MAT_FACTOR_CHOLESKY;
112: }
113: PetscObjectSetName((PetscObject)A, "A");
114: MatViewFromOptions(A, NULL, "-amat_view");
116: MatDuplicate(A, MAT_COPY_VALUES, &F);
117: MatSetOption(F, MAT_SYMMETRIC, symm);
118: /* it seems that the SPD concept in PETSc extends naturally to Hermitian Positive definitess */
119: MatSetOption(F, MAT_HERMITIAN, (PetscBool)(hpd || herm));
120: MatSetOption(F, MAT_SPD, hpd);
121: {
122: PetscInt iftyp = ftyp;
123: PetscOptionsGetEList(NULL, NULL, "-ftype", MatFactorTypes, MAT_FACTOR_NUM_TYPES, &iftyp, NULL);
124: ftyp = (MatFactorType)iftyp;
125: }
126: if (ftyp == MAT_FACTOR_LU) {
127: MatLUFactor(F, NULL, NULL, NULL);
128: } else if (ftyp == MAT_FACTOR_CHOLESKY) {
129: MatCholeskyFactor(F, NULL, NULL);
130: } else if (ftyp == MAT_FACTOR_QR) {
131: MatQRFactor(F, NULL, NULL);
132: } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Factorization %s not supported in this example", MatFactorTypes[ftyp]);
134: for (nsolve = 0; nsolve < 2; nsolve++) {
135: VecSetRandom(x, rand);
136: VecCopy(x, u);
137: if (nsolve) {
138: MatMult(A, x, b);
139: MatSolve(F, b, x);
140: } else {
141: MatMultTranspose(A, x, b);
142: MatSolveTranspose(F, b, x);
143: }
144: /* Check the error */
145: VecAXPY(u, -1.0, x); /* u <- (-1.0)x + u */
146: VecNorm(u, NORM_2, &norm);
147: if (norm > tol) {
148: PetscReal resi;
149: if (nsolve) {
150: MatMult(A, x, u); /* u = A*x */
151: } else {
152: MatMultTranspose(A, x, u); /* u = A*x */
153: }
154: VecAXPY(u, -1.0, b); /* u <- (-1.0)b + u */
155: VecNorm(u, NORM_2, &resi);
156: if (nsolve) {
157: PetscPrintf(PETSC_COMM_SELF, "MatSolve error: Norm of error %g, residual %g\n", (double)norm, (double)resi);
158: } else {
159: PetscPrintf(PETSC_COMM_SELF, "MatSolveTranspose error: Norm of error %g, residual %g\n", (double)norm, (double)resi);
160: }
161: }
162: }
163: MatMatMult(A, C, MAT_REUSE_MATRIX, 2.0, &RHS);
164: MatMatSolve(F, RHS, X);
166: /* Check the error */
167: MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN);
168: MatNorm(X, NORM_FROBENIUS, &norm);
169: if (norm > tol) PetscPrintf(PETSC_COMM_SELF, "MatMatSolve: Norm of error %g\n", (double)norm);
171: /* Free data structures */
172: MatDestroy(&A);
173: MatDestroy(&C);
174: MatDestroy(&F);
175: MatDestroy(&X);
176: MatDestroy(&RHS);
177: PetscRandomDestroy(&rand);
178: VecDestroy(&x);
179: VecDestroy(&b);
180: VecDestroy(&u);
181: PetscFinalize();
182: return 0;
183: }
185: /*TEST
187: testset:
188: output_file: output/ex215.out
189: test:
190: suffix: ns
191: test:
192: suffix: sym
193: args: -symmetric_solve
194: test:
195: suffix: herm
196: args: -hermitian_solve
197: test:
198: suffix: hpd
199: args: -hpd_solve
200: test:
201: suffix: qr
202: args: -ftype qr
204: TEST*/