1 static char help[] = "Tests MatShift(), MatScale(), and MatDiagonalScale() for SHELL and NEST matrices\n\n"; 2 3 #include <petscmat.h> 4 5 typedef struct _n_User *User; 6 struct _n_User { 7 Mat B; 8 }; 9 10 static PetscErrorCode MatView_User(Mat A, PetscViewer viewer) 11 { 12 User user; 13 14 PetscFunctionBegin; 15 PetscCall(MatShellGetContext(A, &user)); 16 PetscCall(MatView(user->B, viewer)); 17 PetscFunctionReturn(PETSC_SUCCESS); 18 } 19 20 static PetscErrorCode MatMult_User(Mat A, Vec X, Vec Y) 21 { 22 User user; 23 24 PetscFunctionBegin; 25 PetscCall(MatShellGetContext(A, &user)); 26 PetscCall(MatMult(user->B, X, Y)); 27 PetscFunctionReturn(PETSC_SUCCESS); 28 } 29 30 static PetscErrorCode MatMultTranspose_User(Mat A, Vec X, Vec Y) 31 { 32 User user; 33 34 PetscFunctionBegin; 35 PetscCall(MatShellGetContext(A, &user)); 36 PetscCall(MatMultTranspose(user->B, X, Y)); 37 PetscFunctionReturn(PETSC_SUCCESS); 38 } 39 40 static PetscErrorCode MatGetDiagonal_User(Mat A, Vec X) 41 { 42 User user; 43 44 PetscFunctionBegin; 45 PetscCall(MatShellGetContext(A, &user)); 46 PetscCall(MatGetDiagonal(user->B, X)); 47 PetscFunctionReturn(PETSC_SUCCESS); 48 } 49 50 static PetscErrorCode TestMatrix(Mat A, Vec X, Vec Y, Vec Z) 51 { 52 Vec W1, W2, diff; 53 Mat E; 54 const char *mattypename; 55 PetscViewer viewer = PETSC_VIEWER_STDOUT_WORLD; 56 PetscScalar diag[2] = {2.9678190300000000e+08, 1.4173141580000000e+09}; 57 PetscScalar multadd[2] = {-6.8966198500000000e+08, -2.0310609940000000e+09}; 58 PetscScalar multtadd[2] = {-9.1052873900000000e+08, -1.8101942400000000e+09}; 59 PetscReal nrm; 60 61 PetscFunctionBegin; 62 PetscCall(PetscObjectGetType((PetscObject)A, &mattypename)); 63 PetscCall(PetscViewerASCIIPrintf(viewer, "\nMatrix of type: %s\n", mattypename)); 64 PetscCall(VecDuplicate(X, &W1)); 65 PetscCall(VecDuplicate(X, &W2)); 66 PetscCall(MatScale(A, 31)); 67 PetscCall(MatShift(A, 37)); 68 PetscCall(MatDiagonalScale(A, X, Y)); 69 PetscCall(MatScale(A, 41)); 70 PetscCall(MatDiagonalScale(A, Y, Z)); 71 PetscCall(MatComputeOperator(A, MATDENSE, &E)); 72 73 PetscCall(MatView(E, viewer)); 74 PetscCall(PetscViewerASCIIPrintf(viewer, "Testing MatMult + MatMultTranspose\n")); 75 PetscCall(MatMult(A, Z, W1)); 76 PetscCall(MatMultTranspose(A, W1, W2)); 77 PetscCall(VecView(W2, viewer)); 78 79 PetscCall(PetscViewerASCIIPrintf(viewer, "Testing MatMultAdd\n")); 80 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, 2, multadd, &diff)); 81 PetscCall(VecSet(W1, -1.0)); 82 PetscCall(MatMultAdd(A, W1, W1, W2)); 83 PetscCall(VecView(W2, viewer)); 84 PetscCall(VecAXPY(W2, -1.0, diff)); 85 PetscCall(VecNorm(W2, NORM_2, &nrm)); 86 #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128) 87 PetscCheck(nrm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatMultAdd(A,x,x,y) produces incorrect result"); 88 #endif 89 90 PetscCall(VecSet(W2, -1.0)); 91 PetscCall(MatMultAdd(A, W1, W2, W2)); 92 PetscCall(VecView(W2, viewer)); 93 PetscCall(VecAXPY(W2, -1.0, diff)); 94 PetscCall(VecNorm(W2, NORM_2, &nrm)); 95 #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128) 96 PetscCheck(nrm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatMultAdd(A,x,y,y) produces incorrect result"); 97 #endif 98 PetscCall(VecDestroy(&diff)); 99 100 PetscCall(PetscViewerASCIIPrintf(viewer, "Testing MatMultTransposeAdd\n")); 101 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, 2, multtadd, &diff)); 102 103 PetscCall(VecSet(W1, -1.0)); 104 PetscCall(MatMultTransposeAdd(A, W1, W1, W2)); 105 PetscCall(VecView(W2, viewer)); 106 PetscCall(VecAXPY(W2, -1.0, diff)); 107 PetscCall(VecNorm(W2, NORM_2, &nrm)); 108 #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128) 109 PetscCheck(nrm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatMultTransposeAdd(A,x,x,y) produces incorrect result"); 110 #endif 111 112 PetscCall(VecSet(W2, -1.0)); 113 PetscCall(MatMultTransposeAdd(A, W1, W2, W2)); 114 PetscCall(VecView(W2, viewer)); 115 PetscCall(VecAXPY(W2, -1.0, diff)); 116 PetscCall(VecNorm(W2, NORM_2, &nrm)); 117 #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128) 118 PetscCheck(nrm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatMultTransposeAdd(A,x,y,y) produces incorrect result"); 119 #endif 120 PetscCall(VecDestroy(&diff)); 121 122 PetscCall(PetscViewerASCIIPrintf(viewer, "Testing MatGetDiagonal\n")); 123 PetscCall(MatGetDiagonal(A, W2)); 124 PetscCall(VecView(W2, viewer)); 125 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, 2, diag, &diff)); 126 PetscCall(VecAXPY(diff, -1.0, W2)); 127 PetscCall(VecNorm(diff, NORM_2, &nrm)); 128 PetscCheck(nrm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatGetDiagonal() produces incorrect result"); 129 PetscCall(VecDestroy(&diff)); 130 131 /* MATSHELL does not support MatDiagonalSet after MatScale */ 132 if (strncmp(mattypename, "shell", 5)) { 133 PetscCall(MatDiagonalSet(A, X, INSERT_VALUES)); 134 PetscCall(MatGetDiagonal(A, W1)); 135 PetscCall(VecView(W1, viewer)); 136 } else { 137 PetscCall(PetscViewerASCIIPrintf(viewer, "MatDiagonalSet not tested on MATSHELL\n")); 138 } 139 140 PetscCall(MatDestroy(&E)); 141 PetscCall(VecDestroy(&W1)); 142 PetscCall(VecDestroy(&W2)); 143 PetscFunctionReturn(PETSC_SUCCESS); 144 } 145 146 int main(int argc, char **args) 147 { 148 const PetscScalar xvals[] = {11, 13}, yvals[] = {17, 19}, zvals[] = {23, 29}; 149 const PetscInt inds[] = {0, 1}; 150 PetscScalar avals[] = {2, 3, 5, 7}; 151 Mat A, S, D[4], N; 152 Vec X, Y, Z; 153 User user; 154 PetscInt i; 155 156 PetscFunctionBeginUser; 157 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 158 PetscCall(MatCreateSeqAIJ(PETSC_COMM_WORLD, 2, 2, 2, NULL, &A)); 159 PetscCall(MatSetUp(A)); 160 PetscCall(MatSetValues(A, 2, inds, 2, inds, avals, INSERT_VALUES)); 161 PetscCall(VecCreateSeq(PETSC_COMM_WORLD, 2, &X)); 162 PetscCall(VecDuplicate(X, &Y)); 163 PetscCall(VecDuplicate(X, &Z)); 164 PetscCall(VecSetValues(X, 2, inds, xvals, INSERT_VALUES)); 165 PetscCall(VecSetValues(Y, 2, inds, yvals, INSERT_VALUES)); 166 PetscCall(VecSetValues(Z, 2, inds, zvals, INSERT_VALUES)); 167 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 168 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 169 PetscCall(VecAssemblyBegin(X)); 170 PetscCall(VecAssemblyBegin(Y)); 171 PetscCall(VecAssemblyBegin(Z)); 172 PetscCall(VecAssemblyEnd(X)); 173 PetscCall(VecAssemblyEnd(Y)); 174 PetscCall(VecAssemblyEnd(Z)); 175 176 PetscCall(PetscNew(&user)); 177 user->B = A; 178 179 PetscCall(MatCreateShell(PETSC_COMM_WORLD, 2, 2, 2, 2, user, &S)); 180 PetscCall(MatSetUp(S)); 181 PetscCall(MatShellSetOperation(S, MATOP_VIEW, (void (*)(void))MatView_User)); 182 PetscCall(MatShellSetOperation(S, MATOP_MULT, (void (*)(void))MatMult_User)); 183 PetscCall(MatShellSetOperation(S, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMultTranspose_User)); 184 PetscCall(MatShellSetOperation(S, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_User)); 185 186 for (i = 0; i < 4; i++) PetscCall(MatCreateSeqDense(PETSC_COMM_WORLD, 1, 1, &avals[i], &D[i])); 187 PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, D, &N)); 188 PetscCall(MatSetUp(N)); 189 190 PetscCall(TestMatrix(S, X, Y, Z)); 191 PetscCall(TestMatrix(A, X, Y, Z)); 192 PetscCall(TestMatrix(N, X, Y, Z)); 193 194 for (i = 0; i < 4; i++) PetscCall(MatDestroy(&D[i])); 195 PetscCall(MatDestroy(&A)); 196 PetscCall(MatDestroy(&S)); 197 PetscCall(MatDestroy(&N)); 198 PetscCall(VecDestroy(&X)); 199 PetscCall(VecDestroy(&Y)); 200 PetscCall(VecDestroy(&Z)); 201 PetscCall(PetscFree(user)); 202 PetscCall(PetscFinalize()); 203 return 0; 204 } 205 206 /*TEST 207 208 test: 209 210 TEST*/ 211