Actual source code: ex93.c
petsc-3.13.6 2020-09-29
1: static char help[] = "Test MatMatMult() and MatPtAP() for AIJ matrices.\n\n";
3: #include <petscmat.h>
5: extern PetscErrorCode testPTAPRectangular(void);
7: int main(int argc,char **argv)
8: {
9: Mat A,B,C,D;
10: PetscScalar a[] ={1.,1.,0.,0.,1.,1.,0.,0.,1.};
11: PetscInt ij[]={0,1,2};
13: PetscReal fill=4.0;
14: PetscMPIInt size,rank;
15: PetscBool isequal;
16: #if defined(PETSC_HAVE_HYPRE)
17: PetscBool test_hypre=PETSC_FALSE;
18: #endif
20: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
21: #if defined(PETSC_HAVE_HYPRE)
22: PetscOptionsGetBool(NULL,NULL,"-test_hypre",&test_hypre,NULL);
23: #endif
24: MPI_Comm_size(PETSC_COMM_WORLD,&size);
25: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
27: MatCreate(PETSC_COMM_WORLD,&A);
28: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,3,3);
29: MatSetType(A,MATAIJ);
30: MatSetFromOptions(A);
31: MatSetUp(A);
32: MatSetOption(A,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);
34: if (!rank) {
35: MatSetValues(A,3,ij,3,ij,a,ADD_VALUES);
36: }
37: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
38: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
40: /* Test MatMatMult() */
41: MatTranspose(A,MAT_INITIAL_MATRIX,&B); /* B = A^T */
42: MatMatMult(B,A,MAT_INITIAL_MATRIX,fill,&C); /* C = B*A */
43: MatMatMult(B,A,MAT_REUSE_MATRIX,fill,&C); /* recompute C=B*A */
44: MatSetOptionsPrefix(C,"C_");
45: MatMatMultEqual(B,A,C,10,&isequal);
46: if (!isequal) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_INCOMP,"MatMatMult: C != B*A");
48: MatMatMult(C,A,MAT_INITIAL_MATRIX,fill,&D); /* D = C*A = (A^T*A)*A */
49: MatMatMult(C,A,MAT_REUSE_MATRIX,fill,&D);
50: MatMatMultEqual(C,A,D,10,&isequal);
51: if (!isequal) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_INCOMP,"MatMatMult: D != C*A");
53: MatDestroy(&B);
54: MatDestroy(&C);
55: MatDestroy(&D);
57: /* Test MatPtAP */
58: MatDuplicate(A,MAT_COPY_VALUES,&B); /* B = A */
59: MatPtAP(A,B,MAT_INITIAL_MATRIX,fill,&C); /* C = B^T*A*B */
60: MatPtAPMultEqual(A,B,C,10,&isequal);
61: if (!isequal) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_INCOMP,"MatPtAP: C != B^T*A*B");
63: /* Repeat MatPtAP to test symbolic/numeric separation for reuse of the symbolic product */
64: MatPtAP(A,B,MAT_REUSE_MATRIX,fill,&C);
65: MatPtAPMultEqual(A,B,C,10,&isequal);
66: if (!isequal) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_INCOMP,"MatPtAP(reuse): C != B^T*A*B");
68: MatDestroy(&C);
69: MatDestroy(&B);
71: if (size == 1) {
72: /* A test contributed by Tobias Neckel <neckel@in.tum.de> */
73: testPTAPRectangular();
75: /* test MatMatTransposeMult(): A*B^T */
76: MatMatTransposeMult(A,A,MAT_INITIAL_MATRIX,fill,&D); /* D = A*A^T */
77: MatScale(A,2.0);
78: MatMatTransposeMult(A,A,MAT_REUSE_MATRIX,fill,&D);
79: MatMatTransposeMultEqual(A,A,D,10,&isequal);
80: if (!isequal) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_INCOMP,"MatMatTranspose: D != A*A^T");
81: }
83: MatDestroy(&A);
84: MatDestroy(&B);
85: MatDestroy(&C);
86: MatDestroy(&D);
87: PetscFinalize();
88: return ierr;
89: }
91: /* a test contributed by Tobias Neckel <neckel@in.tum.de>, 02 Jul 2008 */
92: PetscErrorCode testPTAPRectangular(void)
93: {
94: const int rows = 3,cols = 5;
96: int i;
97: Mat A,P,C;
100: /* set up A */
101: MatCreateSeqAIJ(PETSC_COMM_WORLD, rows, rows,1, NULL, &A);
102: for (i=0; i<rows; i++) {
103: MatSetValue(A, i, i, 1.0, INSERT_VALUES);
104: }
105: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
106: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
108: /* set up P */
109: MatCreateSeqAIJ(PETSC_COMM_WORLD, rows, cols,5, NULL, &P);
110: MatSetValue(P, 0, 0, 1.0, INSERT_VALUES);
111: MatSetValue(P, 0, 1, 2.0, INSERT_VALUES);
112: MatSetValue(P, 0, 2, 0.0, INSERT_VALUES);
114: MatSetValue(P, 0, 3, -1.0, INSERT_VALUES);
116: MatSetValue(P, 1, 0, 0.0, INSERT_VALUES);
117: MatSetValue(P, 1, 1, -1.0, INSERT_VALUES);
118: MatSetValue(P, 1, 2, 1.0, INSERT_VALUES);
120: MatSetValue(P, 2, 0, 3.0, INSERT_VALUES);
121: MatSetValue(P, 2, 1, 0.0, INSERT_VALUES);
122: MatSetValue(P, 2, 2, -3.0, INSERT_VALUES);
124: MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);
125: MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);
127: /* compute C */
128: MatPtAP(A, P, MAT_INITIAL_MATRIX, 1.0, &C);
130: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
131: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
133: /* compare results */
134: /*
135: printf("C:\n");
136: MatView(C,PETSC_VIEWER_STDOUT_WORLD);
138: blitz::Array<double,2> actualC(cols, cols);
139: actualC = 0.0;
140: for (int i=0; i<cols; i++) {
141: for (int j=0; j<cols; j++) {
142: MatGetValues(C, 1, &i, 1, &j, &actualC(i,j));
143: ;
144: }
145: }
146: blitz::Array<double,2> expectedC(cols, cols);
147: expectedC = 0.0;
149: expectedC(0,0) = 10.0;
150: expectedC(0,1) = 2.0;
151: expectedC(0,2) = -9.0;
152: expectedC(0,3) = -1.0;
153: expectedC(1,0) = 2.0;
154: expectedC(1,1) = 5.0;
155: expectedC(1,2) = -1.0;
156: expectedC(1,3) = -2.0;
157: expectedC(2,0) = -9.0;
158: expectedC(2,1) = -1.0;
159: expectedC(2,2) = 10.0;
160: expectedC(2,3) = 0.0;
161: expectedC(3,0) = -1.0;
162: expectedC(3,1) = -2.0;
163: expectedC(3,2) = 0.0;
164: expectedC(3,3) = 1.0;
166: int check = areBlitzArrays2NumericallyEqual(actualC,expectedC);
167: validateEqualsWithParams3(check, -1 , "testPTAPRectangular()", check, actualC(check), expectedC(check));
168: */
170: MatDestroy(&A);
171: MatDestroy(&P);
172: MatDestroy(&C);
173: return(0);
174: }
176: /*TEST
178: test:
180: test:
181: suffix: 2
182: nsize: 2
183: args: -matmatmult_via nonscalable
184: output_file: output/ex93_1.out
186: test:
187: suffix: 3
188: nsize: 2
189: output_file: output/ex93_1.out
191: test:
192: suffix: 4
193: nsize: 2
194: args: -matptap_via scalable
195: output_file: output/ex93_1.out
197: test:
198: suffix: btheap
199: args: -matmatmult_via btheap -matmattransmult_via color
200: output_file: output/ex93_1.out
202: test:
203: suffix: heap
204: args: -matmatmult_via heap
205: output_file: output/ex93_1.out
207: #HYPRE PtAP is broken for complex numbers
208: test:
209: suffix: hypre
210: nsize: 3
211: requires: hypre !complex
212: args: -matmatmult_via hypre -matptap_via hypre -test_hypre
213: output_file: output/ex93_hypre.out
215: test:
216: suffix: llcondensed
217: args: -matmatmult_via llcondensed
218: output_file: output/ex93_1.out
220: test:
221: suffix: scalable
222: args: -matmatmult_via scalable
223: output_file: output/ex93_1.out
225: test:
226: suffix: scalable_fast
227: args: -matmatmult_via scalable_fast
228: output_file: output/ex93_1.out
230: TEST*/