Actual source code: ex109.c
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;
10: PetscRandom r;
11: PetscBool equal=PETSC_FALSE,flg;
12: PetscReal fill = 1.0,norm;
13: PetscMPIInt size;
15: PetscInitialize(&argc,&argv,(char*)0,help);
16: PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
17: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
18: PetscOptionsGetReal(NULL,NULL,"-fill",&fill,NULL);
20: PetscRandomCreate(PETSC_COMM_WORLD,&r);
21: PetscRandomSetFromOptions(r);
23: /* Create a aij matrix A */
24: M = N = m*n;
25: MatCreate(PETSC_COMM_WORLD,&A);
26: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,N);
27: MatSetType(A,MATAIJ);
28: MatSetFromOptions(A);
29: MatMPIAIJSetPreallocation(A,5,NULL,5,NULL);
30: MatSeqAIJSetPreallocation(A,5,NULL);
32: MatGetOwnershipRange(A,&Istart,&Iend);
33: am = Iend - Istart;
34: for (Ii=Istart; Ii<Iend; Ii++) {
35: v = -1.0; i = Ii/n; j = Ii - i*n;
36: if (i>0) {J = Ii - n; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
37: if (i<m-1) {J = Ii + n; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
38: if (j>0) {J = Ii - 1; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
39: if (j<n-1) {J = Ii + 1; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
40: v = 4.0; MatSetValues(A,1,&Ii,1,&Ii,&v,INSERT_VALUES);
41: }
42: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
43: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
45: /* Create a dense matrix B */
46: MatGetLocalSize(A,&am,&an);
47: MatCreate(PETSC_COMM_WORLD,&B);
48: MatSetSizes(B,an,PETSC_DECIDE,PETSC_DECIDE,M);
49: MatSetType(B,MATDENSE);
50: MatSeqDenseSetPreallocation(B,NULL);
51: MatMPIDenseSetPreallocation(B,NULL);
52: MatSetFromOptions(B);
53: MatSetRandom(B,r);
54: PetscRandomDestroy(&r);
55: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
56: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
58: /* Test reuse of user-provided dense C (unassembled) -- not recommended usage */
59: MatCreate(PETSC_COMM_WORLD,&C);
60: MatSetType(C,MATDENSE);
61: MatSetSizes(C,am,PETSC_DECIDE,PETSC_DECIDE,N);
62: MatSetFromOptions(C);
63: MatSetUp(C);
64: MatZeroEntries(C);
65: MatMatMult(A,B,MAT_REUSE_MATRIX,fill,&C);
66: MatNorm(C,NORM_INFINITY,&norm);
67: MatDestroy(&C);
69: /* Test C = A*B (aij*dense) */
70: MatMatMult(A,B,MAT_INITIAL_MATRIX,fill,&C);
71: MatMatMult(A,B,MAT_REUSE_MATRIX,fill,&C);
73: /* Test developer API */
74: MatProductCreate(A,B,NULL,&D);
75: MatProductSetType(D,MATPRODUCT_AB);
76: MatProductSetAlgorithm(D,"default");
77: MatProductSetFill(D,fill);
78: MatProductSetFromOptions(D);
79: MatProductSymbolic(D);
80: for (i=0; i<2; i++) {
81: MatProductNumeric(D);
82: }
83: MatEqual(C,D,&equal);
85: MatDestroy(&D);
87: /* Test D = AT*B (transpose(aij)*dense) */
88: MatCreateTranspose(A,&AT);
89: MatMatMult(AT,B,MAT_INITIAL_MATRIX,fill,&D);
90: MatMatMultEqual(AT,B,D,10,&equal);
92: MatDestroy(&D);
93: MatDestroy(&AT);
95: /* Test D = C*A (dense*aij) */
96: MatMatMult(C,A,MAT_INITIAL_MATRIX,fill,&D);
97: MatMatMult(C,A,MAT_REUSE_MATRIX,fill,&D);
98: MatMatMultEqual(C,A,D,10,&equal);
100: MatDestroy(&D);
102: /* Test D = A*C (aij*dense) */
103: MatMatMult(A,C,MAT_INITIAL_MATRIX,fill,&D);
104: MatMatMult(A,C,MAT_REUSE_MATRIX,fill,&D);
105: MatMatMultEqual(A,C,D,10,&equal);
107: MatDestroy(&D);
109: /* Test D = B*C (dense*dense) */
110: MPI_Comm_size(PETSC_COMM_WORLD,&size);
111: if (size == 1) {
112: MatMatMult(B,C,MAT_INITIAL_MATRIX,fill,&D);
113: MatMatMult(B,C,MAT_REUSE_MATRIX,fill,&D);
114: MatMatMultEqual(B,C,D,10,&equal);
116: MatDestroy(&D);
117: }
119: /* Test D = B*C^T (dense*dense) */
120: MatMatTransposeMult(B,C,MAT_INITIAL_MATRIX,fill,&D);
121: MatMatTransposeMult(B,C,MAT_REUSE_MATRIX,fill,&D);
122: MatMatTransposeMultEqual(B,C,D,10,&equal);
124: MatDestroy(&D);
126: /* Test MatProductCreateWithMat() and reuse C and B for B = A*C */
127: flg = PETSC_FALSE;
128: PetscOptionsHasName(NULL,NULL,"-test_userAPI",&flg);
129: if (flg) {
130: /* user driver */
131: MatMatMult(A,C,MAT_REUSE_MATRIX,fill,&B);
132: } else {
133: /* clear internal data structures related with previous products to avoid circular references */
134: MatProductClear(A);
135: MatProductClear(B);
136: MatProductClear(C);
137: MatProductCreateWithMat(A,C,NULL,B);
138: MatProductSetType(B,MATPRODUCT_AB);
139: MatProductSetFromOptions(B);
140: MatProductSymbolic(B);
141: MatProductNumeric(B);
142: }
144: MatDestroy(&C);
145: MatDestroy(&B);
146: MatDestroy(&A);
147: PetscFinalize();
148: return 0;
149: }
151: /*TEST
153: test:
154: args: -M 10 -N 10
155: output_file: output/ex109.out
157: test:
158: suffix: 2
159: nsize: 3
160: output_file: output/ex109.out
162: test:
163: suffix: 3
164: nsize: 2
165: args: -matmattransmult_mpidense_mpidense_via cyclic
166: output_file: output/ex109.out
168: test:
169: suffix: 4
170: args: -test_userAPI
171: output_file: output/ex109.out
173: test:
174: suffix: 5
175: nsize: 3
176: args: -test_userAPI
177: output_file: output/ex109.out
179: TEST*/