Actual source code: ex2.c
petsc-3.14.0 2020-09-29
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: suffix: 11_B
210: args: -mat_type seqaij -rectB
211: filter: grep -v "Mat Object"
213: test:
214: suffix: 12_B
215: args: -mat_type seqdense -rectB
216: filter: grep -v type | grep -v "Mat Object"
218: test:
219: requires: cuda
220: suffix: 12_B_cuda
221: args: -mat_type seqdensecuda -rectB
222: output_file: output/ex2_12_B.out
223: filter: grep -v type | grep -v "Mat Object"
225: test:
226: suffix: 21
227: args: -mat_type mpiaij
228: filter: grep -v type | grep -v "MPI processes"
230: test:
231: suffix: 22
232: args: -mat_type mpidense
233: filter: grep -v type | grep -v "Mat Object"
235: test:
236: requires: cuda
237: suffix: 22_cuda
238: output_file: output/ex2_22.out
239: args: -mat_type mpidensecuda
240: filter: grep -v type | grep -v "Mat Object"
242: test:
243: suffix: 23
244: nsize: 3
245: args: -mat_type mpiaij
246: filter: grep -v type | grep -v "MPI processes"
248: test:
249: suffix: 24
250: nsize: 3
251: args: -mat_type mpidense
252: filter: grep -v type | grep -v "Mat Object"
254: test:
255: requires: cuda
256: suffix: 24_cuda
257: nsize: 3
258: output_file: output/ex2_24.out
259: args: -mat_type mpidensecuda
260: filter: grep -v type | grep -v "Mat Object"
262: test:
263: suffix: 2_aijcusparse_1
264: args: -mat_type mpiaijcusparse
265: output_file: output/ex2_21.out
266: requires: cuda
267: filter: grep -v type | grep -v "MPI processes"
269: test:
270: suffix: 2_aijcusparse_2
271: nsize: 3
272: args: -mat_type mpiaijcusparse
273: output_file: output/ex2_23.out
274: requires: cuda
275: filter: grep -v type | grep -v "MPI processes"
277: test:
278: suffix: 3
279: nsize: 2
280: args: -mat_type mpiaij -rectA
282: test:
283: suffix: 3_aijcusparse
284: nsize: 2
285: args: -mat_type mpiaijcusparse -rectA
286: requires: cuda
288: test:
289: suffix: 4
290: nsize: 2
291: args: -mat_type mpidense -rectA
292: filter: grep -v type | grep -v "MPI processes"
294: test:
295: requires: cuda
296: suffix: 4_cuda
297: nsize: 2
298: output_file: output/ex2_4.out
299: args: -mat_type mpidensecuda -rectA
300: filter: grep -v type | grep -v "MPI processes"
302: test:
303: suffix: aijcusparse_1
304: args: -mat_type seqaijcusparse -rectA
305: filter: grep -v "Mat Object"
306: output_file: output/ex2_11_A_aijcusparse.out
307: requires: cuda
309: test:
310: suffix: aijcusparse_2
311: args: -mat_type seqaijcusparse -rectB
312: filter: grep -v "Mat Object"
313: output_file: output/ex2_11_B_aijcusparse.out
314: requires: cuda
316: TEST*/