1 static char help[] = "Tests MatSolve(), MatSolveTranspose() and MatMatSolve() with SEQDENSE\n";
2
3 #include <petscmat.h>
4
main(int argc,char ** args)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];
17
18 PetscFunctionBeginUser;
19 PetscCall(PetscInitialize(&argc, &args, NULL, 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));
29
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 }
60
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);
72
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));
81
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));
87
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 */
94
95 /* make a symmetric matrix */
96 if (symm) {
97 Mat AT;
98
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;
107
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"));
115
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]);
133
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));
165
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));
170
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 }
184
185 /*TEST
186
187 testset:
188 output_file: output/empty.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
203
204 TEST*/
205