Actual source code: ex2.c
petsc-3.12.5 2020-03-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, test MatValid() --------- */
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: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
176: /* Free data structures */
177: MatDestroy(&mat);
178: MatDestroy(&tmat);
179: PetscFinalize();
180: return ierr;
181: }
187: /*TEST
189: test:
190: suffix: 11_A
191: args: -mat_type seqaij -rectA
192: filter: grep -v "Mat Object"
194: test:
195: suffix: 12_A
196: args: -mat_type seqdense -rectA
197: filter: grep -v "Mat Object"
199: test:
200: suffix: 11_B
201: args: -mat_type seqaij -rectB
202: filter: grep -v "Mat Object"
204: test:
205: suffix: 12_B
206: args: -mat_type seqdense -rectB
207: filter: grep -v "Mat Object"
209: test:
210: suffix: 21
211: args: -mat_type mpiaij
213: test:
214: suffix: 22
215: args: -mat_type mpidense
217: test:
218: suffix: 23
219: nsize: 3
220: args: -mat_type mpiaij
221: filter: grep -v type | grep -v "MPI processes"
223: test:
224: suffix: 24
225: nsize: 3
226: args: -mat_type mpidense
228: test:
229: suffix: 2_aijcusparse_1
230: args: -mat_type mpiaijcusparse
231: output_file: output/ex2_23.out
232: requires: cuda
233: filter: grep -v type | grep -v "MPI processes"
235: test:
236: suffix: 2_aijcusparse_2
237: nsize: 3
238: args: -mat_type mpiaijcusparse
239: output_file: output/ex2_23.out
240: requires: cuda
241: filter: grep -v type | grep -v "MPI processes"
243: test:
244: suffix: 3
245: nsize: 2
246: args: -mat_type mpiaij -rectA
248: test:
249: suffix: 3_aijcusparse
250: nsize: 2
251: args: -mat_type mpiaijcusparse -rectA
252: requires: cuda
254: test:
255: suffix: 4
256: nsize: 2
257: args: -mat_type mpidense -rectA
259: test:
260: suffix: aijcusparse_1
261: args: -mat_type seqaijcusparse -rectA
262: filter: grep -v "Mat Object"
263: output_file: output/ex2_11_A_aijcusparse.out
264: requires: cuda
266: test:
267: suffix: aijcusparse_2
268: args: -mat_type seqaijcusparse -rectB
269: filter: grep -v "Mat Object"
270: output_file: output/ex2_11_B_aijcusparse.out
271: requires: cuda
273: TEST*/