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];
18: PetscFunctionBeginUser;
19: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
20: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
21: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor test");
22: /* Determine which type of solver we want to test for */
23: herm = PETSC_FALSE;
24: symm = PETSC_FALSE;
25: hpd = PETSC_FALSE;
26: PetscCall(PetscOptionsGetBool(NULL, NULL, "-symmetric_solve", &symm, NULL));
27: PetscCall(PetscOptionsGetBool(NULL, NULL, "-hermitian_solve", &herm, NULL));
28: PetscCall(PetscOptionsGetBool(NULL, NULL, "-hpd_solve", &hpd, NULL));
30: /* Determine file from which we read the matrix A */
31: ftyp = MAT_FACTOR_LU;
32: PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), &data_provided));
33: if (!data_provided) { /* get matrices from PETSc distribution */
34: PetscCall(PetscStrncpy(file, "${PETSC_DIR}/share/petsc/datafiles/matrices/", sizeof(file)));
35: if (hpd) {
36: #if defined(PETSC_USE_COMPLEX)
37: PetscCall(PetscStrlcat(file, "hpd-complex-", sizeof(file)));
38: #else
39: PetscCall(PetscStrlcat(file, "spd-real-", sizeof(file)));
40: #endif
41: ftyp = MAT_FACTOR_CHOLESKY;
42: } else {
43: #if defined(PETSC_USE_COMPLEX)
44: PetscCall(PetscStrlcat(file, "nh-complex-", sizeof(file)));
45: #else
46: PetscCall(PetscStrlcat(file, "ns-real-", sizeof(file)));
47: #endif
48: }
49: #if defined(PETSC_USE_64BIT_INDICES)
50: PetscCall(PetscStrlcat(file, "int64-", sizeof(file)));
51: #else
52: PetscCall(PetscStrlcat(file, "int32-", sizeof(file)));
53: #endif
54: #if defined(PETSC_USE_REAL_SINGLE)
55: PetscCall(PetscStrlcat(file, "float32", sizeof(file)));
56: #else
57: PetscCall(PetscStrlcat(file, "float64", sizeof(file)));
58: #endif
59: }
61: /* Load matrix A */
62: #if defined(PETSC_USE_REAL___FLOAT128)
63: PetscCall(PetscOptionsInsertString(NULL, "-binary_read_double"));
64: #endif
65: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd));
66: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
67: PetscCall(MatLoad(A, fd));
68: PetscCall(PetscViewerDestroy(&fd));
69: PetscCall(MatConvert(A, MATSEQDENSE, MAT_INPLACE_MATRIX, &A));
70: PetscCall(MatGetSize(A, &m, &n));
71: PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")", m, n);
73: /* Create dense matrix C and X; C holds true solution with identical columns */
74: nrhs = 2;
75: PetscCall(PetscOptionsGetInt(NULL, NULL, "-nrhs", &nrhs, NULL));
76: PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
77: PetscCall(MatSetSizes(C, m, PETSC_DECIDE, PETSC_DECIDE, nrhs));
78: PetscCall(MatSetType(C, MATDENSE));
79: PetscCall(MatSetFromOptions(C));
80: PetscCall(MatSetUp(C));
82: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand));
83: PetscCall(PetscRandomSetFromOptions(rand));
84: PetscCall(MatSetRandom(C, rand));
85: PetscCall(MatDuplicate(C, MAT_DO_NOT_COPY_VALUES, &X));
86: PetscCall(MatDuplicate(C, MAT_DO_NOT_COPY_VALUES, &RHS));
88: /* Create vectors */
89: PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
90: PetscCall(VecSetSizes(x, n, PETSC_DECIDE));
91: PetscCall(VecSetFromOptions(x));
92: PetscCall(VecDuplicate(x, &b));
93: PetscCall(VecDuplicate(x, &u)); /* save the true solution */
95: /* make a symmetric matrix */
96: if (symm) {
97: Mat AT;
99: PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &AT));
100: PetscCall(MatAXPY(A, 1.0, AT, SAME_NONZERO_PATTERN));
101: PetscCall(MatDestroy(&AT));
102: ftyp = MAT_FACTOR_CHOLESKY;
103: }
104: /* make an hermitian matrix */
105: if (herm) {
106: Mat AH;
108: PetscCall(MatHermitianTranspose(A, MAT_INITIAL_MATRIX, &AH));
109: PetscCall(MatAXPY(A, 1.0, AH, SAME_NONZERO_PATTERN));
110: PetscCall(MatDestroy(&AH));
111: ftyp = MAT_FACTOR_CHOLESKY;
112: }
113: PetscCall(PetscObjectSetName((PetscObject)A, "A"));
114: PetscCall(MatViewFromOptions(A, NULL, "-amat_view"));
116: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &F));
117: PetscCall(MatSetOption(F, MAT_SYMMETRIC, symm));
118: /* it seems that the SPD concept in PETSc extends naturally to Hermitian Positive definitess */
119: PetscCall(MatSetOption(F, MAT_HERMITIAN, (PetscBool)(hpd || herm)));
120: PetscCall(MatSetOption(F, MAT_SPD, hpd));
121: {
122: PetscInt iftyp = ftyp;
123: PetscCall(PetscOptionsGetEList(NULL, NULL, "-ftype", MatFactorTypes, MAT_FACTOR_NUM_TYPES, &iftyp, NULL));
124: ftyp = (MatFactorType)iftyp;
125: }
126: if (ftyp == MAT_FACTOR_LU) {
127: PetscCall(MatLUFactor(F, NULL, NULL, NULL));
128: } else if (ftyp == MAT_FACTOR_CHOLESKY) {
129: PetscCall(MatCholeskyFactor(F, NULL, NULL));
130: } else if (ftyp == MAT_FACTOR_QR) {
131: PetscCall(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: PetscCall(VecSetRandom(x, rand));
136: PetscCall(VecCopy(x, u));
137: if (nsolve) {
138: PetscCall(MatMult(A, x, b));
139: PetscCall(MatSolve(F, b, x));
140: } else {
141: PetscCall(MatMultTranspose(A, x, b));
142: PetscCall(MatSolveTranspose(F, b, x));
143: }
144: /* Check the error */
145: PetscCall(VecAXPY(u, -1.0, x)); /* u <- (-1.0)x + u */
146: PetscCall(VecNorm(u, NORM_2, &norm));
147: if (norm > tol) {
148: PetscReal resi;
149: if (nsolve) {
150: PetscCall(MatMult(A, x, u)); /* u = A*x */
151: } else {
152: PetscCall(MatMultTranspose(A, x, u)); /* u = A*x */
153: }
154: PetscCall(VecAXPY(u, -1.0, b)); /* u <- (-1.0)b + u */
155: PetscCall(VecNorm(u, NORM_2, &resi));
156: if (nsolve) {
157: PetscCall(PetscPrintf(PETSC_COMM_SELF, "MatSolve error: Norm of error %g, residual %g\n", (double)norm, (double)resi));
158: } else {
159: PetscCall(PetscPrintf(PETSC_COMM_SELF, "MatSolveTranspose error: Norm of error %g, residual %g\n", (double)norm, (double)resi));
160: }
161: }
162: }
163: PetscCall(MatMatMult(A, C, MAT_REUSE_MATRIX, 2.0, &RHS));
164: PetscCall(MatMatSolve(F, RHS, X));
166: /* Check the error */
167: PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN));
168: PetscCall(MatNorm(X, NORM_FROBENIUS, &norm));
169: if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "MatMatSolve: Norm of error %g\n", (double)norm));
171: /* Free data structures */
172: PetscCall(MatDestroy(&A));
173: PetscCall(MatDestroy(&C));
174: PetscCall(MatDestroy(&F));
175: PetscCall(MatDestroy(&X));
176: PetscCall(MatDestroy(&RHS));
177: PetscCall(PetscRandomDestroy(&rand));
178: PetscCall(VecDestroy(&x));
179: PetscCall(VecDestroy(&b));
180: PetscCall(VecDestroy(&u));
181: PetscCall(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*/