Actual source code: ex1.c
1: const char help[] = "Test PCMatSetApplyOperation() and PCMatGetApplyOperation()";
3: #include <petscpc.h>
5: static PetscErrorCode TestVecEquality(Vec x, Vec y)
6: {
7: Vec diff;
8: PetscReal err, scale;
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: }
21: static PetscErrorCode TestMatEquality(Mat x, Mat y)
22: {
23: Mat diff;
24: PetscReal err, scale;
25: PetscInt m, n;
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: }
38: // Test PCMat with a given operation versus a Matrix that encode the same operation relative to the original matrix
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;
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");
51: PetscCall(MatCreateVecs(A, &b, &x));
52: PetscCall(VecDuplicate(x, &x2));
53: PetscCall(VecSetRandom(b, rand));
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));
59: PetscCall(MatMult(A, b, x));
60: PetscCall(PCApply(pc, b, x2));
61: PetscCall(TestVecEquality(x, x2));
63: PetscCall(MatMultTranspose(A, b, x));
64: PetscCall(PCApplyTranspose(pc, b, x2));
65: PetscCall(TestVecEquality(x, x2));
67: PetscCall(MatMatMult(A, B, MAT_INITIAL_MATRIX, PETSC_CURRENT, &X));
68: PetscCall(PCMatApply(pc, B, X2));
69: PetscCall(TestMatEquality(X, X2));
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: }
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;
88: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
89: comm = PETSC_COMM_SELF;
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: }
99: PetscCall(MatCreateSeqDense(comm, n, n, NULL, &A));
100: PetscCall(MatSetRandom(A, rand));
102: PetscCall(PCCreate(comm, &pc));
103: PetscCall(PCSetType(pc, PCMAT));
104: PetscCall(PCSetOperators(pc, A, A));
105: PetscCall(PCSetUp(pc));
107: MatOperation default_op;
108: PetscCall(PCMatGetApplyOperation(pc, &default_op));
109: PetscCheck(default_op == MATOP_MULT, comm, PETSC_ERR_PLIB, "Default operation has changed");
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");
117: PetscCall(TestPCMatVersusMat(pc, A, rand, MATOP_MULT));
119: PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &AT));
120: PetscCall(TestPCMatVersusMat(pc, AT, rand, MATOP_MULT_TRANSPOSE));
122: PetscCall(MatHermitianTranspose(A, MAT_INITIAL_MATRIX, &AH));
123: PetscCall(TestPCMatVersusMat(pc, AH, rand, MATOP_MULT_HERMITIAN_TRANSPOSE));
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));
134: PetscCall(MatTranspose(Ainv, MAT_INITIAL_MATRIX, &AinvT));
135: PetscCall(TestPCMatVersusMat(pc, AinvT, rand, MATOP_SOLVE_TRANSPOSE));
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: }
149: /*TEST
151: test:
152: suffix: 0
154: TEST*/