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 CHKERRQ(MatShellGetContext(A,&user)); 17 CHKERRQ(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 CHKERRQ(MatShellGetContext(A,&user)); 27 CHKERRQ(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 CHKERRQ(MatShellGetContext(A,&user)); 37 CHKERRQ(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 CHKERRQ(MatShellGetContext(A,&user)); 47 CHKERRQ(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 CHKERRQ(PetscObjectGetType((PetscObject)A,&mattypename)); 64 CHKERRQ(PetscViewerASCIIPrintf(viewer,"\nMatrix of type: %s\n",mattypename)); 65 CHKERRQ(VecDuplicate(X,&W1)); 66 CHKERRQ(VecDuplicate(X,&W2)); 67 CHKERRQ(MatScale(A,31)); 68 CHKERRQ(MatShift(A,37)); 69 CHKERRQ(MatDiagonalScale(A,X,Y)); 70 CHKERRQ(MatScale(A,41)); 71 CHKERRQ(MatDiagonalScale(A,Y,Z)); 72 CHKERRQ(MatComputeOperator(A,MATDENSE,&E)); 73 74 CHKERRQ(MatView(E,viewer)); 75 CHKERRQ(PetscViewerASCIIPrintf(viewer,"Testing MatMult + MatMultTranspose\n")); 76 CHKERRQ(MatMult(A,Z,W1)); 77 CHKERRQ(MatMultTranspose(A,W1,W2)); 78 CHKERRQ(VecView(W2,viewer)); 79 80 CHKERRQ(PetscViewerASCIIPrintf(viewer,"Testing MatMultAdd\n")); 81 CHKERRQ(VecCreateSeqWithArray(PETSC_COMM_SELF,1,2,multadd,&diff)); 82 CHKERRQ(VecSet(W1,-1.0)); 83 CHKERRQ(MatMultAdd(A,W1,W1,W2)); 84 CHKERRQ(VecView(W2,viewer)); 85 CHKERRQ(VecAXPY(W2,-1.0,diff)); 86 CHKERRQ(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 CHKERRQ(VecSet(W2,-1.0)); 92 CHKERRQ(MatMultAdd(A,W1,W2,W2)); 93 CHKERRQ(VecView(W2,viewer)); 94 CHKERRQ(VecAXPY(W2,-1.0,diff)); 95 CHKERRQ(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 CHKERRQ(VecDestroy(&diff)); 100 101 CHKERRQ(PetscViewerASCIIPrintf(viewer,"Testing MatMultTransposeAdd\n")); 102 CHKERRQ(VecCreateSeqWithArray(PETSC_COMM_SELF,1,2,multtadd,&diff)); 103 104 CHKERRQ(VecSet(W1,-1.0)); 105 CHKERRQ(MatMultTransposeAdd(A,W1,W1,W2)); 106 CHKERRQ(VecView(W2,viewer)); 107 CHKERRQ(VecAXPY(W2,-1.0,diff)); 108 CHKERRQ(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 CHKERRQ(VecSet(W2,-1.0)); 114 CHKERRQ(MatMultTransposeAdd(A,W1,W2,W2)); 115 CHKERRQ(VecView(W2,viewer)); 116 CHKERRQ(VecAXPY(W2,-1.0,diff)); 117 CHKERRQ(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 CHKERRQ(VecDestroy(&diff)); 122 123 CHKERRQ(PetscViewerASCIIPrintf(viewer,"Testing MatGetDiagonal\n")); 124 CHKERRQ(MatGetDiagonal(A,W2)); 125 CHKERRQ(VecView(W2,viewer)); 126 CHKERRQ(VecCreateSeqWithArray(PETSC_COMM_SELF,1,2,diag,&diff)); 127 CHKERRQ(VecAXPY(diff,-1.0,W2)); 128 CHKERRQ(VecNorm(diff,NORM_2,&nrm)); 129 PetscCheckFalse(nrm > PETSC_SMALL,PETSC_COMM_SELF,PETSC_ERR_PLIB,"MatGetDiagonal() produces incorrect result"); 130 CHKERRQ(VecDestroy(&diff)); 131 132 /* MATSHELL does not support MatDiagonalSet after MatScale */ 133 if (strncmp(mattypename, "shell", 5)) { 134 CHKERRQ(MatDiagonalSet(A,X,INSERT_VALUES)); 135 CHKERRQ(MatGetDiagonal(A,W1)); 136 CHKERRQ(VecView(W1,viewer)); 137 } else { 138 CHKERRQ(PetscViewerASCIIPrintf(viewer,"MatDiagonalSet not tested on MATSHELL\n")); 139 } 140 141 CHKERRQ(MatDestroy(&E)); 142 CHKERRQ(VecDestroy(&W1)); 143 CHKERRQ(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 PetscErrorCode ierr; 157 158 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 159 CHKERRQ(MatCreateSeqAIJ(PETSC_COMM_WORLD,2,2,2,NULL,&A)); 160 CHKERRQ(MatSetUp(A)); 161 CHKERRQ(MatSetValues(A,2,inds,2,inds,avals,INSERT_VALUES)); 162 CHKERRQ(VecCreateSeq(PETSC_COMM_WORLD,2,&X)); 163 CHKERRQ(VecDuplicate(X,&Y)); 164 CHKERRQ(VecDuplicate(X,&Z)); 165 CHKERRQ(VecSetValues(X,2,inds,xvals,INSERT_VALUES)); 166 CHKERRQ(VecSetValues(Y,2,inds,yvals,INSERT_VALUES)); 167 CHKERRQ(VecSetValues(Z,2,inds,zvals,INSERT_VALUES)); 168 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 169 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 170 CHKERRQ(VecAssemblyBegin(X)); 171 CHKERRQ(VecAssemblyBegin(Y)); 172 CHKERRQ(VecAssemblyBegin(Z)); 173 CHKERRQ(VecAssemblyEnd(X)); 174 CHKERRQ(VecAssemblyEnd(Y)); 175 CHKERRQ(VecAssemblyEnd(Z)); 176 177 CHKERRQ(PetscNew(&user)); 178 user->B = A; 179 180 CHKERRQ(MatCreateShell(PETSC_COMM_WORLD,2,2,2,2,user,&S)); 181 CHKERRQ(MatSetUp(S)); 182 CHKERRQ(MatShellSetOperation(S,MATOP_VIEW,(void (*)(void))MatView_User)); 183 CHKERRQ(MatShellSetOperation(S,MATOP_MULT,(void (*)(void))MatMult_User)); 184 CHKERRQ(MatShellSetOperation(S,MATOP_MULT_TRANSPOSE,(void (*)(void))MatMultTranspose_User)); 185 CHKERRQ(MatShellSetOperation(S,MATOP_GET_DIAGONAL,(void (*)(void))MatGetDiagonal_User)); 186 187 for (i=0; i<4; i++) { 188 CHKERRQ(MatCreateSeqDense(PETSC_COMM_WORLD,1,1,&avals[i],&D[i])); 189 } 190 CHKERRQ(MatCreateNest(PETSC_COMM_WORLD,2,NULL,2,NULL,D,&N)); 191 CHKERRQ(MatSetUp(N)); 192 193 CHKERRQ(TestMatrix(S,X,Y,Z)); 194 CHKERRQ(TestMatrix(A,X,Y,Z)); 195 CHKERRQ(TestMatrix(N,X,Y,Z)); 196 197 for (i=0; i<4; i++) CHKERRQ(MatDestroy(&D[i])); 198 CHKERRQ(MatDestroy(&A)); 199 CHKERRQ(MatDestroy(&S)); 200 CHKERRQ(MatDestroy(&N)); 201 CHKERRQ(VecDestroy(&X)); 202 CHKERRQ(VecDestroy(&Y)); 203 CHKERRQ(VecDestroy(&Z)); 204 CHKERRQ(PetscFree(user)); 205 ierr = PetscFinalize(); 206 return ierr; 207 } 208 209 /*TEST 210 211 test: 212 213 TEST*/ 214