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_DEFAULT, &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*/