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*/