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