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, W3, diff; 53 Mat E; 54 const char *mattypename; 55 PetscViewer viewer = PETSC_VIEWER_STDOUT_WORLD; 56 PetscReal nrm; 57 #if defined(PETSC_USE_COMPLEX) 58 const PetscScalar diag[2] = {PetscCMPLX(-6.2902938000000000e+07, 4.5741953400000000e+08), PetscCMPLX(1.0828994620000000e+09, 1.2955916360000000e+09)}; 59 const PetscScalar multadd[2] = {PetscCMPLX(1.4926230300000000e+08, -1.2811063360000000e+09), PetscCMPLX(-1.2985220710000000e+09, -2.1029893020000000e+09)}; 60 const PetscScalar multtadd[2] = {PetscCMPLX(-1.5271967100000000e+08, -1.2648172000000000e+09), PetscCMPLX(-9.9654009700000000e+08, -2.1192784380000000e+09)}; 61 #else 62 const PetscScalar diag[2] = {2.9678190300000000e+08, 1.4173141580000000e+09}; 63 const PetscScalar multadd[2] = {-6.8966198500000000e+08, -2.0310609940000000e+09}; 64 const PetscScalar multtadd[2] = {-9.1052873900000000e+08, -1.8101942400000000e+09}; 65 #endif 66 67 PetscFunctionBegin; 68 PetscCall(PetscObjectGetType((PetscObject)A, &mattypename)); 69 PetscCall(PetscViewerASCIIPrintf(viewer, "\nMatrix of type: %s\n", mattypename)); 70 PetscCall(VecDuplicate(X, &W1)); 71 PetscCall(VecDuplicate(X, &W2)); 72 PetscCall(VecDuplicate(X, &W3)); 73 PetscCall(MatScale(A, 31)); 74 PetscCall(MatShift(A, 37)); 75 PetscCall(MatDiagonalScale(A, X, Y)); 76 PetscCall(MatScale(A, 41)); 77 PetscCall(MatDiagonalScale(A, Y, Z)); 78 PetscCall(MatComputeOperator(A, MATDENSE, &E)); 79 80 PetscCall(MatView(E, viewer)); 81 PetscCall(PetscViewerASCIIPrintf(viewer, "Testing MatMult + MatMultTranspose\n")); 82 PetscCall(MatMult(A, Z, W1)); 83 PetscCall(MatMultTranspose(A, W1, W2)); 84 PetscCall(VecView(W2, viewer)); 85 PetscCall(PetscViewerASCIIPrintf(viewer, "Testing MatMultHermitianTranspose\n")); 86 PetscCall(VecConjugate(W1)); 87 PetscCall(MatMultHermitianTranspose(A, W1, W2)); 88 PetscCall(VecConjugate(W2)); 89 PetscCall(VecView(W2, viewer)); 90 91 PetscCall(PetscViewerASCIIPrintf(viewer, "Testing MatMultAdd\n")); 92 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, 2, multadd, &diff)); 93 PetscCall(VecSet(W1, -1.0)); 94 PetscCall(MatMultAdd(A, W1, W1, W2)); 95 PetscCall(VecView(W2, viewer)); 96 PetscCall(VecAXPY(W2, -1.0, diff)); 97 PetscCall(VecNorm(W2, NORM_2, &nrm)); 98 #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128) 99 PetscCheck(nrm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatMultAdd(A,x,x,y) produces incorrect result"); 100 #endif 101 102 PetscCall(VecSet(W2, -1.0)); 103 PetscCall(MatMultAdd(A, W1, W2, W2)); 104 PetscCall(VecView(W2, viewer)); 105 PetscCall(VecAXPY(W2, -1.0, diff)); 106 PetscCall(VecNorm(W2, NORM_2, &nrm)); 107 #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128) 108 PetscCheck(nrm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatMultAdd(A,x,y,y) produces incorrect result"); 109 #endif 110 PetscCall(VecDestroy(&diff)); 111 112 PetscCall(PetscViewerASCIIPrintf(viewer, "Testing MatMultTransposeAdd\n")); 113 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, 2, multtadd, &diff)); 114 115 PetscCall(VecSet(W1, -1.0)); 116 PetscCall(MatMultTransposeAdd(A, W1, W1, W2)); 117 PetscCall(VecView(W2, viewer)); 118 PetscCall(VecAXPY(W2, -1.0, diff)); 119 PetscCall(VecNorm(W2, NORM_2, &nrm)); 120 #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128) 121 PetscCheck(nrm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatMultTransposeAdd(A,x,x,y) produces incorrect result"); 122 #endif 123 124 PetscCall(VecSet(W2, -1.0)); 125 PetscCall(MatMultTransposeAdd(A, W1, W2, W2)); 126 PetscCall(VecView(W2, viewer)); 127 PetscCall(VecAXPY(W2, -1.0, diff)); 128 PetscCall(VecNorm(W2, NORM_2, &nrm)); 129 #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128) 130 PetscCheck(nrm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatMultTransposeAdd(A,x,y,y) produces incorrect result"); 131 #endif 132 PetscCall(VecDestroy(&diff)); 133 134 PetscCall(PetscViewerASCIIPrintf(viewer, "Testing MatMultHermitianTransposeAdd\n")); 135 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, 2, multtadd, &diff)); 136 137 PetscCall(VecSet(W1, -1.0)); 138 PetscCall(MatMultHermitianTransposeAdd(A, W1, W1, W3)); 139 PetscCall(VecConjugate(W3)); 140 PetscCall(VecView(W3, viewer)); 141 PetscCall(VecAXPY(W3, -1.0, diff)); 142 PetscCall(VecNorm(W3, NORM_2, &nrm)); 143 #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128) 144 PetscCheck(nrm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatMultHermitianTransposeAdd(A,x,x,y) produces incorrect result"); 145 #endif 146 147 PetscCall(VecSet(W3, -1.0)); 148 PetscCall(MatMultHermitianTransposeAdd(A, W1, W3, W3)); 149 PetscCall(VecConjugate(W3)); 150 PetscCall(VecView(W3, viewer)); 151 PetscCall(VecAXPY(W3, -1.0, diff)); 152 PetscCall(VecNorm(W3, NORM_2, &nrm)); 153 #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128) 154 PetscCheck(nrm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatMultHermitianTransposeAdd(A,x,y,y) produces incorrect result"); 155 #endif 156 PetscCall(VecDestroy(&diff)); 157 158 PetscCall(PetscViewerASCIIPrintf(viewer, "Testing MatGetDiagonal\n")); 159 PetscCall(MatGetDiagonal(A, W2)); 160 PetscCall(VecView(W2, viewer)); 161 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, 2, diag, &diff)); 162 PetscCall(VecAXPY(diff, -1.0, W2)); 163 PetscCall(VecNorm(diff, NORM_2, &nrm)); 164 #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128) 165 PetscCheck(nrm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatGetDiagonal() produces incorrect result"); 166 #endif 167 PetscCall(VecDestroy(&diff)); 168 169 /* MATSHELL does not support MatDiagonalSet after MatScale */ 170 if (strncmp(mattypename, "shell", 5)) { 171 PetscCall(MatDiagonalSet(A, X, INSERT_VALUES)); 172 PetscCall(MatGetDiagonal(A, W1)); 173 PetscCall(VecView(W1, viewer)); 174 } else { 175 PetscCall(PetscViewerASCIIPrintf(viewer, "MatDiagonalSet not tested on MATSHELL\n")); 176 } 177 178 PetscCall(MatDestroy(&E)); 179 PetscCall(VecDestroy(&W1)); 180 PetscCall(VecDestroy(&W2)); 181 PetscCall(VecDestroy(&W3)); 182 PetscFunctionReturn(PETSC_SUCCESS); 183 } 184 185 int main(int argc, char **args) 186 { 187 const PetscInt inds[] = {0, 1}; 188 #if defined(PETSC_USE_COMPLEX) 189 const PetscScalar xvals[] = {PetscCMPLX(11, 4), PetscCMPLX(13, 2)}, yvals[] = {PetscCMPLX(17, 3), PetscCMPLX(19, 1)}, zvals[] = {PetscCMPLX(23, 6), PetscCMPLX(29, 2)}; 190 PetscScalar avals[] = {PetscCMPLX(2, 3), PetscCMPLX(3, 5), PetscCMPLX(5, 4), PetscCMPLX(7, 5)}; 191 #else 192 const PetscScalar xvals[] = {11, 13}, yvals[] = {17, 19}, zvals[] = {23, 29}; 193 PetscScalar avals[] = {2, 3, 5, 7}; 194 #endif 195 Mat A, S, D[4], N; 196 Vec X, Y, Z; 197 User user; 198 PetscInt i; 199 200 PetscFunctionBeginUser; 201 PetscCall(PetscInitialize(&argc, &args, NULL, help)); 202 PetscCall(MatCreateSeqAIJ(PETSC_COMM_WORLD, 2, 2, 2, NULL, &A)); 203 PetscCall(MatSetUp(A)); 204 PetscCall(VecCreateSeq(PETSC_COMM_WORLD, 2, &X)); 205 PetscCall(VecDuplicate(X, &Y)); 206 PetscCall(VecDuplicate(X, &Z)); 207 PetscCall(MatSetValues(A, 2, inds, 2, inds, avals, INSERT_VALUES)); 208 PetscCall(VecSetValues(X, 2, inds, xvals, INSERT_VALUES)); 209 PetscCall(VecSetValues(Y, 2, inds, yvals, INSERT_VALUES)); 210 PetscCall(VecSetValues(Z, 2, inds, zvals, INSERT_VALUES)); 211 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 212 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 213 PetscCall(VecAssemblyBegin(X)); 214 PetscCall(VecAssemblyBegin(Y)); 215 PetscCall(VecAssemblyBegin(Z)); 216 PetscCall(VecAssemblyEnd(X)); 217 PetscCall(VecAssemblyEnd(Y)); 218 PetscCall(VecAssemblyEnd(Z)); 219 220 PetscCall(PetscNew(&user)); 221 user->B = A; 222 223 PetscCall(MatCreateShell(PETSC_COMM_WORLD, 2, 2, 2, 2, user, &S)); 224 PetscCall(MatSetUp(S)); 225 PetscCall(MatShellSetOperation(S, MATOP_VIEW, (PetscErrorCodeFn *)MatView_User)); 226 PetscCall(MatShellSetOperation(S, MATOP_MULT, (PetscErrorCodeFn *)MatMult_User)); 227 PetscCall(MatShellSetOperation(S, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)MatMultTranspose_User)); 228 PetscCall(MatShellSetOperation(S, MATOP_GET_DIAGONAL, (PetscErrorCodeFn *)MatGetDiagonal_User)); 229 230 for (i = 0; i < 4; i++) PetscCall(MatCreateSeqDense(PETSC_COMM_WORLD, 1, 1, &avals[i], &D[i])); 231 PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, D, &N)); 232 PetscCall(MatSetUp(N)); 233 234 PetscCall(TestMatrix(S, X, Y, Z)); 235 PetscCall(TestMatrix(A, X, Y, Z)); 236 PetscCall(TestMatrix(N, X, Y, Z)); 237 238 for (i = 0; i < 4; i++) PetscCall(MatDestroy(&D[i])); 239 PetscCall(MatDestroy(&A)); 240 PetscCall(MatDestroy(&S)); 241 PetscCall(MatDestroy(&N)); 242 PetscCall(VecDestroy(&X)); 243 PetscCall(VecDestroy(&Y)); 244 PetscCall(VecDestroy(&Z)); 245 PetscCall(PetscFree(user)); 246 PetscCall(PetscFinalize()); 247 return 0; 248 } 249 250 /*TEST 251 252 testset: 253 test: 254 suffix: 1 255 requires:!complex 256 test: 257 suffix: 2 258 requires: complex 259 260 TEST*/ 261