Actual source code: ex88.c


  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;

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

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

 24:   MatShellGetContext(A,&user);
 25:   MatMult(user->B,X,Y);
 26:   return 0;
 27: }

 29: static PetscErrorCode MatMultTranspose_User(Mat A,Vec X,Vec Y)
 30: {
 31:   User           user;

 33:   MatShellGetContext(A,&user);
 34:   MatMultTranspose(user->B,X,Y);
 35:   return 0;
 36: }

 38: static PetscErrorCode MatGetDiagonal_User(Mat A,Vec X)
 39: {
 40:   User           user;

 42:   MatShellGetContext(A,&user);
 43:   MatGetDiagonal(user->B,X);
 44:   return 0;
 45: }

 47: static PetscErrorCode TestMatrix(Mat A,Vec X,Vec Y,Vec Z)
 48: {
 49:   Vec            W1,W2,diff;
 50:   Mat            E;
 51:   const char     *mattypename;
 52:   PetscViewer    viewer = PETSC_VIEWER_STDOUT_WORLD;
 53:   PetscScalar    diag[2]     = { 2.9678190300000000e+08, 1.4173141580000000e+09};
 54:   PetscScalar    multadd[2]  = {-6.8966198500000000e+08,-2.0310609940000000e+09};
 55:   PetscScalar    multtadd[2] = {-9.1052873900000000e+08,-1.8101942400000000e+09};
 56:   PetscReal      nrm;

 58:   PetscObjectGetType((PetscObject)A,&mattypename);
 59:   PetscViewerASCIIPrintf(viewer,"\nMatrix of type: %s\n",mattypename);
 60:   VecDuplicate(X,&W1);
 61:   VecDuplicate(X,&W2);
 62:   MatScale(A,31);
 63:   MatShift(A,37);
 64:   MatDiagonalScale(A,X,Y);
 65:   MatScale(A,41);
 66:   MatDiagonalScale(A,Y,Z);
 67:   MatComputeOperator(A,MATDENSE,&E);

 69:   MatView(E,viewer);
 70:   PetscViewerASCIIPrintf(viewer,"Testing MatMult + MatMultTranspose\n");
 71:   MatMult(A,Z,W1);
 72:   MatMultTranspose(A,W1,W2);
 73:   VecView(W2,viewer);

 75:   PetscViewerASCIIPrintf(viewer,"Testing MatMultAdd\n");
 76:   VecCreateSeqWithArray(PETSC_COMM_SELF,1,2,multadd,&diff);
 77:   VecSet(W1,-1.0);
 78:   MatMultAdd(A,W1,W1,W2);
 79:   VecView(W2,viewer);
 80:   VecAXPY(W2,-1.0,diff);
 81:   VecNorm(W2,NORM_2,&nrm);
 82: #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
 84: #endif

 86:   VecSet(W2,-1.0);
 87:   MatMultAdd(A,W1,W2,W2);
 88:   VecView(W2,viewer);
 89:   VecAXPY(W2,-1.0,diff);
 90:   VecNorm(W2,NORM_2,&nrm);
 91: #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
 93: #endif
 94:   VecDestroy(&diff);

 96:   PetscViewerASCIIPrintf(viewer,"Testing MatMultTransposeAdd\n");
 97:   VecCreateSeqWithArray(PETSC_COMM_SELF,1,2,multtadd,&diff);

 99:   VecSet(W1,-1.0);
100:   MatMultTransposeAdd(A,W1,W1,W2);
101:   VecView(W2,viewer);
102:   VecAXPY(W2,-1.0,diff);
103:   VecNorm(W2,NORM_2,&nrm);
104: #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
106: #endif

108:   VecSet(W2,-1.0);
109:   MatMultTransposeAdd(A,W1,W2,W2);
110:   VecView(W2,viewer);
111:   VecAXPY(W2,-1.0,diff);
112:   VecNorm(W2,NORM_2,&nrm);
113: #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
115: #endif
116:   VecDestroy(&diff);

118:   PetscViewerASCIIPrintf(viewer,"Testing MatGetDiagonal\n");
119:   MatGetDiagonal(A,W2);
120:   VecView(W2,viewer);
121:   VecCreateSeqWithArray(PETSC_COMM_SELF,1,2,diag,&diff);
122:   VecAXPY(diff,-1.0,W2);
123:   VecNorm(diff,NORM_2,&nrm);
125:   VecDestroy(&diff);

127:   /* MATSHELL does not support MatDiagonalSet after MatScale */
128:   if (strncmp(mattypename, "shell", 5)) {
129:     MatDiagonalSet(A,X,INSERT_VALUES);
130:     MatGetDiagonal(A,W1);
131:     VecView(W1,viewer);
132:   } else {
133:     PetscViewerASCIIPrintf(viewer,"MatDiagonalSet not tested on MATSHELL\n");
134:   }

136:   MatDestroy(&E);
137:   VecDestroy(&W1);
138:   VecDestroy(&W2);
139:   return 0;
140: }

142: int main(int argc,char **args)
143: {
144:   const PetscScalar xvals[] = {11,13},yvals[] = {17,19},zvals[] = {23,29};
145:   const PetscInt    inds[]  = {0,1};
146:   PetscScalar       avals[] = {2,3,5,7};
147:   Mat               A,S,D[4],N;
148:   Vec               X,Y,Z;
149:   User              user;
150:   PetscInt          i;

152:   PetscInitialize(&argc,&args,(char*)0,help);
153:   MatCreateSeqAIJ(PETSC_COMM_WORLD,2,2,2,NULL,&A);
154:   MatSetUp(A);
155:   MatSetValues(A,2,inds,2,inds,avals,INSERT_VALUES);
156:   VecCreateSeq(PETSC_COMM_WORLD,2,&X);
157:   VecDuplicate(X,&Y);
158:   VecDuplicate(X,&Z);
159:   VecSetValues(X,2,inds,xvals,INSERT_VALUES);
160:   VecSetValues(Y,2,inds,yvals,INSERT_VALUES);
161:   VecSetValues(Z,2,inds,zvals,INSERT_VALUES);
162:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
163:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
164:   VecAssemblyBegin(X);
165:   VecAssemblyBegin(Y);
166:   VecAssemblyBegin(Z);
167:   VecAssemblyEnd(X);
168:   VecAssemblyEnd(Y);
169:   VecAssemblyEnd(Z);

171:   PetscNew(&user);
172:   user->B = A;

174:   MatCreateShell(PETSC_COMM_WORLD,2,2,2,2,user,&S);
175:   MatSetUp(S);
176:   MatShellSetOperation(S,MATOP_VIEW,(void (*)(void))MatView_User);
177:   MatShellSetOperation(S,MATOP_MULT,(void (*)(void))MatMult_User);
178:   MatShellSetOperation(S,MATOP_MULT_TRANSPOSE,(void (*)(void))MatMultTranspose_User);
179:   MatShellSetOperation(S,MATOP_GET_DIAGONAL,(void (*)(void))MatGetDiagonal_User);

181:   for (i=0; i<4; i++) {
182:     MatCreateSeqDense(PETSC_COMM_WORLD,1,1,&avals[i],&D[i]);
183:   }
184:   MatCreateNest(PETSC_COMM_WORLD,2,NULL,2,NULL,D,&N);
185:   MatSetUp(N);

187:   TestMatrix(S,X,Y,Z);
188:   TestMatrix(A,X,Y,Z);
189:   TestMatrix(N,X,Y,Z);

191:   for (i=0; i<4; i++) MatDestroy(&D[i]);
192:   MatDestroy(&A);
193:   MatDestroy(&S);
194:   MatDestroy(&N);
195:   VecDestroy(&X);
196:   VecDestroy(&Y);
197:   VecDestroy(&Z);
198:   PetscFree(user);
199:   PetscFinalize();
200:   return 0;
201: }

203: /*TEST

205:    test:

207: TEST*/