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