Actual source code: ex109.c
petsc-3.14.6 2021-03-30
1: static char help[] = "Test MatMatMult() for AIJ and Dense matrices.\n\n";
3: #include <petscmat.h>
5: int main(int argc,char **argv)
6: {
7: Mat A,B,C,D,AT;
8: PetscInt i,M,N,Istart,Iend,n=7,j,J,Ii,m=8,am,an;
9: PetscScalar v;
11: PetscRandom r;
12: PetscBool equal=PETSC_FALSE,flg;
13: PetscReal fill = 1.0,norm;
14: PetscMPIInt size;
16: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
17: PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
18: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
19: PetscOptionsGetReal(NULL,NULL,"-fill",&fill,NULL);
21: PetscRandomCreate(PETSC_COMM_WORLD,&r);
22: PetscRandomSetFromOptions(r);
24: /* Create a aij matrix A */
25: M = N = m*n;
26: MatCreate(PETSC_COMM_WORLD,&A);
27: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,N);
28: MatSetType(A,MATAIJ);
29: MatSetFromOptions(A);
30: MatMPIAIJSetPreallocation(A,5,NULL,5,NULL);
31: MatSeqAIJSetPreallocation(A,5,NULL);
33: MatGetOwnershipRange(A,&Istart,&Iend);
34: am = Iend - Istart;
35: for (Ii=Istart; Ii<Iend; Ii++) {
36: v = -1.0; i = Ii/n; j = Ii - i*n;
37: if (i>0) {J = Ii - n; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
38: if (i<m-1) {J = Ii + n; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
39: if (j>0) {J = Ii - 1; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
40: if (j<n-1) {J = Ii + 1; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
41: v = 4.0; MatSetValues(A,1,&Ii,1,&Ii,&v,INSERT_VALUES);
42: }
43: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
44: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
46: /* Create a dense matrix B */
47: MatGetLocalSize(A,&am,&an);
48: MatCreate(PETSC_COMM_WORLD,&B);
49: MatSetSizes(B,an,PETSC_DECIDE,PETSC_DECIDE,M);
50: MatSetType(B,MATDENSE);
51: MatSeqDenseSetPreallocation(B,NULL);
52: MatMPIDenseSetPreallocation(B,NULL);
53: MatSetFromOptions(B);
54: MatSetRandom(B,r);
55: PetscRandomDestroy(&r);
56: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
57: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
59: /* Test reuse of user-provided dense C (unassembled) -- not recommended usage */
60: MatCreate(PETSC_COMM_WORLD,&C);
61: MatSetType(C,MATDENSE);
62: MatSetSizes(C,am,PETSC_DECIDE,PETSC_DECIDE,N);
63: MatSetFromOptions(C);
64: MatSetUp(C);
65: MatZeroEntries(C);
66: MatMatMult(A,B,MAT_REUSE_MATRIX,fill,&C);
67: MatNorm(C,NORM_INFINITY,&norm);
68: MatDestroy(&C);
70: /* Test C = A*B (aij*dense) */
71: MatMatMult(A,B,MAT_INITIAL_MATRIX,fill,&C);
72: MatMatMult(A,B,MAT_REUSE_MATRIX,fill,&C);
74: /* Test developer API */
75: MatProductCreate(A,B,NULL,&D);
76: MatProductSetType(D,MATPRODUCT_AB);
77: MatProductSetAlgorithm(D,"default");
78: MatProductSetFill(D,fill);
79: MatProductSetFromOptions(D);
80: MatProductSymbolic(D);
81: for (i=0; i<2; i++) {
82: MatProductNumeric(D);
83: }
84: MatEqual(C,D,&equal);
85: if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"C != D");
86: MatDestroy(&D);
88: /* Test D = AT*B (transpose(aij)*dense) */
89: MatCreateTranspose(A,&AT);
90: MatMatMult(AT,B,MAT_INITIAL_MATRIX,fill,&D);
91: MatMatMultEqual(AT,B,D,10,&equal);
92: if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"D != AT*B (transpose(aij)*dense)");
93: MatDestroy(&D);
94: MatDestroy(&AT);
96: /* Test D = C*A (dense*aij) */
97: MatMatMult(C,A,MAT_INITIAL_MATRIX,fill,&D);
98: MatMatMult(C,A,MAT_REUSE_MATRIX,fill,&D);
99: MatMatMultEqual(C,A,D,10,&equal);
100: if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"D != C*A (dense*aij)");
101: MatDestroy(&D);
103: /* Test D = A*C (aij*dense) */
104: MatMatMult(A,C,MAT_INITIAL_MATRIX,fill,&D);
105: MatMatMult(A,C,MAT_REUSE_MATRIX,fill,&D);
106: MatMatMultEqual(A,C,D,10,&equal);
107: if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"D != A*C (aij*dense)");
108: MatDestroy(&D);
110: /* Test D = B*C (dense*dense) */
111: MPI_Comm_size(PETSC_COMM_WORLD,&size);
112: if (size == 1) {
113: MatMatMult(B,C,MAT_INITIAL_MATRIX,fill,&D);
114: MatMatMult(B,C,MAT_REUSE_MATRIX,fill,&D);
115: MatMatMultEqual(B,C,D,10,&equal);
116: if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"D != B*C (dense*dense)");
117: MatDestroy(&D);
118: }
120: /* Test D = B*C^T (dense*dense) */
121: MatMatTransposeMult(B,C,MAT_INITIAL_MATRIX,fill,&D);
122: MatMatTransposeMult(B,C,MAT_REUSE_MATRIX,fill,&D);
123: MatMatTransposeMultEqual(B,C,D,10,&equal);
124: if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"D != B*C^T (dense*dense)");
125: MatDestroy(&D);
127: /* Test MatProductCreateWithMat() and reuse C and B for B = A*C */
128: flg = PETSC_FALSE;
129: PetscOptionsHasName(NULL,NULL,"-test_userAPI",&flg);
130: if (flg) {
131: /* user driver */
132: MatMatMult(A,C,MAT_REUSE_MATRIX,fill,&B);
133: } else {
134: /* clear internal data structures related with previous products to avoid circular references */
135: MatProductClear(A);
136: MatProductClear(B);
137: MatProductClear(C);
138: MatProductCreateWithMat(A,C,NULL,B);
139: MatProductSetType(B,MATPRODUCT_AB);
140: MatProductSetFromOptions(B);
141: MatProductSymbolic(B);
142: MatProductNumeric(B);
143: }
145: MatDestroy(&C);
146: MatDestroy(&B);
147: MatDestroy(&A);
148: PetscFinalize();
149: return ierr;
150: }
153: /*TEST
155: test:
156: args: -M 10 -N 10
157: output_file: output/ex109.out
159: test:
160: suffix: 2
161: nsize: 3
162: output_file: output/ex109.out
164: test:
165: suffix: 3
166: nsize: 2
167: args: -matmattransmult_mpidense_mpidense_via cyclic
168: output_file: output/ex109.out
170: test:
171: suffix: 4
172: args: -test_userAPI
173: output_file: output/ex109.out
175: test:
176: suffix: 5
177: nsize: 3
178: args: -test_userAPI
179: output_file: output/ex109.out
181: TEST*/