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: /* 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,C,C1;
 17:   Vec            x,v1,v2;
 18:   PetscViewer    viewer;
 20:   PetscMPIInt    size,rank;
 21:   PetscInt       i,m,n,j,*idxn,M,N,nzp;
 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_MatMatMultTr=PETSC_TRUE;
 28:   Vec            v3,v4,v5;
 29:   PetscInt       pm,pn,pM,pN;
 30:   PetscBool      Test_MatPtAP=PETSC_TRUE;
 31:   MatInfo        info;
 32: 
 33:   PetscInitialize(&argc,&args,(char *)0,help);
 34:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 35:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

 37:   PetscOptionsGetReal(PETSC_NULL,"-fill",&fill,PETSC_NULL);

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

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

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

 63:   MatGetSize(B,&M,&N);
 64:   nzp  = (PetscInt)(0.1*M);
 65:   PetscMalloc((nzp+1)*(sizeof(PetscInt)+sizeof(PetscScalar)),&idxn);
 66:   a    = (PetscScalar*)(idxn + nzp);
 67: 
 68:   /* Create vectors v1 and v2 that are compatible with A_save */
 69:   VecCreate(PETSC_COMM_WORLD,&v1);
 70:   MatGetLocalSize(A_save,&m,PETSC_NULL);
 71:   VecSetSizes(v1,m,PETSC_DECIDE);
 72:   VecSetFromOptions(v1);
 73:   VecDuplicate(v1,&v2);

 75:   PetscRandomCreate(PETSC_COMM_WORLD,&rdm);
 76:   PetscRandomSetFromOptions(rdm);
 77:   PetscOptionsGetReal(PETSC_NULL,"-fill",&fill,PETSC_NULL);

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

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

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

102:     norm = 0.0;
103:     for (i=0; i<10; i++) {
104:       VecSetRandom(x,rdm);
105:       MatMult(B,x,v1);
106:       MatMult(A,v1,v2);  /* v2 = A*B*x */
107:       MatMult(C,x,v1);   /* v1 = C*x   */
108:       VecNorm(v1,NORM_2,&norm_abs);
109:       VecAXPY(v1,none,v2);
110:       VecNorm(v1,NORM_2,&norm_tmp);
111:       norm_tmp /= norm_abs;
112:       if (norm_tmp > norm) norm = norm_tmp;
113:     }
114:     if (norm >= tol) {
115:       PetscPrintf(PETSC_COMM_SELF,"Error: MatMatMult(), |v1 - v2|: %G\n",norm);
116:     }

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

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

127:   /* Test MatMatMultTranspose() */
128:   /*----------------------------*/
129:   if (size>1) Test_MatMatMultTr = PETSC_FALSE;
130:   if (Test_MatMatMultTr){
131:     PetscInt PN;
132:     PN   = M/2;
133:     nzp  = 5; /* num of nonzeros in each row of P */
134:     MatCreate(PETSC_COMM_WORLD,&P);
135:     MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,M,PN);
136:     MatSetType(P,MATAIJ);
137:     MatSeqAIJSetPreallocation(P,nzp,PETSC_NULL);
138:     MatMPIAIJSetPreallocation(P,nzp,PETSC_NULL,nzp,PETSC_NULL);
139:     for (i=0; i<nzp; i++){
140:       PetscRandomGetValue(rdm,&a[i]);
141:     }
142:     for (i=0; i<M; 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);
151: 
152:     MatMatMultTranspose(P,B,MAT_INITIAL_MATRIX,fill,&C);
153:     MatSetOptionsPrefix(C,"matmatmulttr_"); /* enable '-matmatmulttr_' for matrix C */
154:     MatGetInfo(C,MAT_GLOBAL_SUM,&info);
155:     PetscPrintf(PETSC_COMM_WORLD,"MatMatMultTranspose: nz_allocated = %g; nz_used = %g; nz_unneeded = %g\n",info.nz_allocated,info.nz_used, info.nz_unneeded);

157:     /* Test MAT_REUSE_MATRIX - reuse symbolic C */
158:     alpha=1.0;
159:     for (i=0; i<2; i++){
160:       alpha -=0.1;
161:       MatScale(P,alpha);
162:       MatMatMultTranspose(P,B,MAT_REUSE_MATRIX,fill,&C);
163:     }

165:     /* Create vector x, v5 that are compatible with B */
166:     VecCreate(PETSC_COMM_WORLD,&x);
167:     MatGetLocalSize(B,&m,&n);
168:     VecSetSizes(x,n,PETSC_DECIDE);
169:     VecSetFromOptions(x);

171:     VecCreate(PETSC_COMM_WORLD,&v5);
172:     VecSetSizes(v5,m,PETSC_DECIDE);
173:     VecSetFromOptions(v5);
174: 
175:     MatGetLocalSize(P,PETSC_NULL,&n);
176:     VecCreate(PETSC_COMM_WORLD,&v3);
177:     VecSetSizes(v3,n,PETSC_DECIDE);
178:     VecSetFromOptions(v3);
179:     VecDuplicate(v3,&v4);

181:     norm = 0.0;
182:     for (i=0; i<10; i++) {
183:       VecSetRandom(x,rdm);
184:       MatMult(B,x,v5);            /* v5 = B*x   */
185:       MatMultTranspose(P,v5,v3);  /* v3 = Pt*B*x */
186:       MatMult(C,x,v4);            /* v4 = C*x   */
187:       VecNorm(v4,NORM_2,&norm_abs);
188:       VecAXPY(v4,none,v3);
189:       VecNorm(v4,NORM_2,&norm_tmp);
190:       norm_tmp /= norm_abs;
191:       if (norm_tmp > norm) norm = norm_tmp;
192:     }
193:     if (norm >= tol) {
194:       PetscPrintf(PETSC_COMM_SELF,"Error: MatMatMultTr(), |v3 - v4|: %G\n",norm);
195:     }
196:     MatDestroy(&P);
197:     MatDestroy(&C);
198:     VecDestroy(&v3);
199:     VecDestroy(&v4);
200:     VecDestroy(&v5);
201:     VecDestroy(&x);
202:   }

204:   /* Test MatPtAP() */
205:   /*----------------------*/
206:   if (Test_MatPtAP){
207:     PetscInt PN;
208:     MatDuplicate(A_save,MAT_COPY_VALUES,&A);
209:     MatGetSize(A,&M,&N);
210:     MatGetLocalSize(A,&m,&n);
211:     /* PetscPrintf(PETSC_COMM_SELF,"[%d] A: %d,%d, %d,%d\n",rank,m,n,M,N); */

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

233:     /* MatView(P,PETSC_VIEWER_STDOUT_WORLD); */
234:     MatGetSize(P,&pM,&pN);
235:     MatGetLocalSize(P,&pm,&pn);
236:     /* PetscPrintf(PETSC_COMM_SELF," [%d] P, %d, %d, %d,%d\n",rank,pm,pn,pM,pN); */
237:     MatPtAP(A,P,MAT_INITIAL_MATRIX,fill,&C);

239:     /* Test MAT_REUSE_MATRIX - reuse symbolic C */
240:     alpha=1.0;
241:     for (i=0; i<2; i++){
242:       alpha -=0.1;
243:       MatScale(A,alpha);
244:       MatPtAP(A,P,MAT_REUSE_MATRIX,fill,&C);
245:     }

247:     /* Create vector x that is compatible with P */
248:     VecCreate(PETSC_COMM_WORLD,&x);
249:     MatGetLocalSize(P,&m,&n);
250:     VecSetSizes(x,n,PETSC_DECIDE);
251:     VecSetFromOptions(x);
252: 
253:     VecCreate(PETSC_COMM_WORLD,&v3);
254:     VecSetSizes(v3,n,PETSC_DECIDE);
255:     VecSetFromOptions(v3);
256:     VecDuplicate(v3,&v4);

258:     norm = 0.0;
259:     for (i=0; i<10; i++) {
260:       VecSetRandom(x,rdm);
261:       MatMult(P,x,v1);
262:       MatMult(A,v1,v2);  /* v2 = A*P*x */

264:       MatMultTranspose(P,v2,v3); /* v3 = Pt*A*P*x */
265:       MatMult(C,x,v4);           /* v3 = C*x   */
266:       VecNorm(v4,NORM_2,&norm_abs);
267:       VecAXPY(v4,none,v3);
268:       VecNorm(v4,NORM_2,&norm_tmp);
269:       norm_tmp /= norm_abs;
270:       if (norm_tmp > norm) norm = norm_tmp;
271:     }
272:     if (norm >= tol) {
273:       PetscPrintf(PETSC_COMM_SELF,"Error: MatPtAP(), |v1 - v2|: %G\n",norm);
274:     }
275: 
276:     MatDestroy(&A);
277:     MatDestroy(&P);
278:     MatDestroy(&C);
279:     VecDestroy(&v3);
280:     VecDestroy(&v4);
281:     VecDestroy(&x);
282:   } /* if (Test_MatPtAP) */

284:   /* Destroy objects */
285:   VecDestroy(&v1);
286:   VecDestroy(&v2);
287:   PetscRandomDestroy(&rdm);
288:   PetscFree(idxn);
289: 
290:   MatDestroy(&A_save);
291:   MatDestroy(&B);

293:   PetscPreLoadEnd();
294:   PetscFinalize();

296:   return 0;
297: }