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(0); 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(0); 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(0); 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(0); 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 PetscCheckFalse(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 PetscCheckFalse(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 PetscCheckFalse(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 PetscCheckFalse(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 PetscCheckFalse(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(0); 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 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++) { 187 PetscCall(MatCreateSeqDense(PETSC_COMM_WORLD,1,1,&avals[i],&D[i])); 188 } 189 PetscCall(MatCreateNest(PETSC_COMM_WORLD,2,NULL,2,NULL,D,&N)); 190 PetscCall(MatSetUp(N)); 191 192 PetscCall(TestMatrix(S,X,Y,Z)); 193 PetscCall(TestMatrix(A,X,Y,Z)); 194 PetscCall(TestMatrix(N,X,Y,Z)); 195 196 for (i=0; i<4; i++) PetscCall(MatDestroy(&D[i])); 197 PetscCall(MatDestroy(&A)); 198 PetscCall(MatDestroy(&S)); 199 PetscCall(MatDestroy(&N)); 200 PetscCall(VecDestroy(&X)); 201 PetscCall(VecDestroy(&Y)); 202 PetscCall(VecDestroy(&Z)); 203 PetscCall(PetscFree(user)); 204 PetscCall(PetscFinalize()); 205 return 0; 206 } 207 208 /*TEST 209 210 test: 211 212 TEST*/ 213