Actual source code: ex2.c


  2: static char help[] = "Tests MatTranspose(), MatNorm(), MatAXPY() and MatAYPX().\n\n";

  4: #include <petscmat.h>

  6: static PetscErrorCode TransposeAXPY(Mat C,PetscScalar alpha,Mat mat,PetscErrorCode (*f)(Mat,Mat*))
  7: {
  8:   Mat            D,E,F,G;
  9:   MatType        mtype;

 11:   MatGetType(mat,&mtype);
 12:   if (f == MatCreateTranspose) {
 13:     PetscPrintf(PETSC_COMM_WORLD,"\nMatAXPY:  (C^T)^T = (C^T)^T + alpha * A, C=A, SAME_NONZERO_PATTERN\n");
 14:   } else {
 15:     PetscPrintf(PETSC_COMM_WORLD,"\nMatAXPY:  (C^H)^H = (C^H)^H + alpha * A, C=A, SAME_NONZERO_PATTERN\n");
 16:   }
 17:   MatDuplicate(mat,MAT_COPY_VALUES,&C);
 18:   f(C,&D);
 19:   f(D,&E);
 20:   MatAXPY(E,alpha,mat,SAME_NONZERO_PATTERN);
 21:   MatConvert(E,mtype,MAT_INPLACE_MATRIX,&E);
 22:   MatView(E,PETSC_VIEWER_STDOUT_WORLD);
 23:   MatDestroy(&E);
 24:   MatDestroy(&D);
 25:   MatDestroy(&C);
 26:   if (f == MatCreateTranspose) {
 27:     PetscPrintf(PETSC_COMM_WORLD,"MatAXPY:  C = C + alpha * (A^T)^T, C=A, SAME_NONZERO_PATTERN\n");
 28:   } else {
 29:     PetscPrintf(PETSC_COMM_WORLD,"MatAXPY:  C = C + alpha * (A^H)^H, C=A, SAME_NONZERO_PATTERN\n");
 30:   }
 31:   if (f == MatCreateTranspose) {
 32:     MatTranspose(mat,MAT_INITIAL_MATRIX,&D);
 33:   } else {
 34:     MatHermitianTranspose(mat,MAT_INITIAL_MATRIX,&D);
 35:   }
 36:   f(D,&E);
 37:   MatDuplicate(mat,MAT_COPY_VALUES,&C);
 38:   MatAXPY(C,alpha,E,SAME_NONZERO_PATTERN);
 39:   MatView(C,PETSC_VIEWER_STDOUT_WORLD);
 40:   MatDestroy(&E);
 41:   MatDestroy(&D);
 42:   MatDestroy(&C);
 43:   if (f == MatCreateTranspose) {
 44:     PetscPrintf(PETSC_COMM_WORLD,"MatAXPY:  (C^T)^T = (C^T)^T + alpha * (A^T)^T, C=A, SAME_NONZERO_PATTERN\n");
 45:   } else {
 46:     PetscPrintf(PETSC_COMM_WORLD,"MatAXPY:  (C^H)^H = (C^H)^H + alpha * (A^H)^H, C=A, SAME_NONZERO_PATTERN\n");
 47:   }
 48:   MatDuplicate(mat,MAT_COPY_VALUES,&C);
 49:   f(C,&D);
 50:   f(D,&E);
 51:   f(mat,&F);
 52:   f(F,&G);
 53:   MatAXPY(E,alpha,G,SAME_NONZERO_PATTERN);
 54:   MatConvert(E,mtype,MAT_INPLACE_MATRIX,&E);
 55:   MatView(E,PETSC_VIEWER_STDOUT_WORLD);
 56:   MatDestroy(&G);
 57:   MatDestroy(&F);
 58:   MatDestroy(&E);
 59:   MatDestroy(&D);
 60:   MatDestroy(&C);

 62:   /* Call f on a matrix that does not implement the transposition */
 63:   PetscPrintf(PETSC_COMM_WORLD,"MatAXPY:  Now without the transposition operation\n");
 64:   MatConvert(mat,MATSHELL,MAT_INITIAL_MATRIX,&C);
 65:   f(C,&D);
 66:   f(D,&E);
 67:   /* XXX cannot use MAT_INPLACE_MATRIX, it leaks mat */
 68:   MatConvert(E,mtype,MAT_INITIAL_MATRIX,&F);
 69:   MatAXPY(F,alpha,mat,SAME_NONZERO_PATTERN);
 70:   MatView(F,PETSC_VIEWER_STDOUT_WORLD);
 71:   MatDestroy(&F);
 72:   MatDestroy(&E);
 73:   MatDestroy(&D);
 74:   MatDestroy(&C);
 75:   return 0;
 76: }

 78: int main(int argc,char **argv)
 79: {
 80:   Mat            mat,tmat = 0;
 81:   PetscInt       m = 7,n,i,j,rstart,rend,rect = 0;
 82:   PetscMPIInt    size,rank;
 83:   PetscBool      flg;
 84:   PetscScalar    v, alpha;
 85:   PetscReal      normf,normi,norm1;

 87:   PetscInitialize(&argc,&argv,(char*)0,help);
 88:   PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_COMMON);
 89:   PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
 90:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 91:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 92:   n    = m;
 93:   PetscOptionsHasName(NULL,NULL,"-rectA",&flg);
 94:   if (flg) {n += 2; rect = 1;}
 95:   PetscOptionsHasName(NULL,NULL,"-rectB",&flg);
 96:   if (flg) {n -= 2; rect = 1;}

 98:   /* ------- Assemble matrix --------- */
 99:   MatCreate(PETSC_COMM_WORLD,&mat);
100:   MatSetSizes(mat,PETSC_DECIDE,PETSC_DECIDE,m,n);
101:   MatSetFromOptions(mat);
102:   MatSetUp(mat);
103:   MatGetOwnershipRange(mat,&rstart,&rend);
104:   for (i=rstart; i<rend; i++) {
105:     for (j=0; j<n; j++) {
106:       v    = 10.0*i+j+1.0;
107:       MatSetValues(mat,1,&i,1,&j,&v,INSERT_VALUES);
108:     }
109:   }
110:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
111:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);

113:   /* ----------------- Test MatNorm()  ----------------- */
114:   MatNorm(mat,NORM_FROBENIUS,&normf);
115:   MatNorm(mat,NORM_1,&norm1);
116:   MatNorm(mat,NORM_INFINITY,&normi);
117:   PetscPrintf(PETSC_COMM_WORLD,"original A: Frobenious norm = %g, one norm = %g, infinity norm = %g\n",(double)normf,(double)norm1,(double)normi);
118:   MatView(mat,PETSC_VIEWER_STDOUT_WORLD);

120:   /* --------------- Test MatTranspose()  -------------- */
121:   PetscOptionsHasName(NULL,NULL,"-in_place",&flg);
122:   if (!rect && flg) {
123:     MatTranspose(mat,MAT_REUSE_MATRIX,&mat);   /* in-place transpose */
124:     tmat = mat;
125:     mat  = NULL;
126:   } else { /* out-of-place transpose */
127:     MatTranspose(mat,MAT_INITIAL_MATRIX,&tmat);
128:   }

130:   /* ----------------- Test MatNorm()  ----------------- */
131:   /* Print info about transpose matrix */
132:   MatNorm(tmat,NORM_FROBENIUS,&normf);
133:   MatNorm(tmat,NORM_1,&norm1);
134:   MatNorm(tmat,NORM_INFINITY,&normi);
135:   PetscPrintf(PETSC_COMM_WORLD,"B = A^T: Frobenious norm = %g, one norm = %g, infinity norm = %g\n",(double)normf,(double)norm1,(double)normi);
136:   MatView(tmat,PETSC_VIEWER_STDOUT_WORLD);

138:   /* ----------------- Test MatAXPY(), MatAYPX()  ----------------- */
139:   if (mat && !rect) {
140:     alpha = 1.0;
141:     PetscOptionsGetScalar(NULL,NULL,"-alpha",&alpha,NULL);
142:     PetscPrintf(PETSC_COMM_WORLD,"MatAXPY:  B = B + alpha * A\n");
143:     MatAXPY(tmat,alpha,mat,DIFFERENT_NONZERO_PATTERN);
144:     MatView(tmat,PETSC_VIEWER_STDOUT_WORLD);

146:     PetscPrintf(PETSC_COMM_WORLD,"MatAYPX:  B = alpha*B + A\n");
147:     MatAYPX(tmat,alpha,mat,DIFFERENT_NONZERO_PATTERN);
148:     MatView(tmat,PETSC_VIEWER_STDOUT_WORLD);
149:   }

151:   {
152:     Mat C;
153:     alpha = 1.0;
154:     PetscPrintf(PETSC_COMM_WORLD,"MatAXPY:  C = C + alpha * A, C=A, SAME_NONZERO_PATTERN\n");
155:     MatDuplicate(mat,MAT_COPY_VALUES,&C);
156:     MatAXPY(C,alpha,mat,SAME_NONZERO_PATTERN);
157:     MatView(C,PETSC_VIEWER_STDOUT_WORLD);
158:     MatDestroy(&C);
159:     TransposeAXPY(C,alpha,mat,MatCreateTranspose);
160:     TransposeAXPY(C,alpha,mat,MatCreateHermitianTranspose);
161:   }

163:   {
164:     Mat matB;
165:     /* get matB that has nonzeros of mat in all even numbers of row and col */
166:     MatCreate(PETSC_COMM_WORLD,&matB);
167:     MatSetSizes(matB,PETSC_DECIDE,PETSC_DECIDE,m,n);
168:     MatSetFromOptions(matB);
169:     MatSetUp(matB);
170:     MatGetOwnershipRange(matB,&rstart,&rend);
171:     if (rstart % 2 != 0) rstart++;
172:     for (i=rstart; i<rend; i += 2) {
173:       for (j=0; j<n; j += 2) {
174:         v    = 10.0*i+j+1.0;
175:         MatSetValues(matB,1,&i,1,&j,&v,INSERT_VALUES);
176:       }
177:     }
178:     MatAssemblyBegin(matB,MAT_FINAL_ASSEMBLY);
179:     MatAssemblyEnd(matB,MAT_FINAL_ASSEMBLY);
180:     PetscPrintf(PETSC_COMM_WORLD," A: original matrix:\n");
181:     MatView(mat,PETSC_VIEWER_STDOUT_WORLD);
182:     PetscPrintf(PETSC_COMM_WORLD," B(a subset of A):\n");
183:     MatView(matB,PETSC_VIEWER_STDOUT_WORLD);
184:     PetscPrintf(PETSC_COMM_WORLD,"MatAXPY:  B = B + alpha * A, SUBSET_NONZERO_PATTERN\n");
185:     MatAXPY(mat,alpha,matB,SUBSET_NONZERO_PATTERN);
186:     MatView(mat,PETSC_VIEWER_STDOUT_WORLD);
187:     MatDestroy(&matB);
188:   }

190:   /* Test MatZeroRows */
191:   j = rstart - 1;
192:   if (j < 0) j = m-1;
193:   MatZeroRows(mat,1,&j,0.0,NULL,NULL);
194:   MatView(mat,PETSC_VIEWER_STDOUT_WORLD);

196:   PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
197:   /* Free data structures */
198:   MatDestroy(&mat);
199:   MatDestroy(&tmat);
200:   PetscFinalize();
201:   return 0;
202: }

204: /*TEST

206:    test:
207:       suffix: 11_A
208:       args: -mat_type seqaij -rectA
209:       filter: grep -v "Mat Object"

211:    test:
212:       suffix: 12_A
213:       args: -mat_type seqdense -rectA
214:       filter: grep -v type | grep -v "Mat Object"

216:    test:
217:       requires: cuda
218:       suffix: 12_A_cuda
219:       args: -mat_type seqdensecuda -rectA
220:       output_file: output/ex2_12_A.out
221:       filter: grep -v type | grep -v "Mat Object"

223:    test:
224:       requires: kokkos_kernels
225:       suffix: 12_A_kokkos
226:       args: -mat_type aijkokkos -rectA
227:       output_file: output/ex2_12_A.out
228:       filter: grep -v type | grep -v "Mat Object"

230:    test:
231:       suffix: 11_B
232:       args: -mat_type seqaij -rectB
233:       filter: grep -v "Mat Object"

235:    test:
236:       suffix: 12_B
237:       args: -mat_type seqdense -rectB
238:       filter: grep -v type | grep -v "Mat Object"

240:    test:
241:       requires: cuda
242:       suffix: 12_B_cuda
243:       args: -mat_type seqdensecuda -rectB
244:       output_file: output/ex2_12_B.out
245:       filter: grep -v type | grep -v "Mat Object"

247:    test:
248:       requires: kokkos_kernels
249:       suffix: 12_B_kokkos
250:       args: -mat_type aijkokkos -rectB
251:       output_file: output/ex2_12_B.out
252:       filter: grep -v type | grep -v "Mat Object"

254:    test:
255:       suffix: 21
256:       args: -mat_type mpiaij
257:       filter: grep -v type | grep -v "MPI processes"

259:    test:
260:       suffix: 22
261:       args: -mat_type mpidense
262:       filter: grep -v type | grep -v "Mat Object"

264:    test:
265:       requires: cuda
266:       suffix: 22_cuda
267:       output_file: output/ex2_22.out
268:       args: -mat_type mpidensecuda
269:       filter: grep -v type | grep -v "Mat Object"

271:    test:
272:       requires: kokkos_kernels
273:       suffix: 22_kokkos
274:       output_file: output/ex2_22.out
275:       args: -mat_type aijkokkos
276:       filter: grep -v type | grep -v "Mat Object"

278:    test:
279:       suffix: 23
280:       nsize: 3
281:       args: -mat_type mpiaij
282:       filter: grep -v type | grep -v "MPI processes"

284:    test:
285:       suffix: 24
286:       nsize: 3
287:       args: -mat_type mpidense
288:       filter: grep -v type | grep -v "Mat Object"

290:    test:
291:       requires: cuda
292:       suffix: 24_cuda
293:       nsize: 3
294:       output_file: output/ex2_24.out
295:       args: -mat_type mpidensecuda
296:       filter: grep -v type | grep -v "Mat Object"

298:    test:
299:       suffix: 2_aijcusparse_1
300:       args: -mat_type mpiaijcusparse
301:       output_file: output/ex2_21.out
302:       requires: cuda
303:       filter: grep -v type | grep -v "MPI processes"

305:    test:
306:       suffix: 2_aijkokkos_1
307:       args: -mat_type aijkokkos
308:       output_file: output/ex2_21.out
309:       requires: kokkos_kernels
310:       filter: grep -v type | grep -v "MPI processes"

312:    test:
313:       suffix: 2_aijcusparse_2
314:       nsize: 3
315:       args: -mat_type mpiaijcusparse
316:       output_file: output/ex2_23.out
317:       requires: cuda
318:       filter: grep -v type | grep -v "MPI processes"

320:    test:
321:       suffix: 2_aijkokkos_2
322:       nsize: 3
323:       args: -mat_type aijkokkos
324:       output_file: output/ex2_23.out
325:       # Turn off hip due to intermittent CI failures on hip.txcorp.com. Should re-enable this test when the machine is upgraded.
326:       requires: !sycl !hip kokkos_kernels
327:       filter: grep -v type | grep -v "MPI processes"

329:    test:
330:       suffix: 3
331:       nsize: 2
332:       args: -mat_type mpiaij -rectA

334:    test:
335:       suffix: 3_aijcusparse
336:       nsize: 2
337:       args: -mat_type mpiaijcusparse -rectA
338:       requires: cuda

340:    test:
341:       suffix: 4
342:       nsize: 2
343:       args: -mat_type mpidense -rectA
344:       filter: grep -v type | grep -v "MPI processes"

346:    test:
347:       requires: cuda
348:       suffix: 4_cuda
349:       nsize: 2
350:       output_file: output/ex2_4.out
351:       args: -mat_type mpidensecuda -rectA
352:       filter: grep -v type | grep -v "MPI processes"

354:    test:
355:       suffix: aijcusparse_1
356:       args: -mat_type seqaijcusparse -rectA
357:       filter: grep -v "Mat Object"
358:       output_file: output/ex2_11_A_aijcusparse.out
359:       requires: cuda

361:    test:
362:       suffix: aijcusparse_2
363:       args: -mat_type seqaijcusparse -rectB
364:       filter: grep -v "Mat Object"
365:       output_file: output/ex2_11_B_aijcusparse.out
366:       requires: cuda

368: TEST*/