Actual source code: ex94.c

petsc-3.4.2 2013-07-02
  2: static char help[] = "Tests sequential and parallel MatMatMult() and MatPtAP(), MatTransposeMatMult(), sequential MatMatTransposeMult(), MatRARt()\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: /* Example of usage:
  6:    ./ex94 -f0 <A_binary> -f1 <B_binary> -matmatmult_mat_view ::ascii_info -matmatmulttr_mat_view
  7:    mpiexec -n 3 ./ex94 -f0 medium -f1 medium -f2 arco1 -f3 arco1 -matmatmult_mat_view
  8: */

 10: #include <petscmat.h>

 14: int main(int argc,char **args)
 15: {
 16:   Mat            A,A_save,B,P,R,C,C1;
 17:   Vec            x,v1,v2,v3,v4;
 18:   PetscViewer    viewer;
 20:   PetscMPIInt    size,rank;
 21:   PetscInt       i,m,n,j,*idxn,M,N,nzp,rstart,rend;
 22:   PetscReal      norm,norm_abs,norm_tmp,tol=1.e-8,fill=4.0;
 23:   PetscRandom    rdm;
 24:   char           file[4][128];
 25:   PetscBool      flg,preload = PETSC_TRUE;
 26:   PetscScalar    *a,rval,alpha,none = -1.0;
 27:   PetscBool      Test_MatMatMult=PETSC_TRUE,Test_MatMatTr=PETSC_TRUE,Test_MatPtAP=PETSC_TRUE,Test_MatRARt=PETSC_TRUE,Test_MatMatMatMult=PETSC_TRUE;
 28:   PetscInt       pm,pn,pM,pN;
 29:   MatInfo        info;

 31:   PetscInitialize(&argc,&args,(char*)0,help);
 32:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 33:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

 35:   PetscOptionsGetReal(NULL,"-fill",&fill,NULL);

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

 50:   PetscPreLoadBegin(preload,"Load system");
 51:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[2*PetscPreLoadIt],FILE_MODE_READ,&viewer);
 52:   MatCreate(PETSC_COMM_WORLD,&A_save);
 53:   MatLoad(A_save,viewer);
 54:   PetscViewerDestroy(&viewer);

 56:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[2*PetscPreLoadIt+1],FILE_MODE_READ,&viewer);
 57:   MatCreate(PETSC_COMM_WORLD,&B);
 58:   MatLoad(B,viewer);
 59:   PetscViewerDestroy(&viewer);

 61:   MatGetSize(B,&M,&N);
 62:   nzp  = PetscMax((PetscInt)(0.1*M),5);
 63:   PetscMalloc((nzp+1)*(sizeof(PetscInt)+sizeof(PetscScalar)),&idxn);
 64:   a    = (PetscScalar*)(idxn + nzp);

 66:   /* Create vectors v1 and v2 that are compatible with A_save */
 67:   VecCreate(PETSC_COMM_WORLD,&v1);
 68:   MatGetLocalSize(A_save,&m,NULL);
 69:   VecSetSizes(v1,m,PETSC_DECIDE);
 70:   VecSetFromOptions(v1);
 71:   VecDuplicate(v1,&v2);

 73:   PetscRandomCreate(PETSC_COMM_WORLD,&rdm);
 74:   PetscRandomSetFromOptions(rdm);
 75:   PetscOptionsGetReal(NULL,"-fill",&fill,NULL);

 77:   /* Test MatMatMult() */
 78:   /*-------------------*/
 79:   if (Test_MatMatMult) {
 80:     MatDuplicate(A_save,MAT_COPY_VALUES,&A);
 81:     MatMatMult(A,B,MAT_INITIAL_MATRIX,fill,&C);
 82:     MatSetOptionsPrefix(C,"matmatmult_"); /* enable option '-matmatmult_' for matrix C */
 83:     MatGetInfo(C,MAT_GLOBAL_SUM,&info);
 84:     /* PetscPrintf(PETSC_COMM_WORLD,"MatMatMult: nz_allocated = %g; nz_used = %g; nz_unneeded = %g\n",info.nz_allocated,info.nz_used, info.nz_unneeded); */

 86:     /* Test MAT_REUSE_MATRIX - reuse symbolic C */
 87:     alpha=1.0;
 88:     for (i=0; i<2; i++) {
 89:       alpha -=0.1;
 90:       MatScale(A,alpha);
 91:       MatMatMult(A,B,MAT_REUSE_MATRIX,fill,&C);
 92:     }

 94:     /* Create vector x that is compatible with B */
 95:     VecCreate(PETSC_COMM_WORLD,&x);
 96:     MatGetLocalSize(B,NULL,&n);
 97:     VecSetSizes(x,n,PETSC_DECIDE);
 98:     VecSetFromOptions(x);

100:     norm = 0.0;
101:     for (i=0; i<10; i++) {
102:       VecSetRandom(x,rdm);
103:       MatMult(B,x,v1);
104:       MatMult(A,v1,v2);  /* v2 = A*B*x */
105:       MatMult(C,x,v1);   /* v1 = C*x   */
106:       VecNorm(v1,NORM_2,&norm_abs);
107:       VecAXPY(v1,none,v2);
108:       VecNorm(v1,NORM_2,&norm_tmp);

110:       norm_tmp /= norm_abs;
111:       if (norm_tmp > norm) norm = norm_tmp;
112:     }
113:     if (norm >= tol) {
114:       PetscPrintf(PETSC_COMM_SELF,"Error: MatMatMult(), |v1 - v2|: %G\n",norm);
115:     }

117:     VecDestroy(&x);
118:     MatDestroy(&A);

120:     /* Test MatDuplicate() of C */
121:     MatDuplicate(C,MAT_COPY_VALUES,&C1);
122:     MatDestroy(&C1);
123:     MatDestroy(&C);
124:   } /* if (Test_MatMatMult) */

126:   /* Test MatTransposeMatMult() and MatMatTransposeMult() */
127:   /*------------------------------------------------------*/
128:   if (Test_MatMatTr) {
129:     /* Create P */
130:     PetscInt PN,rstart,rend;
131:     PN   = M/2;
132:     nzp  = 5; /* num of nonzeros in each row of P */
133:     MatCreate(PETSC_COMM_WORLD,&P);
134:     MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,M,PN);
135:     MatSetType(P,MATAIJ);
136:     MatSeqAIJSetPreallocation(P,nzp,NULL);
137:     MatMPIAIJSetPreallocation(P,nzp,NULL,nzp,NULL);
138:     MatGetOwnershipRange(P,&rstart,&rend);
139:     for (i=0; i<nzp; i++) {
140:       PetscRandomGetValue(rdm,&a[i]);
141:     }
142:     for (i=rstart; i<rend; i++) {
143:       for (j=0; j<nzp; j++) {
144:         PetscRandomGetValue(rdm,&rval);
145:         idxn[j] = (PetscInt)(PetscRealPart(rval)*PN);
146:       }
147:       MatSetValues(P,1,&i,nzp,idxn,a,ADD_VALUES);
148:     }
149:     MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
150:     MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);

152:     /* Create R = P^T */
153:     MatTranspose(P,MAT_INITIAL_MATRIX,&R);

155:     { /* Test R = P^T, C1 = R*B */
156:       MatMatMult(R,B,MAT_INITIAL_MATRIX,fill,&C1);
157:       MatTranspose(P,MAT_REUSE_MATRIX,&R);
158:       MatMatMult(R,B,MAT_REUSE_MATRIX,fill,&C1);
159:       MatDestroy(&C1);
160:     }

162:     /* C = P^T*B */
163:     MatTransposeMatMult(P,B,MAT_INITIAL_MATRIX,fill,&C);
164:     MatGetInfo(C,MAT_GLOBAL_SUM,&info);

166:     /* Test MAT_REUSE_MATRIX - reuse symbolic C */
167:     MatTransposeMatMult(P,B,MAT_REUSE_MATRIX,fill,&C);

169:     /* Compare P^T*B and R*B */
170:     MatMatMult(R,B,MAT_INITIAL_MATRIX,fill,&C1);
171:     MatEqual(C,C1,&flg);
172:     if (!flg) {
173:       /* Check norm of C1 = (-1.0)*C + C1 */
174:       PetscReal nrm;
175:       MatAXPY(C1,-1.0,C,DIFFERENT_NONZERO_PATTERN);
176:       MatNorm(C1,NORM_INFINITY,&nrm);
177:       if (nrm > 1.e-14) {
178:         PetscPrintf(PETSC_COMM_WORLD,"Error in MatTransposeMatMult(): %g\n",nrm);
179:       }
180:     }
181:     MatDestroy(&C1);
182:     MatDestroy(&C);

184:     /* C = B*R^T */
185:     if (size == 1) {
186:       MatMatTransposeMult(B,R,MAT_INITIAL_MATRIX,fill,&C);
187:       MatSetOptionsPrefix(C,"matmatmulttr_"); /* enable '-matmatmulttr_' for matrix C */
188:       MatGetInfo(C,MAT_GLOBAL_SUM,&info);

190:       /* Test MAT_REUSE_MATRIX - reuse symbolic C */
191:       MatMatTransposeMult(B,R,MAT_REUSE_MATRIX,fill,&C);

193:       /* Check */
194:       MatMatMult(B,P,MAT_INITIAL_MATRIX,fill,&C1);
195:       MatEqual(C,C1,&flg);
196:       if (!flg) {
197:         PetscPrintf(PETSC_COMM_WORLD,"Error in MatMatTransposeMult()\n");
198:       }
199:       MatDestroy(&C1);
200:       MatDestroy(&C);
201:     }
202:     MatDestroy(&P);
203:     MatDestroy(&R);
204:   }

206:   /* Test MatPtAP() */
207:   /*----------------------*/
208:   if (Test_MatPtAP) {
209:     PetscInt PN;
210:     Mat      Cdup;

212:     MatDuplicate(A_save,MAT_COPY_VALUES,&A);
213:     MatGetSize(A,&M,&N);
214:     MatGetLocalSize(A,&m,&n);
215:     /* PetscPrintf(PETSC_COMM_SELF,"[%d] A: %d,%d, %d,%d\n",rank,m,n,M,N); */

217:     PN   = M/2;
218:     nzp  = (PetscInt)(0.1*PN+1); /* num of nozeros in each row of P */
219:     MatCreate(PETSC_COMM_WORLD,&P);
220:     MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,N,PN);
221:     MatSetType(P,MATAIJ);
222:     MatSeqAIJSetPreallocation(P,nzp,NULL);
223:     MatMPIAIJSetPreallocation(P,nzp,NULL,nzp,NULL);
224:     for (i=0; i<nzp; i++) {
225:       PetscRandomGetValue(rdm,&a[i]);
226:     }
227:     MatGetOwnershipRange(P,&rstart,&rend);
228:     for (i=rstart; i<rend; i++) {
229:       for (j=0; j<nzp; j++) {
230:         PetscRandomGetValue(rdm,&rval);
231:         idxn[j] = (PetscInt)(PetscRealPart(rval)*PN);
232:       }
233:       MatSetValues(P,1,&i,nzp,idxn,a,ADD_VALUES);
234:     }
235:     MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
236:     MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);

238:     /* MatView(P,PETSC_VIEWER_STDOUT_WORLD); */
239:     MatGetSize(P,&pM,&pN);
240:     MatGetLocalSize(P,&pm,&pn);
241:     MatPtAP(A,P,MAT_INITIAL_MATRIX,fill,&C);
242:     /* if (!rank) {PetscPrintf(PETSC_COMM_SELF," MatPtAP() is done, P, %d, %d, %d,%d\n",pm,pn,pM,pN);} */

244:     /* Test MAT_REUSE_MATRIX - reuse symbolic C */
245:     alpha=1.0;
246:     for (i=0; i<2; i++) {
247:       alpha -=0.1;
248:       MatScale(A,alpha);
249:       MatPtAP(A,P,MAT_REUSE_MATRIX,fill,&C);
250:     }

252:     /* Test MatDuplicate() */
253:     MatDuplicate(C,MAT_COPY_VALUES,&Cdup);
254:     MatDestroy(&Cdup);

256:     if (size>1) Test_MatRARt = PETSC_FALSE;
257:     /* Test MatRARt() */
258:     if (Test_MatRARt) {
259:       Mat       R, RARt;
260:       PetscBool equal;
261:       MatTranspose(P,MAT_INITIAL_MATRIX,&R);
262:       MatRARt(A,R,MAT_INITIAL_MATRIX,2.0,&RARt);
263:       MatEqual(C,RARt,&equal);
264:       if (!equal) {
265:         PetscReal norm;
266:         MatAXPY(RARt,-1.0,C,DIFFERENT_NONZERO_PATTERN); /* RARt = -RARt + C */
267:         MatNorm(RARt,NORM_FROBENIUS,&norm);
268:         if (norm > 1.e-12) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"|PtAP - RARt| = %g",norm);
269:       }
270:       MatDestroy(&R);
271:       MatDestroy(&RARt);
272:     }

274:     if (Test_MatMatMatMult && size == 1) {
275:       Mat       R, RAP;
276:       PetscBool equal;
277:       MatTranspose(P,MAT_INITIAL_MATRIX,&R);
278:       MatMatMatMult(R,A,P,MAT_INITIAL_MATRIX,2.0,&RAP);
279:       MatEqual(C,RAP,&equal);
280:       if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"PtAP != RAP");
281:       MatDestroy(&R);
282:       MatDestroy(&RAP);
283:     }

285:     /* Create vector x that is compatible with P */
286:     VecCreate(PETSC_COMM_WORLD,&x);
287:     MatGetLocalSize(P,&m,&n);
288:     VecSetSizes(x,n,PETSC_DECIDE);
289:     VecSetFromOptions(x);

291:     VecCreate(PETSC_COMM_WORLD,&v3);
292:     VecSetSizes(v3,n,PETSC_DECIDE);
293:     VecSetFromOptions(v3);
294:     VecDuplicate(v3,&v4);

296:     norm = 0.0;
297:     for (i=0; i<10; i++) {
298:       VecSetRandom(x,rdm);
299:       MatMult(P,x,v1);
300:       MatMult(A,v1,v2);  /* v2 = A*P*x */

302:       MatMultTranspose(P,v2,v3); /* v3 = Pt*A*P*x */
303:       MatMult(C,x,v4);           /* v3 = C*x   */
304:       VecNorm(v4,NORM_2,&norm_abs);
305:       VecAXPY(v4,none,v3);
306:       VecNorm(v4,NORM_2,&norm_tmp);

308:       norm_tmp /= norm_abs;
309:       if (norm_tmp > norm) norm = norm_tmp;
310:     }
311:     if (norm >= tol) {
312:       PetscPrintf(PETSC_COMM_SELF,"Error: MatPtAP(), |v1 - v2|: %G\n",norm);
313:     }

315:     MatDestroy(&A);
316:     MatDestroy(&P);
317:     MatDestroy(&C);
318:     VecDestroy(&v3);
319:     VecDestroy(&v4);
320:     VecDestroy(&x);
321:   }

323:   /* Destroy objects */
324:   VecDestroy(&v1);
325:   VecDestroy(&v2);
326:   PetscRandomDestroy(&rdm);
327:   PetscFree(idxn);

329:   MatDestroy(&A_save);
330:   MatDestroy(&B);

332:   PetscPreLoadEnd();
333:   PetscFinalize();
334:   return 0;
335: }