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, 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));
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 (PetscDefined(USE_COMPLEX)) PetscCall(PetscStrlcat(file, "hpd-complex-", sizeof(file)));
37: else PetscCall(PetscStrlcat(file, "spd-real-", sizeof(file)));
38: ftyp = MAT_FACTOR_CHOLESKY;
39: } else {
40: if (PetscDefined(USE_COMPLEX)) PetscCall(PetscStrlcat(file, "nh-complex-", sizeof(file)));
41: else PetscCall(PetscStrlcat(file, "ns-real-", sizeof(file)));
42: }
43: if (PetscDefined(USE_64BIT_INDICES)) PetscCall(PetscStrlcat(file, "int64-", sizeof(file)));
44: else PetscCall(PetscStrlcat(file, "int32-", sizeof(file)));
45: if (PetscDefined(USE_REAL_SINGLE)) PetscCall(PetscStrlcat(file, "float32", sizeof(file)));
46: else PetscCall(PetscStrlcat(file, "float64", sizeof(file)));
47: }
49: /* Load matrix A */
50: if (PetscDefined(USE_REAL___FLOAT128)) PetscCall(PetscOptionsInsertString(NULL, "-binary_read_double"));
51: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd));
52: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
53: PetscCall(MatLoad(A, fd));
54: PetscCall(PetscViewerDestroy(&fd));
55: PetscCall(MatConvert(A, MATSEQDENSE, MAT_INPLACE_MATRIX, &A));
56: PetscCall(MatGetSize(A, &m, &n));
57: PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")", m, n);
59: /* Create dense matrix C and X; C holds true solution with identical columns */
60: nrhs = 2;
61: PetscCall(PetscOptionsGetInt(NULL, NULL, "-nrhs", &nrhs, NULL));
62: PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
63: PetscCall(MatSetSizes(C, m, PETSC_DECIDE, PETSC_DECIDE, nrhs));
64: PetscCall(MatSetType(C, MATDENSE));
65: PetscCall(MatSetFromOptions(C));
66: PetscCall(MatSetUp(C));
68: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand));
69: PetscCall(PetscRandomSetFromOptions(rand));
70: PetscCall(MatSetRandom(C, rand));
71: PetscCall(MatDuplicate(C, MAT_DO_NOT_COPY_VALUES, &X));
72: PetscCall(MatDuplicate(C, MAT_DO_NOT_COPY_VALUES, &RHS));
74: /* Create vectors */
75: PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
76: PetscCall(VecSetSizes(x, n, PETSC_DECIDE));
77: PetscCall(VecSetFromOptions(x));
78: PetscCall(VecDuplicate(x, &b));
79: PetscCall(VecDuplicate(x, &u)); /* save the true solution */
81: /* make a symmetric matrix */
82: if (symm) {
83: Mat AT;
85: PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &AT));
86: PetscCall(MatAXPY(A, 1.0, AT, SAME_NONZERO_PATTERN));
87: PetscCall(MatDestroy(&AT));
88: ftyp = MAT_FACTOR_CHOLESKY;
89: }
90: /* make an hermitian matrix */
91: if (herm) {
92: Mat AH;
94: PetscCall(MatHermitianTranspose(A, MAT_INITIAL_MATRIX, &AH));
95: PetscCall(MatAXPY(A, 1.0, AH, SAME_NONZERO_PATTERN));
96: PetscCall(MatDestroy(&AH));
97: ftyp = MAT_FACTOR_CHOLESKY;
98: }
99: PetscCall(PetscObjectSetName((PetscObject)A, "A"));
100: PetscCall(MatViewFromOptions(A, NULL, "-amat_view"));
102: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &F));
103: PetscCall(MatSetOption(F, MAT_SYMMETRIC, symm));
104: /* it seems that the SPD concept in PETSc extends naturally to Hermitian Positive definitess */
105: PetscCall(MatSetOption(F, MAT_HERMITIAN, (PetscBool)(hpd || herm)));
106: PetscCall(MatSetOption(F, MAT_SPD, hpd));
107: {
108: PetscInt iftyp = ftyp;
109: PetscCall(PetscOptionsGetEList(NULL, NULL, "-ftype", MatFactorTypes, MAT_FACTOR_NUM_TYPES, &iftyp, NULL));
110: ftyp = (MatFactorType)iftyp;
111: }
112: if (ftyp == MAT_FACTOR_LU) {
113: PetscCall(MatLUFactor(F, NULL, NULL, NULL));
114: } else if (ftyp == MAT_FACTOR_CHOLESKY) {
115: PetscCall(MatCholeskyFactor(F, NULL, NULL));
116: } else if (ftyp == MAT_FACTOR_QR) {
117: PetscCall(MatQRFactor(F, NULL, NULL));
118: } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Factorization %s not supported in this example", MatFactorTypes[ftyp]);
120: for (nsolve = 0; nsolve < 2; nsolve++) {
121: PetscCall(VecSetRandom(x, rand));
122: PetscCall(VecCopy(x, u));
123: if (nsolve) {
124: PetscCall(MatMult(A, x, b));
125: PetscCall(MatSolve(F, b, x));
126: } else {
127: PetscCall(MatMultTranspose(A, x, b));
128: PetscCall(MatSolveTranspose(F, b, x));
129: }
130: /* Check the error */
131: PetscCall(VecAXPY(u, -1.0, x)); /* u <- (-1.0)x + u */
132: PetscCall(VecNorm(u, NORM_2, &norm));
133: if (norm > tol) {
134: PetscReal resi;
135: if (nsolve) {
136: PetscCall(MatMult(A, x, u)); /* u = A*x */
137: } else {
138: PetscCall(MatMultTranspose(A, x, u)); /* u = A*x */
139: }
140: PetscCall(VecAXPY(u, -1.0, b)); /* u <- (-1.0)b + u */
141: PetscCall(VecNorm(u, NORM_2, &resi));
142: if (nsolve) {
143: PetscCall(PetscPrintf(PETSC_COMM_SELF, "MatSolve error: Norm of error %g, residual %g\n", (double)norm, (double)resi));
144: } else {
145: PetscCall(PetscPrintf(PETSC_COMM_SELF, "MatSolveTranspose error: Norm of error %g, residual %g\n", (double)norm, (double)resi));
146: }
147: }
148: }
149: PetscCall(MatMatMult(A, C, MAT_REUSE_MATRIX, 2.0, &RHS));
150: PetscCall(MatMatSolve(F, RHS, X));
152: /* Check the error */
153: PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN));
154: PetscCall(MatNorm(X, NORM_FROBENIUS, &norm));
155: if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "MatMatSolve: Norm of error %g\n", (double)norm));
157: /* Free data structures */
158: PetscCall(MatDestroy(&A));
159: PetscCall(MatDestroy(&C));
160: PetscCall(MatDestroy(&F));
161: PetscCall(MatDestroy(&X));
162: PetscCall(MatDestroy(&RHS));
163: PetscCall(PetscRandomDestroy(&rand));
164: PetscCall(VecDestroy(&x));
165: PetscCall(VecDestroy(&b));
166: PetscCall(VecDestroy(&u));
167: PetscCall(PetscFinalize());
168: return 0;
169: }
171: /*TEST
173: testset:
174: output_file: output/empty.out
175: test:
176: suffix: ns
177: test:
178: suffix: sym
179: args: -symmetric_solve
180: test:
181: suffix: herm
182: args: -hermitian_solve
183: test:
184: suffix: hpd
185: args: -hpd_solve
186: test:
187: suffix: qr
188: args: -ftype qr
190: TEST*/