Actual source code: ex88.c

petsc-3.10.5 2019-03-28
Report Typos and Errors

  2: static char help[] = "Tests MatShift(), MatScale(), and MatDiagonalScale() for SHELL and NEST matrices\n\n";

  4:  #include <petscmat.h>

  6: typedef struct _n_User *User;
  7: struct _n_User {
  8:   Mat B;
  9: };

 11: static PetscErrorCode MatView_User(Mat A,PetscViewer viewer)
 12: {
 13:   User           user;

 17:   MatShellGetContext(A,&user);
 18:   MatView(user->B,viewer);
 19:   return(0);
 20: }

 22: static PetscErrorCode MatMult_User(Mat A,Vec X,Vec Y)
 23: {
 24:   User           user;

 28:   MatShellGetContext(A,&user);
 29:   MatMult(user->B,X,Y);
 30:   return(0);
 31: }

 33: static PetscErrorCode MatMultTranspose_User(Mat A,Vec X,Vec Y)
 34: {
 35:   User           user;

 39:   MatShellGetContext(A,&user);
 40:   MatMultTranspose(user->B,X,Y);
 41:   return(0);
 42: }

 44: static PetscErrorCode MatGetDiagonal_User(Mat A,Vec X)
 45: {
 46:   User           user;

 50:   MatShellGetContext(A,&user);
 51:   MatGetDiagonal(user->B,X);
 52:   return(0);
 53: }

 55: static PetscErrorCode TestMatrix(Mat A,Vec X,Vec Y,Vec Z)
 56: {
 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,norm1,norminf;

 68:   PetscObjectGetType((PetscObject)A,&mattypename);
 69:   PetscViewerASCIIPrintf(viewer,"\nMatrix of type: %s\n",mattypename);
 70:   VecDuplicate(X,&W1);
 71:   VecDuplicate(X,&W2);
 72:   MatScale(A,31);
 73:   MatShift(A,37);
 74:   MatDiagonalScale(A,X,Y);
 75:   MatScale(A,41);
 76:   MatDiagonalScale(A,Y,Z);
 77:   MatComputeExplicitOperator(A,&E);

 79:   MatView(E,viewer);
 80:   PetscViewerASCIIPrintf(viewer,"Testing MatMult + MatMultTranspose\n");
 81:   MatMult(A,Z,W1);
 82:   MatMultTranspose(A,W1,W2);
 83:   VecView(W2,viewer);

 85:   PetscViewerASCIIPrintf(viewer,"Testing MatMultAdd\n");
 86:   VecCreateSeqWithArray(PETSC_COMM_SELF,1,2,multadd,&diff);
 87:   VecSet(W1,-1.0);
 88:   MatMultAdd(A,W1,W1,W2);
 89:   VecView(W2,viewer);
 90:   VecAXPY(W2,-1.0,diff);
 91:   VecNorm(W2,NORM_2,&nrm);
 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

 96:   VecSet(W2,-1.0);
 97:   MatMultAdd(A,W1,W2,W2);
 98:   VecView(W2,viewer);
 99:   VecAXPY(W2,-1.0,diff);
100:   VecNorm(W2,NORM_2,&nrm);
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:   VecDestroy(&diff);

106:   PetscViewerASCIIPrintf(viewer,"Testing MatMultTranposeAdd\n");
107:   VecCreateSeqWithArray(PETSC_COMM_SELF,1,2,multtadd,&diff);

109:   VecSet(W1,-1.0);
110:   MatMultTransposeAdd(A,W1,W1,W2);
111:   VecView(W2,viewer);
112:   VecAXPY(W2,-1.0,diff);
113:   VecNorm(W2,NORM_2,&nrm);
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

118:   VecSet(W2,-1.0);
119:   MatMultTransposeAdd(A,W1,W2,W2);
120:   VecView(W2,viewer);
121:   VecAXPY(W2,-1.0,diff);
122:   VecNorm(W2,NORM_2,&nrm);
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:   VecDestroy(&diff);

128:   PetscViewerASCIIPrintf(viewer,"Testing MatGetDiagonal\n");
129:   MatGetDiagonal(A,W2);
130:   VecView(W2,viewer);
131:   VecCreateSeqWithArray(PETSC_COMM_SELF,1,2,diag,&diff);
132:   VecAXPY(diff,-1.0,W2);
133:   VecNorm(diff,NORM_2,&nrm);
134:   if (nrm > PETSC_SMALL) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"MatGetDiagonal() produces incorrect result");
135:   VecDestroy(&diff);
136:   MatNorm(A,NORM_1,&norm1);
137:   MatNorm(A,NORM_INFINITY,&norminf);
138:   PetscViewerASCIIPrintf(viewer,"Norm_1: %g, Norm_infty %g\n",(double)norm1,(double)norminf);
139:   /* MATSHELL does not support MatDiagonalSet after MatScale */
140:   if (strncmp(mattypename, "shell", 5)) {
141:     MatDiagonalSet(A,X,INSERT_VALUES);
142:     MatGetDiagonal(A,W1);
143:     VecView(W1,viewer);
144:   } else {
145:     PetscViewerASCIIPrintf(viewer,"MatDiagonalSet not tested on MATSHELL\n");
146:   }

148:   MatDestroy(&E);
149:   VecDestroy(&W1);
150:   VecDestroy(&W2);
151:   return(0);
152: }

154: int main(int argc,char **args)
155: {
156:   const PetscScalar xvals[] = {11,13},yvals[] = {17,19},zvals[] = {23,29};
157:   const PetscInt    inds[]  = {0,1};
158:   PetscScalar       avals[] = {2,3,5,7};
159:   Mat               A,S,D[4],N;
160:   Vec               X,Y,Z;
161:   User              user;
162:   PetscInt          i;
163:   PetscErrorCode    ierr;

165:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
166:   MatCreateSeqAIJ(PETSC_COMM_WORLD,2,2,2,NULL,&A);
167:   MatSetUp(A);
168:   MatSetValues(A,2,inds,2,inds,avals,INSERT_VALUES);
169:   VecCreateSeq(PETSC_COMM_WORLD,2,&X);
170:   VecDuplicate(X,&Y);
171:   VecDuplicate(X,&Z);
172:   VecSetValues(X,2,inds,xvals,INSERT_VALUES);
173:   VecSetValues(Y,2,inds,yvals,INSERT_VALUES);
174:   VecSetValues(Z,2,inds,zvals,INSERT_VALUES);
175:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
176:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
177:   VecAssemblyBegin(X);
178:   VecAssemblyBegin(Y);
179:   VecAssemblyBegin(Z);
180:   VecAssemblyEnd(X);
181:   VecAssemblyEnd(Y);
182:   VecAssemblyEnd(Z);

184:   PetscNew(&user);
185:   user->B = A;

187:   MatCreateShell(PETSC_COMM_WORLD,2,2,2,2,user,&S);
188:   MatSetUp(S);
189:   MatShellSetOperation(S,MATOP_VIEW,(void (*)(void))MatView_User);
190:   MatShellSetOperation(S,MATOP_MULT,(void (*)(void))MatMult_User);
191:   MatShellSetOperation(S,MATOP_MULT_TRANSPOSE,(void (*)(void))MatMultTranspose_User);
192:   MatShellSetOperation(S,MATOP_GET_DIAGONAL,(void (*)(void))MatGetDiagonal_User);

194:   for (i=0; i<4; i++) {
195:     MatCreateSeqDense(PETSC_COMM_WORLD,1,1,&avals[i],&D[i]);
196:   }
197:   MatCreateNest(PETSC_COMM_WORLD,2,NULL,2,NULL,D,&N);
198:   MatSetUp(N);

200:   TestMatrix(S,X,Y,Z);
201:   TestMatrix(A,X,Y,Z);
202:   TestMatrix(N,X,Y,Z);

204:   for (i=0; i<4; i++) {MatDestroy(&D[i]);}
205:   MatDestroy(&A);
206:   MatDestroy(&S);
207:   MatDestroy(&N);
208:   VecDestroy(&X);
209:   VecDestroy(&Y);
210:   VecDestroy(&Z);
211:   PetscFree(user);
212:   PetscFinalize();
213:   return ierr;
214: }


217: /*TEST

219:    test:

221: TEST*/