Actual source code: ex109.c
petsc-3.7.3 2016-08-01
1: static char help[] = "Test MatMatMult() for AIJ and Dense matrices.\n\n";
3: #include <petscmat.h>
7: int main(int argc,char **argv)
8: {
9: Mat A,B,C,D;
10: PetscInt i,M,N,Istart,Iend,n=7,j,J,Ii,m=8,am,an;
11: PetscScalar v;
13: PetscRandom r;
14: PetscBool equal=PETSC_FALSE;
15: PetscReal fill = 1.0;
16: PetscMPIInt size;
18: PetscInitialize(&argc,&argv,(char*)0,help);
19: PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
20: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
21: PetscOptionsGetReal(NULL,NULL,"-fill",&fill,NULL);
23: PetscRandomCreate(PETSC_COMM_WORLD,&r);
24: PetscRandomSetFromOptions(r);
26: /* Create a aij matrix A */
27: M = N = m*n;
28: MatCreate(PETSC_COMM_WORLD,&A);
29: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,N);
30: MatSetType(A,MATAIJ);
31: MatSetFromOptions(A);
32: MatMPIAIJSetPreallocation(A,5,NULL,5,NULL);
33: MatSeqAIJSetPreallocation(A,5,NULL);
35: MatGetOwnershipRange(A,&Istart,&Iend);
36: am = Iend - Istart;
37: for (Ii=Istart; Ii<Iend; Ii++) {
38: v = -1.0; i = Ii/n; j = Ii - i*n;
39: if (i>0) {J = Ii - n; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
40: if (i<m-1) {J = Ii + n; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
41: if (j>0) {J = Ii - 1; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
42: if (j<n-1) {J = Ii + 1; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
43: v = 4.0; MatSetValues(A,1,&Ii,1,&Ii,&v,INSERT_VALUES);
44: }
45: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
46: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
48: /* Create a dense matrix B */
49: MatGetLocalSize(A,&am,&an);
50: MatCreate(PETSC_COMM_WORLD,&B);
51: MatSetSizes(B,an,PETSC_DECIDE,PETSC_DECIDE,M);
52: MatSetType(B,MATDENSE);
53: MatSeqDenseSetPreallocation(B,NULL);
54: MatMPIDenseSetPreallocation(B,NULL);
55: MatSetFromOptions(B);
56: MatSetRandom(B,r);
57: PetscRandomDestroy(&r);
58: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
59: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
61: /* Test C = A*B (aij*dense) */
62: MatMatMult(A,B,MAT_INITIAL_MATRIX,fill,&C);
63: MatMatMult(A,B,MAT_REUSE_MATRIX,fill,&C);
65: MatMatMultSymbolic(A,B,fill,&D);
66: for (i=0; i<2; i++) {
67: MatMatMultNumeric(A,B,D);
68: }
69: MatEqual(C,D,&equal);
70: if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"C != D");
71: MatDestroy(&D);
73: /* Test D = C*A (dense*aij) */
74: MatMatMult(C,A,MAT_INITIAL_MATRIX,fill,&D);
75: MatMatMult(C,A,MAT_REUSE_MATRIX,fill,&D);
76: MatDestroy(&D);
78: /* Test D = A*C (aij*dense) */
79: MatMatMult(A,C,MAT_INITIAL_MATRIX,fill,&D);
80: MatMatMult(A,C,MAT_REUSE_MATRIX,fill,&D);
81: MatDestroy(&D);
83: /* Test D = B*C (dense*dense) */
84: MPI_Comm_size(PETSC_COMM_WORLD,&size);
85: if (size == 1) {
86: MatMatMult(B,C,MAT_INITIAL_MATRIX,fill,&D);
87: MatMatMult(B,C,MAT_REUSE_MATRIX,fill,&D);
88: MatDestroy(&D);
89: }
91: MatDestroy(&C);
92: MatDestroy(&B);
93: MatDestroy(&A);
94: PetscFinalize();
95: return(0);
96: }