Actual source code: ex94.c

petsc-3.3-p6 2013-02-11
  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_info -matmatmulttr_mat_view_info
  7:    mpiexec -n 3 ./ex94 -f0 medium -f1 medium -f2 arco1 -f3 arco1 -matmatmult_mat_view_info
  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;
 28:   PetscInt       pm,pn,pM,pN;
 29:   MatInfo        info;
 30: 
 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(PETSC_NULL,"-fill",&fill,PETSC_NULL);

 37:   /*  Load the matrices A_save and B */
 38:   PetscOptionsGetString(PETSC_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(PETSC_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(PETSC_NULL,"-f2",file[2],128,&flg);
 43:   if (!flg) {
 44:     preload = PETSC_FALSE;
 45:   } else {
 46:     PetscOptionsGetString(PETSC_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  = (PetscInt)(0.1*M);
 63:   PetscMalloc((nzp+1)*(sizeof(PetscInt)+sizeof(PetscScalar)),&idxn);
 64:   a    = (PetscScalar*)(idxn + nzp);
 65: 
 66:   /* Create vectors v1 and v2 that are compatible with A_save */
 67:   VecCreate(PETSC_COMM_WORLD,&v1);
 68:   MatGetLocalSize(A_save,&m,PETSC_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(PETSC_NULL,"-fill",&fill,PETSC_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,PETSC_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);
109:       norm_tmp /= norm_abs;
110:       if (norm_tmp > norm) norm = norm_tmp;
111:     }
112:     if (norm >= tol) {
113:       PetscPrintf(PETSC_COMM_SELF,"Error: MatMatMult(), |v1 - v2|: %G\n",norm);
114:     }

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

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

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

151:     /* Create R = P^T */
152:     MatTranspose(P,MAT_INITIAL_MATRIX,&R);
153: 
154:     /* C = P^T*B */
155:     MatTransposeMatMult(P,B,MAT_INITIAL_MATRIX,fill,&C);
156:     MatGetInfo(C,MAT_GLOBAL_SUM,&info);
157: 
158:     /* Test MAT_REUSE_MATRIX - reuse symbolic C */
159:     MatTransposeMatMult(P,B,MAT_REUSE_MATRIX,fill,&C);
160: 
161:     /* Compare P^T*B and R*B */
162:     MatMatMult(R,B,MAT_INITIAL_MATRIX,fill,&C1);
163:     MatEqual(C,C1,&flg);
164:     if (!flg){
165:       /* Check norm of C1 = (-1.0)*C + C1 */
166:       PetscReal nrm;
167:       MatAXPY(C1,-1.0,C,DIFFERENT_NONZERO_PATTERN);
168:       MatNorm(C1,NORM_INFINITY,&nrm);
169:       if (nrm > 1.e-14){
170:         PetscPrintf(PETSC_COMM_WORLD,"Error in MatTransposeMatMult(): %g\n",nrm);
171:       }
172:     }
173:     MatDestroy(&C1);
174:     MatDestroy(&C);

176:     /* C = B*R^T */
177:     if (size == 1){
178:       MatMatTransposeMult(B,R,MAT_INITIAL_MATRIX,fill,&C);
179:       MatSetOptionsPrefix(C,"matmatmulttr_"); /* enable '-matmatmulttr_' for matrix C */
180:       MatGetInfo(C,MAT_GLOBAL_SUM,&info);
181: 
182:       /* Test MAT_REUSE_MATRIX - reuse symbolic C */
183:       MatMatTransposeMult(B,R,MAT_REUSE_MATRIX,fill,&C);
184: 
185:       /* Check */
186:       MatMatMult(B,P,MAT_INITIAL_MATRIX,fill,&C1);
187:       MatEqual(C,C1,&flg);
188:       if (!flg){
189:         PetscPrintf(PETSC_COMM_WORLD,"Error in MatMatTransposeMult()\n");
190:       }
191:       MatDestroy(&C1);
192:       MatDestroy(&C);
193:     }
194:     MatDestroy(&P);
195:     MatDestroy(&R);
196:   }

198:   /* Test MatPtAP() */
199:   /*----------------------*/
200:   if (Test_MatPtAP){
201:     PetscInt PN;
202:     MatDuplicate(A_save,MAT_COPY_VALUES,&A);
203:     MatGetSize(A,&M,&N);
204:     MatGetLocalSize(A,&m,&n);
205:     /* PetscPrintf(PETSC_COMM_SELF,"[%d] A: %d,%d, %d,%d\n",rank,m,n,M,N); */

207:     PN   = M/2;
208:     nzp  = (PetscInt)(0.1*PN+1); /* num of nozeros in each row of P */
209:     MatCreate(PETSC_COMM_WORLD,&P);
210:     MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,N,PN);
211:     MatSetType(P,MATAIJ);
212:     MatSeqAIJSetPreallocation(P,nzp,PETSC_NULL);
213:     MatMPIAIJSetPreallocation(P,nzp,PETSC_NULL,nzp,PETSC_NULL);
214:     for (i=0; i<nzp; i++){
215:       PetscRandomGetValue(rdm,&a[i]);
216:     }
217:     MatGetOwnershipRange(P,&rstart,&rend);
218:     for (i=rstart; i<rend; i++){
219:       for (j=0; j<nzp; j++){
220:         PetscRandomGetValue(rdm,&rval);
221:         idxn[j] = (PetscInt)(PetscRealPart(rval)*PN);
222:       }
223:       MatSetValues(P,1,&i,nzp,idxn,a,ADD_VALUES);
224:     }
225:     MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
226:     MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);

228:     /* MatView(P,PETSC_VIEWER_STDOUT_WORLD); */
229:     MatGetSize(P,&pM,&pN);
230:     MatGetLocalSize(P,&pm,&pn);
231:     MatPtAP(A,P,MAT_INITIAL_MATRIX,fill,&C);
232:     /* if (!rank){PetscPrintf(PETSC_COMM_SELF," MatPtAP() is done, P, %d, %d, %d,%d\n",pm,pn,pM,pN);} */

234:     /* Test MAT_REUSE_MATRIX - reuse symbolic C */
235:     alpha=1.0;
236:     for (i=0; i<2; i++){
237:       alpha -=0.1;
238:       MatScale(A,alpha);
239:       MatPtAP(A,P,MAT_REUSE_MATRIX,fill,&C);
240:     }

242:     /* Test MatDuplicate() */
243:     Mat Cdup;
244:     MatDuplicate(C,MAT_COPY_VALUES,&Cdup);
245:     MatDestroy(&Cdup);

247:     if (size>1) Test_MatRARt = PETSC_FALSE;
248:     /* Test MatRARt() */
249:     if (Test_MatRARt){
250:       Mat       R, RARt;
251:       PetscBool equal;
252:       MatTranspose(P,MAT_INITIAL_MATRIX,&R);
253:       MatRARt(A,R,MAT_INITIAL_MATRIX,2.0,&RARt);
254:       MatEqual(C,RARt,&equal);
255:       if (!equal) {
256:         PetscReal norm;
257:         MatAXPY(RARt,-1.0,C,DIFFERENT_NONZERO_PATTERN); /* RARt = -RARt + C */
258:         MatNorm(RARt,NORM_FROBENIUS,&norm);
259:         if (norm > 1.e-12) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"|PtAP - RARt| = %g",norm);
260:       }
261:       MatDestroy(&R);
262:       MatDestroy(&RARt);
263:     }

265:     /* Create vector x that is compatible with P */
266:     VecCreate(PETSC_COMM_WORLD,&x);
267:     MatGetLocalSize(P,&m,&n);
268:     VecSetSizes(x,n,PETSC_DECIDE);
269:     VecSetFromOptions(x);
270: 
271:     VecCreate(PETSC_COMM_WORLD,&v3);
272:     VecSetSizes(v3,n,PETSC_DECIDE);
273:     VecSetFromOptions(v3);
274:     VecDuplicate(v3,&v4);

276:     norm = 0.0;
277:     for (i=0; i<10; i++) {
278:       VecSetRandom(x,rdm);
279:       MatMult(P,x,v1);
280:       MatMult(A,v1,v2);  /* v2 = A*P*x */

282:       MatMultTranspose(P,v2,v3); /* v3 = Pt*A*P*x */
283:       MatMult(C,x,v4);           /* v3 = C*x   */
284:       VecNorm(v4,NORM_2,&norm_abs);
285:       VecAXPY(v4,none,v3);
286:       VecNorm(v4,NORM_2,&norm_tmp);
287:       norm_tmp /= norm_abs;
288:       if (norm_tmp > norm) norm = norm_tmp;
289:     }
290:     if (norm >= tol) {
291:       PetscPrintf(PETSC_COMM_SELF,"Error: MatPtAP(), |v1 - v2|: %G\n",norm);
292:     }
293: 
294:     MatDestroy(&A);
295:     MatDestroy(&P);
296:     MatDestroy(&C);
297:     VecDestroy(&v3);
298:     VecDestroy(&v4);
299:     VecDestroy(&x);
300:   }

302:   /* Destroy objects */
303:   VecDestroy(&v1);
304:   VecDestroy(&v2);
305:   PetscRandomDestroy(&rdm);
306:   PetscFree(idxn);
307: 
308:   MatDestroy(&A_save);
309:   MatDestroy(&B);

311:   PetscPreLoadEnd();
312:   PetscFinalize();
313:   return 0;
314: }