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