1 const char help[] = "Test PCMatSetApplyOperation() and PCMatGetApplyOperation()";
2
3 #include <petscpc.h>
4
TestVecEquality(Vec x,Vec y)5 static PetscErrorCode TestVecEquality(Vec x, Vec y)
6 {
7 Vec diff;
8 PetscReal err, scale;
9
10 PetscFunctionBegin;
11 PetscCall(VecDuplicate(x, &diff));
12 PetscCall(VecCopy(x, diff));
13 PetscCall(VecAXPY(diff, -1.0, y));
14 PetscCall(VecNorm(diff, NORM_INFINITY, &err));
15 PetscCall(VecNorm(x, NORM_INFINITY, &scale));
16 PetscCheck(err <= PETSC_SMALL * scale, PetscObjectComm((PetscObject)x), PETSC_ERR_PLIB, "PC operation does not match Mat operation");
17 PetscCall(VecDestroy(&diff));
18 PetscFunctionReturn(PETSC_SUCCESS);
19 }
20
TestMatEquality(Mat x,Mat y)21 static PetscErrorCode TestMatEquality(Mat x, Mat y)
22 {
23 Mat diff;
24 PetscReal err, scale;
25 PetscInt m, n;
26
27 PetscFunctionBegin;
28 PetscCall(MatGetSize(x, &m, &n));
29 PetscCall(MatDuplicate(x, MAT_COPY_VALUES, &diff));
30 PetscCall(MatAXPY(diff, -1.0, y, SAME_NONZERO_PATTERN));
31 PetscCall(MatNorm(diff, NORM_FROBENIUS, &err));
32 PetscCall(MatNorm(x, NORM_FROBENIUS, &scale));
33 PetscCheck(err < PETSC_SMALL * m * n * scale, PetscObjectComm((PetscObject)x), PETSC_ERR_PLIB, "PC operation does not match Mat operation");
34 PetscCall(MatDestroy(&diff));
35 PetscFunctionReturn(PETSC_SUCCESS);
36 }
37
38 // Test PCMat with a given operation versus a Matrix that encode the same operation relative to the original matrix
TestPCMatVersusMat(PC pc,Mat A,PetscRandom rand,MatOperation op)39 static PetscErrorCode TestPCMatVersusMat(PC pc, Mat A, PetscRandom rand, MatOperation op)
40 {
41 Vec b, x, x2;
42 Mat B, X, X2;
43 MatOperation op2;
44
45 PetscFunctionBegin;
46 PetscCall(PCMatSetApplyOperation(pc, op));
47 PetscCall(PCView(pc, PETSC_VIEWER_STDOUT_SELF));
48 PetscCall(PCMatGetApplyOperation(pc, &op2));
49 PetscCheck(op == op2, PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "Input and output MatOperation differ");
50
51 PetscCall(MatCreateVecs(A, &b, &x));
52 PetscCall(VecDuplicate(x, &x2));
53 PetscCall(VecSetRandom(b, rand));
54
55 PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &B));
56 PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &X2));
57 PetscCall(MatSetRandom(B, rand));
58
59 PetscCall(MatMult(A, b, x));
60 PetscCall(PCApply(pc, b, x2));
61 PetscCall(TestVecEquality(x, x2));
62
63 PetscCall(MatMultTranspose(A, b, x));
64 PetscCall(PCApplyTranspose(pc, b, x2));
65 PetscCall(TestVecEquality(x, x2));
66
67 PetscCall(MatMatMult(A, B, MAT_INITIAL_MATRIX, PETSC_CURRENT, &X));
68 PetscCall(PCMatApply(pc, B, X2));
69 PetscCall(TestMatEquality(X, X2));
70
71 PetscCall(MatDestroy(&X2));
72 PetscCall(MatDestroy(&X));
73 PetscCall(MatDestroy(&B));
74 PetscCall(VecDestroy(&x2));
75 PetscCall(VecDestroy(&x));
76 PetscCall(VecDestroy(&b));
77 PetscFunctionReturn(PETSC_SUCCESS);
78 }
79
main(int argc,char ** argv)80 int main(int argc, char **argv)
81 {
82 PetscInt n = 10;
83 Mat A, AT, AH, II, Ainv, AinvT;
84 MPI_Comm comm;
85 PC pc;
86 PetscRandom rand;
87
88 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
89 comm = PETSC_COMM_SELF;
90
91 PetscCall(PetscRandomCreate(comm, &rand));
92 if (PetscDefined(USE_COMPLEX)) {
93 PetscScalar i = PetscSqrtScalar(-1.0);
94 PetscCall(PetscRandomSetInterval(rand, -1.0 - i, 1.0 + i));
95 } else {
96 PetscCall(PetscRandomSetInterval(rand, -1.0, 1.0));
97 }
98
99 PetscCall(MatCreateSeqDense(comm, n, n, NULL, &A));
100 PetscCall(MatSetRandom(A, rand));
101
102 PetscCall(PCCreate(comm, &pc));
103 PetscCall(PCSetType(pc, PCMAT));
104 PetscCall(PCSetOperators(pc, A, A));
105 PetscCall(PCSetUp(pc));
106
107 MatOperation default_op;
108 PetscCall(PCMatGetApplyOperation(pc, &default_op));
109 PetscCheck(default_op == MATOP_MULT, comm, PETSC_ERR_PLIB, "Default operation has changed");
110
111 // Test setting an invalid operation
112 PetscCall(PetscPushErrorHandler(PetscReturnErrorHandler, NULL));
113 PetscErrorCode ierr = PCMatSetApplyOperation(pc, MATOP_SET_VALUES);
114 PetscCall(PetscPopErrorHandler());
115 PetscCheck(ierr == PETSC_ERR_ARG_INCOMP, comm, PETSC_ERR_PLIB, "Wrong error message for unsupported MatOperation");
116
117 PetscCall(TestPCMatVersusMat(pc, A, rand, MATOP_MULT));
118
119 PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &AT));
120 PetscCall(TestPCMatVersusMat(pc, AT, rand, MATOP_MULT_TRANSPOSE));
121
122 PetscCall(MatHermitianTranspose(A, MAT_INITIAL_MATRIX, &AH));
123 PetscCall(TestPCMatVersusMat(pc, AH, rand, MATOP_MULT_HERMITIAN_TRANSPOSE));
124
125 PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &II));
126 PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &Ainv));
127 PetscCall(MatZeroEntries(II));
128 PetscCall(MatShift(II, 1.0));
129 PetscCall(MatLUFactor(A, NULL, NULL, NULL));
130 PetscCall(MatMatSolve(A, II, Ainv));
131 PetscCall(PCSetOperators(pc, A, A));
132 PetscCall(TestPCMatVersusMat(pc, Ainv, rand, MATOP_SOLVE));
133
134 PetscCall(MatTranspose(Ainv, MAT_INITIAL_MATRIX, &AinvT));
135 PetscCall(TestPCMatVersusMat(pc, AinvT, rand, MATOP_SOLVE_TRANSPOSE));
136
137 PetscCall(PCDestroy(&pc));
138 PetscCall(MatDestroy(&AinvT));
139 PetscCall(MatDestroy(&Ainv));
140 PetscCall(MatDestroy(&II));
141 PetscCall(MatDestroy(&AH));
142 PetscCall(MatDestroy(&AT));
143 PetscCall(MatDestroy(&A));
144 PetscCall(PetscRandomDestroy(&rand));
145 PetscCall(PetscFinalize());
146 return 0;
147 }
148
149 /*TEST
150
151 test:
152 suffix: 0
153
154 TEST*/
155