Actual source code: ex2.c
petsc-3.14.6 2021-03-30
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;
12: if (f == MatCreateTranspose) {
13: PetscPrintf(PETSC_COMM_WORLD,"MatAXPY: (C^T)^T = (C^T)^T + alpha * A, C=A, SAME_NONZERO_PATTERN\n");
14: } else {
15: PetscPrintf(PETSC_COMM_WORLD,"MatAXPY: (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: MatView(C,PETSC_VIEWER_STDOUT_WORLD);
22: MatDestroy(&E);
23: MatDestroy(&D);
24: MatDestroy(&C);
25: if (f == MatCreateTranspose) {
26: PetscPrintf(PETSC_COMM_WORLD,"MatAXPY: C = C + alpha * (A^T)^T, C=A, SAME_NONZERO_PATTERN\n");
27: } else {
28: PetscPrintf(PETSC_COMM_WORLD,"MatAXPY: C = C + alpha * (A^H)^H, C=A, SAME_NONZERO_PATTERN\n");
29: }
30: MatDuplicate(mat,MAT_COPY_VALUES,&C);
31: /* MATTRANSPOSE should have a MatTranspose_Transpose or MatTranspose_HT implementation */
32: if (f == MatCreateTranspose) {
33: MatTranspose(mat,MAT_INITIAL_MATRIX,&D);
34: } else {
35: MatHermitianTranspose(mat,MAT_INITIAL_MATRIX,&D);
36: }
37: f(D,&E);
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: MatView(C,PETSC_VIEWER_STDOUT_WORLD);
55: MatDestroy(&G);
56: MatDestroy(&F);
57: MatDestroy(&E);
58: MatDestroy(&D);
59: MatDestroy(&C);
60: return(0);
61: }
63: int main(int argc,char **argv)
64: {
65: Mat mat,tmat = 0;
66: PetscInt m = 7,n,i,j,rstart,rend,rect = 0;
68: PetscMPIInt size,rank;
69: PetscBool flg;
70: PetscScalar v, alpha;
71: PetscReal normf,normi,norm1;
73: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
74: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_COMMON);
75: PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
76: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
77: MPI_Comm_size(PETSC_COMM_WORLD,&size);
78: n = m;
79: PetscOptionsHasName(NULL,NULL,"-rectA",&flg);
80: if (flg) {n += 2; rect = 1;}
81: PetscOptionsHasName(NULL,NULL,"-rectB",&flg);
82: if (flg) {n -= 2; rect = 1;}
84: /* ------- Assemble matrix --------- */
85: MatCreate(PETSC_COMM_WORLD,&mat);
86: MatSetSizes(mat,PETSC_DECIDE,PETSC_DECIDE,m,n);
87: MatSetFromOptions(mat);
88: MatSetUp(mat);
89: MatGetOwnershipRange(mat,&rstart,&rend);
90: for (i=rstart; i<rend; i++) {
91: for (j=0; j<n; j++) {
92: v = 10.0*i+j;
93: MatSetValues(mat,1,&i,1,&j,&v,INSERT_VALUES);
94: }
95: }
96: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
97: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
99: /* ----------------- Test MatNorm() ----------------- */
100: MatNorm(mat,NORM_FROBENIUS,&normf);
101: MatNorm(mat,NORM_1,&norm1);
102: MatNorm(mat,NORM_INFINITY,&normi);
103: PetscPrintf(PETSC_COMM_WORLD,"original A: Frobenious norm = %g, one norm = %g, infinity norm = %g\n",(double)normf,(double)norm1,(double)normi);
104: MatView(mat,PETSC_VIEWER_STDOUT_WORLD);
106: /* --------------- Test MatTranspose() -------------- */
107: PetscOptionsHasName(NULL,NULL,"-in_place",&flg);
108: if (!rect && flg) {
109: MatTranspose(mat,MAT_REUSE_MATRIX,&mat); /* in-place transpose */
110: tmat = mat; mat = 0;
111: } else { /* out-of-place transpose */
112: MatTranspose(mat,MAT_INITIAL_MATRIX,&tmat);
113: }
115: /* ----------------- Test MatNorm() ----------------- */
116: /* Print info about transpose matrix */
117: MatNorm(tmat,NORM_FROBENIUS,&normf);
118: MatNorm(tmat,NORM_1,&norm1);
119: MatNorm(tmat,NORM_INFINITY,&normi);
120: PetscPrintf(PETSC_COMM_WORLD,"B = A^T: Frobenious norm = %g, one norm = %g, infinity norm = %g\n",(double)normf,(double)norm1,(double)normi);
121: MatView(tmat,PETSC_VIEWER_STDOUT_WORLD);
123: /* ----------------- Test MatAXPY(), MatAYPX() ----------------- */
124: if (mat && !rect) {
125: alpha = 1.0;
126: PetscOptionsGetScalar(NULL,NULL,"-alpha",&alpha,NULL);
127: PetscPrintf(PETSC_COMM_WORLD,"MatAXPY: B = B + alpha * A\n");
128: MatAXPY(tmat,alpha,mat,DIFFERENT_NONZERO_PATTERN);
129: MatView(tmat,PETSC_VIEWER_STDOUT_WORLD);
131: PetscPrintf(PETSC_COMM_WORLD,"MatAYPX: B = alpha*B + A\n");
132: MatAYPX(tmat,alpha,mat,DIFFERENT_NONZERO_PATTERN);
133: MatView(tmat,PETSC_VIEWER_STDOUT_WORLD);
134: }
136: {
137: Mat C;
138: alpha = 1.0;
139: PetscPrintf(PETSC_COMM_WORLD,"MatAXPY: C = C + alpha * A, C=A, SAME_NONZERO_PATTERN\n");
140: MatDuplicate(mat,MAT_COPY_VALUES,&C);
141: MatAXPY(C,alpha,mat,SAME_NONZERO_PATTERN);
142: MatView(C,PETSC_VIEWER_STDOUT_WORLD);
143: MatDestroy(&C);
144: TransposeAXPY(C,alpha,mat,MatCreateTranspose);
145: TransposeAXPY(C,alpha,mat,MatCreateHermitianTranspose);
146: }
148: {
149: Mat matB;
150: /* get matB that has nonzeros of mat in all even numbers of row and col */
151: MatCreate(PETSC_COMM_WORLD,&matB);
152: MatSetSizes(matB,PETSC_DECIDE,PETSC_DECIDE,m,n);
153: MatSetFromOptions(matB);
154: MatSetUp(matB);
155: MatGetOwnershipRange(matB,&rstart,&rend);
156: if (rstart % 2 != 0) rstart++;
157: for (i=rstart; i<rend; i += 2) {
158: for (j=0; j<n; j += 2) {
159: v = 10.0*i+j;
160: MatSetValues(matB,1,&i,1,&j,&v,INSERT_VALUES);
161: }
162: }
163: MatAssemblyBegin(matB,MAT_FINAL_ASSEMBLY);
164: MatAssemblyEnd(matB,MAT_FINAL_ASSEMBLY);
165: PetscPrintf(PETSC_COMM_WORLD," A: original matrix:\n");
166: MatView(mat,PETSC_VIEWER_STDOUT_WORLD);
167: PetscPrintf(PETSC_COMM_WORLD," B(a subset of A):\n");
168: MatView(matB,PETSC_VIEWER_STDOUT_WORLD);
169: PetscPrintf(PETSC_COMM_WORLD,"MatAXPY: B = B + alpha * A, SUBSET_NONZERO_PATTERN\n");
170: MatAXPY(mat,alpha,matB,SUBSET_NONZERO_PATTERN);
171: MatView(mat,PETSC_VIEWER_STDOUT_WORLD);
172: MatDestroy(&matB);
173: }
175: /* Test MatZeroRows */
176: j = rstart - 1;
177: if (j < 0) j = m-1;
178: MatZeroRows(mat,1,&j,0.0,NULL,NULL);
179: MatView(mat,PETSC_VIEWER_STDOUT_WORLD);
181: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
182: /* Free data structures */
183: MatDestroy(&mat);
184: MatDestroy(&tmat);
185: PetscFinalize();
186: return ierr;
187: }
189: /*TEST
191: test:
192: suffix: 11_A
193: args: -mat_type seqaij -rectA
194: filter: grep -v "Mat Object"
196: test:
197: suffix: 12_A
198: args: -mat_type seqdense -rectA
199: filter: grep -v type | grep -v "Mat Object"
201: test:
202: requires: cuda
203: suffix: 12_A_cuda
204: args: -mat_type seqdensecuda -rectA
205: output_file: output/ex2_12_A.out
206: filter: grep -v type | grep -v "Mat Object"
208: test:
209: requires: kokkos_kernels
210: suffix: 12_A_kokkos
211: args: -mat_type seqaijkokkos -rectA
212: output_file: output/ex2_12_A.out
213: filter: grep -v type | grep -v "Mat Object"
215: test:
216: suffix: 11_B
217: args: -mat_type seqaij -rectB
218: filter: grep -v "Mat Object"
220: test:
221: suffix: 12_B
222: args: -mat_type seqdense -rectB
223: filter: grep -v type | grep -v "Mat Object"
225: test:
226: requires: cuda
227: suffix: 12_B_cuda
228: args: -mat_type seqdensecuda -rectB
229: output_file: output/ex2_12_B.out
230: filter: grep -v type | grep -v "Mat Object"
232: test:
233: requires: kokkos_kernels
234: suffix: 12_B_kokkos
235: args: -mat_type seqaijkokkos -rectB
236: output_file: output/ex2_12_B.out
237: filter: grep -v type | grep -v "Mat Object"
239: test:
240: suffix: 21
241: args: -mat_type mpiaij
242: filter: grep -v type | grep -v "MPI processes"
244: test:
245: suffix: 22
246: args: -mat_type mpidense
247: filter: grep -v type | grep -v "Mat Object"
249: test:
250: requires: cuda
251: suffix: 22_cuda
252: output_file: output/ex2_22.out
253: args: -mat_type mpidensecuda
254: filter: grep -v type | grep -v "Mat Object"
256: test:
257: requires: kokkos_kernels
258: suffix: 22_kokkos
259: output_file: output/ex2_22.out
260: args: -mat_type mpiaijkokkos
261: filter: grep -v type | grep -v "Mat Object"
263: test:
264: suffix: 23
265: nsize: 3
266: args: -mat_type mpiaij
267: filter: grep -v type | grep -v "MPI processes"
269: test:
270: suffix: 24
271: nsize: 3
272: args: -mat_type mpidense
273: filter: grep -v type | grep -v "Mat Object"
275: test:
276: requires: cuda
277: suffix: 24_cuda
278: nsize: 3
279: output_file: output/ex2_24.out
280: args: -mat_type mpidensecuda
281: filter: grep -v type | grep -v "Mat Object"
283: test:
284: suffix: 2_aijcusparse_1
285: args: -mat_type mpiaijcusparse
286: output_file: output/ex2_21.out
287: requires: cuda
288: filter: grep -v type | grep -v "MPI processes"
290: test:
291: suffix: 2_aijkokkos_1
292: args: -mat_type mpiaijkokkos
293: output_file: output/ex2_21.out
294: requires: kokkos_kernels
295: filter: grep -v type | grep -v "MPI processes"
297: test:
298: suffix: 2_aijcusparse_2
299: nsize: 3
300: args: -mat_type mpiaijcusparse
301: output_file: output/ex2_23.out
302: requires: cuda
303: filter: grep -v type | grep -v "MPI processes"
305: test:
306: suffix: 2_aijkokkos_2
307: nsize: 3
308: args: -mat_type mpiaijkokkos
309: output_file: output/ex2_23.out
310: requires: kokkos_kernels
311: filter: grep -v type | grep -v "MPI processes"
313: test:
314: suffix: 3
315: nsize: 2
316: args: -mat_type mpiaij -rectA
318: test:
319: suffix: 3_aijcusparse
320: nsize: 2
321: args: -mat_type mpiaijcusparse -rectA
322: requires: cuda
324: test:
325: suffix: 4
326: nsize: 2
327: args: -mat_type mpidense -rectA
328: filter: grep -v type | grep -v "MPI processes"
330: test:
331: requires: cuda
332: suffix: 4_cuda
333: nsize: 2
334: output_file: output/ex2_4.out
335: args: -mat_type mpidensecuda -rectA
336: filter: grep -v type | grep -v "MPI processes"
338: test:
339: suffix: aijcusparse_1
340: args: -mat_type seqaijcusparse -rectA
341: filter: grep -v "Mat Object"
342: output_file: output/ex2_11_A_aijcusparse.out
343: requires: cuda
345: test:
346: suffix: aijcusparse_2
347: args: -mat_type seqaijcusparse -rectB
348: filter: grep -v "Mat Object"
349: output_file: output/ex2_11_B_aijcusparse.out
350: requires: cuda
352: TEST*/