Actual source code: ex94.c

  2: static char help[] = "Tests sequential and parallel MatMatMult() and MatPtAP(), sequential MatMatMultTranspose()\n\
  3: Input arguments are:\n\
  4:   -f0 <input_file> -f1 <input_file> -f2 <input_file> -f3 <input_file> : file to load\n\n";
  5: /* e.g., mpiexec -n 3 ./ex94 -f0 $D/medium -f1 $D/medium -f2 $D/arco1 -f3 $D/arco1 */

 7:  #include petscmat.h

 11: int main(int argc,char **args)
 12: {
 13:   Mat            A,A_save,B,P,C,C1;
 14:   Vec            x,v1,v2;
 15:   PetscViewer    viewer;
 17:   PetscMPIInt    size,rank;
 18:   PetscInt       i,m,n,j,*idxn,M,N,nzp;
 19:   PetscReal      norm,norm_tmp,tol=1.e-8,fill=4.0;
 20:   PetscRandom    rdm;
 21:   char           file[4][128];
 22:   PetscTruth     flg,preload = PETSC_TRUE;
 23:   PetscScalar    *a,rval,alpha,none = -1.0;
 24:   PetscTruth     Test_MatMatMult=PETSC_TRUE,Test_MatMatMultTr=PETSC_TRUE;
 25:   Vec            v3,v4,v5;
 26:   PetscInt       pm,pn,pM,pN;
 27:   PetscTruth     Test_MatPtAP=PETSC_TRUE;
 28: 
 29:   PetscInitialize(&argc,&args,(char *)0,help);
 30:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 31:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

 33:   /*  Load the matrices A and B */
 34:   PetscOptionsGetString(PETSC_NULL,"-f0",file[0],127,&flg);
 35:   if (!flg) SETERRQ(1,"Must indicate a file name for small matrix A with the -f0 option.");
 36:   PetscOptionsGetString(PETSC_NULL,"-f1",file[1],127,&flg);
 37:   if (!flg) SETERRQ(1,"Must indicate a file name for small matrix B with the -f1 option.");
 38:   PetscOptionsGetString(PETSC_NULL,"-f2",file[2],127,&flg);
 39:   if (!flg) {
 40:     preload = PETSC_FALSE;
 41:   } else {
 42:     PetscOptionsGetString(PETSC_NULL,"-f3",file[3],127,&flg);
 43:     if (!flg) SETERRQ(1,"Must indicate a file name for test matrix B with the -f3 option.");
 44:   }

 46:   PreLoadBegin(preload,"Load system");
 47:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[2*PreLoadIt],FILE_MODE_READ,&viewer);
 48:   MatLoad(viewer,MATAIJ,&A_save);
 49:   PetscViewerDestroy(viewer);

 51:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[2*PreLoadIt+1],FILE_MODE_READ,&viewer);
 52:   MatLoad(viewer,MATAIJ,&B);
 53:   PetscViewerDestroy(viewer);

 55:   MatGetSize(B,&M,&N);
 56:   nzp  = (PetscInt)(0.1*M);
 57:   PetscMalloc((nzp+1)*(sizeof(PetscInt)+sizeof(PetscScalar)),&idxn);
 58:   a    = (PetscScalar*)(idxn + nzp);
 59: 
 60:   /* Create vectors v1 and v2 that are compatible with A_save */
 61:   VecCreate(PETSC_COMM_WORLD,&v1);
 62:   MatGetLocalSize(A_save,&m,PETSC_NULL);
 63:   VecSetSizes(v1,m,PETSC_DECIDE);
 64:   VecSetFromOptions(v1);
 65:   VecDuplicate(v1,&v2);

 67:   PetscRandomCreate(PETSC_COMM_WORLD,&rdm);
 68:   PetscRandomSetFromOptions(rdm);
 69:   PetscOptionsGetReal(PETSC_NULL,"-fill",&fill,PETSC_NULL);

 71:   /* Test MatMatMult() */
 72:   /*-------------------*/
 73:   if (Test_MatMatMult){
 74:     MatDuplicate(A_save,MAT_COPY_VALUES,&A);
 75:     MatMatMult(A,B,MAT_INITIAL_MATRIX,fill,&C);

 77:     /* Test MAT_REUSE_MATRIX - reuse symbolic C */
 78:     alpha=1.0;
 79:     for (i=0; i<2; i++){
 80:       alpha -=0.1;
 81:       MatScale(A,alpha);
 82:       MatMatMult(A,B,MAT_REUSE_MATRIX,fill,&C);
 83:     }

 85:     /* Create vector x that is compatible with B */
 86:     VecCreate(PETSC_COMM_WORLD,&x);
 87:     MatGetLocalSize(B,PETSC_NULL,&n);
 88:     VecSetSizes(x,n,PETSC_DECIDE);
 89:     VecSetFromOptions(x);

 91:     norm = 0.0;
 92:     for (i=0; i<10; i++) {
 93:       VecSetRandom(x,rdm);
 94:       MatMult(B,x,v1);
 95:       MatMult(A,v1,v2);  /* v2 = A*B*x */
 96:       MatMult(C,x,v1);   /* v1 = C*x   */
 97:       VecAXPY(v1,none,v2);
 98:       VecNorm(v1,NORM_2,&norm_tmp);
 99:       if (norm_tmp > norm) norm = norm_tmp;
100:     }
101:     if (norm >= tol) {
102:       PetscPrintf(PETSC_COMM_SELF,"Error: MatMatMult(), |v1 - v2|: %G\n",norm);
103:     }

105:     VecDestroy(x);
106:     MatDestroy(A);

108:     /* Test MatDuplicate() of C */
109:     MatDuplicate(C,MAT_COPY_VALUES,&C1);
110:     MatDestroy(C1);
111:     MatDestroy(C);
112:   } /* if (Test_MatMatMult) */

114:   /* Test MatMatMultTranspose() */
115:   /*----------------------------*/
116:   if (size>1) Test_MatMatMultTr = PETSC_FALSE;
117:   if (Test_MatMatMultTr){
118:     PetscInt PN;
119:     /* MatGetSize(B,&M,&N); */
120:     PN   = M/2;
121:     nzp  = 5; /* num of nonzeros in each row of P */
122:     MatCreate(PETSC_COMM_WORLD,&P);
123:     MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,M,PN);
124:     MatSetType(P,MATAIJ);
125:     MatSeqAIJSetPreallocation(P,nzp,PETSC_NULL);
126:     MatMPIAIJSetPreallocation(P,nzp,PETSC_NULL,nzp,PETSC_NULL);
127:     for (i=0; i<nzp; i++){
128:       PetscRandomGetValue(rdm,&a[i]);
129:     }
130:     for (i=0; i<M; i++){
131:       for (j=0; j<nzp; j++){
132:         PetscRandomGetValue(rdm,&rval);
133:         idxn[j] = (PetscInt)(PetscRealPart(rval)*PN);
134:       }
135:       MatSetValues(P,1,&i,nzp,idxn,a,ADD_VALUES);
136:     }
137:     MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
138:     MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);
139: 
140:     MatMatMultTranspose(P,B,MAT_INITIAL_MATRIX,fill,&C);

142:     /* Test MAT_REUSE_MATRIX - reuse symbolic C */
143:     alpha=1.0;
144:     for (i=0; i<2; i++){
145:       alpha -=0.1;
146:       MatScale(P,alpha);
147:       MatMatMultTranspose(P,B,MAT_REUSE_MATRIX,fill,&C);
148:     }

150:     /* Create vector x, v5 that are compatible with B */
151:     VecCreate(PETSC_COMM_WORLD,&x);
152:     MatGetLocalSize(B,&m,&n);
153:     VecSetSizes(x,n,PETSC_DECIDE);
154:     VecSetFromOptions(x);

156:     VecCreate(PETSC_COMM_WORLD,&v5);
157:     VecSetSizes(v5,m,PETSC_DECIDE);
158:     VecSetFromOptions(v5);
159: 
160:     MatGetLocalSize(P,PETSC_NULL,&n);
161:     VecCreate(PETSC_COMM_WORLD,&v3);
162:     VecSetSizes(v3,n,PETSC_DECIDE);
163:     VecSetFromOptions(v3);
164:     VecDuplicate(v3,&v4);

166:     norm = 0.0;
167:     for (i=0; i<10; i++) {
168:       VecSetRandom(x,rdm);
169:       MatMult(B,x,v5);            /* v5 = B*x   */
170:       MatMultTranspose(P,v5,v3);  /* v3 = Pt*B*x */
171:       MatMult(C,x,v4);            /* v4 = C*x   */
172:       VecAXPY(v4,none,v3);
173:       VecNorm(v4,NORM_2,&norm_tmp);
174:       if (norm_tmp > norm) norm = norm_tmp;
175:     }
176:     if (norm >= tol) {
177:       PetscPrintf(PETSC_COMM_SELF,"Error: MatMatMultTr(), |v3 - v4|: %G\n",norm);
178:     }
179:     MatDestroy(P);
180:     MatDestroy(C);
181:     VecDestroy(v3);
182:     VecDestroy(v4);
183:     VecDestroy(v5);
184:     VecDestroy(x);
185:   }

187:   /* Test MatPtAP() */
188:   /*----------------------*/
189:   if (Test_MatPtAP){
190:     PetscInt PN;
191:     MatDuplicate(A_save,MAT_COPY_VALUES,&A);
192:     MatGetSize(A,&M,&N);
193:     MatGetLocalSize(A,&m,&n);
194:     /* PetscPrintf(PETSC_COMM_SELF,"[%d] A: %d,%d, %d,%d\n",rank,m,n,M,N); */

196:     PN   = M/2;
197:     nzp  = (PetscInt)(0.1*PN); /* num of nozeros in each row of P */
198:     MatCreate(PETSC_COMM_WORLD,&P);
199:     MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,N,PN);
200:     MatSetType(P,MATAIJ);
201:     MatSeqAIJSetPreallocation(P,nzp,PETSC_NULL);
202:     MatMPIAIJSetPreallocation(P,nzp,PETSC_NULL,nzp,PETSC_NULL);
203:     for (i=0; i<nzp; i++){
204:       PetscRandomGetValue(rdm,&a[i]);
205:     }
206:     for (i=0; i<M; i++){
207:       for (j=0; j<nzp; j++){
208:         PetscRandomGetValue(rdm,&rval);
209:         idxn[j] = (PetscInt)(PetscRealPart(rval)*PN);
210:       }
211:       MatSetValues(P,1,&i,nzp,idxn,a,ADD_VALUES);
212:     }
213:     MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
214:     MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);

216:     /* MatView(P,PETSC_VIEWER_STDOUT_WORLD); */
217:     MatGetSize(P,&pM,&pN);
218:     MatGetLocalSize(P,&pm,&pn);
219:     /* PetscPrintf(PETSC_COMM_SELF," [%d] P, %d, %d, %d,%d\n",rank,pm,pn,pM,pN); */
220:     MatPtAP(A,P,MAT_INITIAL_MATRIX,fill,&C);

222:     /* Test MAT_REUSE_MATRIX - reuse symbolic C */
223:     alpha=1.0;
224:     for (i=0; i<2; i++){
225:       alpha -=0.1;
226:       MatScale(A,alpha);
227:       MatPtAP(A,P,MAT_REUSE_MATRIX,fill,&C);
228:     }

230:     /* Create vector x that is compatible with P */
231:     VecCreate(PETSC_COMM_WORLD,&x);
232:     MatGetLocalSize(P,&m,&n);
233:     VecSetSizes(x,n,PETSC_DECIDE);
234:     VecSetFromOptions(x);
235: 
236:     VecCreate(PETSC_COMM_WORLD,&v3);
237:     VecSetSizes(v3,n,PETSC_DECIDE);
238:     VecSetFromOptions(v3);
239:     VecDuplicate(v3,&v4);

241:     norm = 0.0;
242:     for (i=0; i<10; i++) {
243:       VecSetRandom(x,rdm);
244:       MatMult(P,x,v1);
245:       MatMult(A,v1,v2);  /* v2 = A*P*x */

247:       MatMultTranspose(P,v2,v3); /* v3 = Pt*A*P*x */
248:       MatMult(C,x,v4);           /* v3 = C*x   */
249:       VecAXPY(v4,none,v3);
250:       VecNorm(v4,NORM_2,&norm_tmp);
251:       if (norm_tmp > norm) norm = norm_tmp;
252:     }
253:     if (norm >= tol) {
254:       PetscPrintf(PETSC_COMM_SELF,"Error: MatPtAP(), |v1 - v2|: %G\n",norm);
255:     }
256: 
257:     MatDestroy(A);
258:     MatDestroy(P);
259:     MatDestroy(C);
260:     VecDestroy(v3);
261:     VecDestroy(v4);
262:     VecDestroy(x);
263:   } /* if (Test_MatPtAP) */

265:   /* Destroy objects */
266:   VecDestroy(v1);
267:   VecDestroy(v2);
268:   PetscRandomDestroy(rdm);
269:   PetscFree(idxn);
270: 
271:   MatDestroy(A_save);
272:   MatDestroy(B);

274:   PreLoadEnd();
275:   PetscFinalize();

277:   return 0;
278: }