Actual source code: ex267.c
1: static char help[] = "Test different MatSolve routines with MATTRANSPOSEVIRTUAL.\n\n";
3: #include <petscmat.h>
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
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));
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: }
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));
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: }
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));
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: }
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));
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));
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: }
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
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");
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));
150: PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &aij, MATSEQAIJ, MATSEQBAIJ, ""));
151: #if defined(PETSC_USE_COMPLEX)
152: PetscCall(MatSetOption(A, MAT_HERMITIAN, PETSC_TRUE));
153: #else
154: PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE));
155: #endif
157: PetscCall(MatCreateTranspose(A, &At));
158: PetscCall(MatCreateHermitianTranspose(A, &Aht));
160: PetscCall(TestMatrix("LU T", At, nrhs, inplace, PETSC_FALSE));
161: PetscCall(TestMatrix("LU HT", Aht, nrhs, inplace, PETSC_FALSE));
162: if (!aij) {
163: PetscCall(TestMatrix("Chol T", At, nrhs, inplace, PETSC_TRUE));
164: PetscCall(TestMatrix("Chol HT", Aht, nrhs, inplace, PETSC_TRUE));
165: }
167: /* Make the matrix non-Hermitian */
168: PetscCall(MatSetValue(A, 0, 1, -5.0, INSERT_VALUES));
169: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
170: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
171: #if defined(PETSC_USE_COMPLEX)
172: PetscCall(MatSetOption(A, MAT_HERMITIAN, PETSC_FALSE));
173: #else
174: PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_FALSE));
175: #endif
177: PetscCall(TestMatrix("LU T nonsym", At, nrhs, inplace, PETSC_FALSE));
178: PetscCall(TestMatrix("LU HT nonsym", Aht, nrhs, inplace, PETSC_FALSE));
180: PetscCall(MatDestroy(&A));
181: PetscCall(MatDestroy(&At));
182: PetscCall(MatDestroy(&Aht));
183: PetscCall(PetscFinalize());
184: return 0;
185: }
187: /*TEST
189: test:
190: suffix: 1
191: args: -inplace {{0 1}} -mat_type {{aij dense}}
192: output_file: output/empty.out
194: TEST*/