1 static char help[] = "Test different MatSolve routines with MATTRANSPOSEVIRTUAL.\n\n";
2
3 #include <petscmat.h>
4
TestMatrix(const char * test,Mat A,PetscInt nrhs,PetscBool inplace,PetscBool chol)5 PetscErrorCode TestMatrix(const char *test, Mat A, PetscInt nrhs, PetscBool inplace, PetscBool chol)
6 {
7 Mat F, RHS, X, C1;
8 Vec b, x, y, f;
9 IS perm, iperm;
10 PetscInt n, i;
11 PetscReal norm, tol = 1000 * PETSC_MACHINE_EPSILON;
12 PetscBool ht;
13 #if defined(PETSC_USE_COMPLEX)
14 PetscScalar v1 = PetscCMPLX(1.0, -0.1), v2 = PetscCMPLX(-1.0, 0.1);
15 #else
16 PetscScalar v1 = 1.0, v2 = -1.0;
17 #endif
18
19 PetscFunctionBegin;
20 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATHERMITIANTRANSPOSEVIRTUAL, &ht));
21 PetscCall(MatCreateVecs(A, &f, &b));
22 PetscCall(MatCreateVecs(A, &x, &y));
23 PetscCall(VecSet(b, v1));
24 PetscCall(VecSet(y, v2));
25
26 PetscCall(MatGetOrdering(A, MATORDERINGND, &perm, &iperm));
27 if (!inplace) {
28 if (!chol) {
29 PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_LU, &F));
30 PetscCall(MatLUFactorSymbolic(F, A, perm, iperm, NULL));
31 PetscCall(MatLUFactorNumeric(F, A, NULL));
32 } else { /* Cholesky */
33 PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &F));
34 PetscCall(MatCholeskyFactorSymbolic(F, A, perm, NULL));
35 PetscCall(MatCholeskyFactorNumeric(F, A, NULL));
36 }
37 } else { /* Test inplace factorization */
38 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &F));
39 if (!chol) PetscCall(MatLUFactor(F, perm, iperm, NULL));
40 else PetscCall(MatCholeskyFactor(F, perm, NULL));
41 }
42
43 /* MatSolve */
44 PetscCall(MatSolve(F, b, x));
45 PetscCall(MatMult(A, x, f));
46 PetscCall(VecAXPY(f, -1.0, b));
47 PetscCall(VecNorm(f, NORM_2, &norm));
48 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%12s MatSolve : Error of norm %g\n", test, (double)norm));
49
50 /* MatSolveTranspose */
51 if (!ht) {
52 PetscCall(MatSolveTranspose(F, b, x));
53 PetscCall(MatMultTranspose(A, x, f));
54 PetscCall(VecAXPY(f, -1.0, b));
55 PetscCall(VecNorm(f, NORM_2, &norm));
56 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%12s MatSolveTranspose : Error of norm %g\n", test, (double)norm));
57 }
58
59 /* MatSolveAdd */
60 PetscCall(MatSolveAdd(F, b, y, x));
61 PetscCall(MatMult(A, y, f));
62 PetscCall(VecScale(f, -1.0));
63 PetscCall(MatMultAdd(A, x, f, f));
64 PetscCall(VecAXPY(f, -1.0, b));
65 PetscCall(VecNorm(f, NORM_2, &norm));
66 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%12s MatSolveAdd : Error of norm %g\n", test, (double)norm));
67
68 /* MatSolveTransposeAdd */
69 if (!ht) {
70 PetscCall(MatSolveTransposeAdd(F, b, y, x));
71 PetscCall(MatMultTranspose(A, y, f));
72 PetscCall(VecScale(f, -1.0));
73 PetscCall(MatMultTransposeAdd(A, x, f, f));
74 PetscCall(VecAXPY(f, -1.0, b));
75 PetscCall(VecNorm(f, NORM_2, &norm));
76 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%12s MatSolveTransposeAdd : Error of norm %g\n", test, (double)norm));
77 }
78
79 /* MatMatSolve */
80 PetscCall(MatGetSize(A, &n, NULL));
81 PetscCall(MatCreate(PETSC_COMM_WORLD, &RHS));
82 PetscCall(MatSetSizes(RHS, PETSC_DECIDE, PETSC_DECIDE, n, nrhs));
83 PetscCall(MatSetType(RHS, MATSEQDENSE));
84 PetscCall(MatSetUp(RHS));
85 for (i = 0; i < nrhs; i++) PetscCall(MatSetValue(RHS, i, i, 1.0, INSERT_VALUES));
86 PetscCall(MatAssemblyBegin(RHS, MAT_FINAL_ASSEMBLY));
87 PetscCall(MatAssemblyEnd(RHS, MAT_FINAL_ASSEMBLY));
88 PetscCall(MatDuplicate(RHS, MAT_DO_NOT_COPY_VALUES, &X));
89
90 if (!ht) {
91 PetscCall(MatMatSolve(F, RHS, X));
92 PetscCall(MatMatMult(A, X, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &C1));
93 PetscCall(MatAXPY(C1, -1.0, RHS, SAME_NONZERO_PATTERN));
94 PetscCall(MatNorm(C1, NORM_FROBENIUS, &norm));
95 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%12s MatMatSolve : Error of norm %g\n", test, (double)norm));
96 PetscCall(MatDestroy(&C1));
97
98 PetscCall(MatMatSolveTranspose(F, RHS, X));
99 PetscCall(MatTransposeMatMult(A, X, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &C1));
100 PetscCall(MatAXPY(C1, -1.0, RHS, SAME_NONZERO_PATTERN));
101 PetscCall(MatNorm(C1, NORM_FROBENIUS, &norm));
102 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%12s MatMatSolveTranspose : Error of norm %g\n", test, (double)norm));
103 PetscCall(MatDestroy(&C1));
104 }
105 PetscCall(VecDestroy(&b));
106 PetscCall(VecDestroy(&x));
107 PetscCall(VecDestroy(&f));
108 PetscCall(VecDestroy(&y));
109 PetscCall(ISDestroy(&perm));
110 PetscCall(ISDestroy(&iperm));
111 PetscCall(MatDestroy(&F));
112 PetscCall(MatDestroy(&RHS));
113 PetscCall(MatDestroy(&X));
114 PetscFunctionReturn(PETSC_SUCCESS);
115 }
116
main(int argc,char ** args)117 int main(int argc, char **args)
118 {
119 PetscMPIInt size;
120 Mat A, At, Aht;
121 PetscInt i, n = 8, nrhs = 2;
122 PetscBool aij, inplace = PETSC_FALSE;
123 #if defined(PETSC_USE_COMPLEX)
124 PetscScalar a = PetscCMPLX(-1.0, 0.5);
125 #else
126 PetscScalar a = -1.0;
127 #endif
128
129 PetscFunctionBeginUser;
130 PetscCall(PetscInitialize(&argc, &args, NULL, help));
131 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
132 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only");
133 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
134 PetscCall(PetscOptionsGetInt(NULL, NULL, "-nrhs", &nrhs, NULL));
135 PetscCall(PetscOptionsGetBool(NULL, NULL, "-inplace", &inplace, NULL));
136 PetscCheck(nrhs <= n, PETSC_COMM_WORLD, PETSC_ERR_ARG_SIZ, "Must have nrhs <= n");
137
138 /* Hermitian matrix */
139 PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
140 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n));
141 PetscCall(MatSetFromOptions(A));
142 for (i = 0; i < n; i++) {
143 if (i > 0) PetscCall(MatSetValue(A, i, i - 1, a, INSERT_VALUES));
144 if (i < n - 1) PetscCall(MatSetValue(A, i, i + 1, PetscConj(a), INSERT_VALUES));
145 PetscCall(MatSetValue(A, i, i, 2.0, INSERT_VALUES));
146 }
147 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
148 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
149
150 PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &aij, MATSEQAIJ, MATSEQBAIJ, ""));
151 PetscCall(MatSetOption(A, PetscDefined(USE_COMPLEX) ? MAT_HERMITIAN : MAT_SYMMETRIC, PETSC_TRUE));
152
153 PetscCall(MatCreateTranspose(A, &At));
154 PetscCall(MatCreateHermitianTranspose(A, &Aht));
155
156 PetscCall(TestMatrix("LU T", At, nrhs, inplace, PETSC_FALSE));
157 PetscCall(TestMatrix("LU HT", Aht, nrhs, inplace, PETSC_FALSE));
158 if (!aij) {
159 PetscCall(TestMatrix("Chol T", At, nrhs, inplace, PETSC_TRUE));
160 PetscCall(TestMatrix("Chol HT", Aht, nrhs, inplace, PETSC_TRUE));
161 }
162
163 /* Make the matrix non-Hermitian */
164 PetscCall(MatSetValue(A, 0, 1, -5.0, INSERT_VALUES));
165 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
166 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
167 #if defined(PETSC_USE_COMPLEX)
168 PetscCall(MatSetOption(A, MAT_HERMITIAN, PETSC_FALSE));
169 #else
170 PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_FALSE));
171 #endif
172
173 PetscCall(TestMatrix("LU T nonsym", At, nrhs, inplace, PETSC_FALSE));
174 PetscCall(TestMatrix("LU HT nonsym", Aht, nrhs, inplace, PETSC_FALSE));
175
176 PetscCall(MatDestroy(&A));
177 PetscCall(MatDestroy(&At));
178 PetscCall(MatDestroy(&Aht));
179 PetscCall(PetscFinalize());
180 return 0;
181 }
182
183 /*TEST
184
185 test:
186 suffix: 1
187 args: -inplace {{0 1}} -mat_type {{aij dense}}
188 output_file: output/empty.out
189
190 TEST*/
191