1 const char help[] = "Test PCMatSetApplyOperation() and PCMatGetApplyOperation()"; 2 3 #include <petscpc.h> 4 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 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 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_DEFAULT, &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 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