Actual source code: ex89.c
1: static char help[] ="Tests MatPtAP() for MPIMAIJ and MPIAIJ \n ";
3: #include <petscdmda.h>
5: int main(int argc,char **argv)
6: {
7: DM coarsedm,finedm;
8: PetscMPIInt size,rank;
9: PetscInt M,N,Z,i,nrows;
10: PetscScalar one = 1.0;
11: PetscReal fill=2.0;
12: Mat A,P,C;
13: PetscScalar *array,alpha;
14: PetscBool Test_3D=PETSC_FALSE,flg;
15: const PetscInt *ia,*ja;
16: PetscInt dof;
17: MPI_Comm comm;
19: PetscInitialize(&argc,&argv,NULL,help);
20: comm = PETSC_COMM_WORLD;
21: MPI_Comm_rank(comm,&rank);
22: MPI_Comm_size(comm,&size);
23: M = 10; N = 10; Z = 10;
24: dof = 10;
26: PetscOptionsGetBool(NULL,NULL,"-test_3D",&Test_3D,NULL);
27: PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);
28: PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);
29: PetscOptionsGetInt(NULL,NULL,"-Z",&Z,NULL);
30: /* Set up distributed array for fine grid */
31: if (!Test_3D) {
32: DMDACreate2d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,M,N,PETSC_DECIDE,PETSC_DECIDE,dof,1,NULL,NULL,&coarsedm);
33: } else {
34: DMDACreate3d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,M,N,Z,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,dof,1,NULL,NULL,NULL,&coarsedm);
35: }
36: DMSetFromOptions(coarsedm);
37: DMSetUp(coarsedm);
39: /* This makes sure the coarse DMDA has the same partition as the fine DMDA */
40: DMRefine(coarsedm,PetscObjectComm((PetscObject)coarsedm),&finedm);
42: /*------------------------------------------------------------*/
43: DMSetMatType(finedm,MATAIJ);
44: DMCreateMatrix(finedm,&A);
46: /* set val=one to A */
47: if (size == 1) {
48: MatGetRowIJ(A,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);
49: if (flg) {
50: MatSeqAIJGetArray(A,&array);
51: for (i=0; i<ia[nrows]; i++) array[i] = one;
52: MatSeqAIJRestoreArray(A,&array);
53: }
54: MatRestoreRowIJ(A,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);
55: } else {
56: Mat AA,AB;
57: MatMPIAIJGetSeqAIJ(A,&AA,&AB,NULL);
58: MatGetRowIJ(AA,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);
59: if (flg) {
60: MatSeqAIJGetArray(AA,&array);
61: for (i=0; i<ia[nrows]; i++) array[i] = one;
62: MatSeqAIJRestoreArray(AA,&array);
63: }
64: MatRestoreRowIJ(AA,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);
65: MatGetRowIJ(AB,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);
66: if (flg) {
67: MatSeqAIJGetArray(AB,&array);
68: for (i=0; i<ia[nrows]; i++) array[i] = one;
69: MatSeqAIJRestoreArray(AB,&array);
70: }
71: MatRestoreRowIJ(AB,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);
72: }
73: /* Create interpolation between the fine and coarse grids */
74: DMCreateInterpolation(coarsedm,finedm,&P,NULL);
76: /* Test P^T * A * P - MatPtAP() */
77: /*------------------------------*/
78: /* (1) Developer API */
79: MatProductCreate(A,P,NULL,&C);
80: MatProductSetType(C,MATPRODUCT_PtAP);
81: MatProductSetAlgorithm(C,"allatonce");
82: MatProductSetFill(C,PETSC_DEFAULT);
83: MatProductSetFromOptions(C);
84: MatProductSymbolic(C);
85: MatProductNumeric(C);
86: MatProductNumeric(C); /* Test reuse of symbolic C */
88: { /* Test MatProductView() */
89: PetscViewer viewer;
90: PetscViewerASCIIOpen(comm,NULL, &viewer);
91: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
92: MatProductView(C,viewer);
93: PetscViewerPopFormat(viewer);
94: PetscViewerDestroy(&viewer);
95: }
97: MatPtAPMultEqual(A,P,C,10,&flg);
99: MatDestroy(&C);
101: /* (2) User API */
102: MatPtAP(A,P,MAT_INITIAL_MATRIX,fill,&C);
103: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
104: alpha=1.0;
105: for (i=0; i<1; i++) {
106: alpha -= 0.1;
107: MatScale(A,alpha);
108: MatPtAP(A,P,MAT_REUSE_MATRIX,fill,&C);
109: }
111: /* Free intermediate data structures created for reuse of C=Pt*A*P */
112: MatProductClear(C);
114: MatPtAPMultEqual(A,P,C,10,&flg);
117: MatDestroy(&C);
118: MatDestroy(&A);
119: MatDestroy(&P);
120: DMDestroy(&finedm);
121: DMDestroy(&coarsedm);
122: PetscFinalize();
123: return 0;
124: }
126: /*TEST
128: test:
129: args: -M 10 -N 10 -Z 10
130: output_file: output/ex89_1.out
132: test:
133: suffix: allatonce
134: nsize: 4
135: args: -M 10 -N 10 -Z 10
136: output_file: output/ex89_2.out
138: test:
139: suffix: allatonce_merged
140: nsize: 4
141: args: -M 10 -M 5 -M 10 -mat_product_algorithm allatonce_merged
142: output_file: output/ex89_3.out
144: test:
145: suffix: nonscalable_3D
146: nsize: 4
147: args: -M 10 -M 5 -M 10 -test_3D 1 -mat_product_algorithm nonscalable
148: output_file: output/ex89_4.out
150: test:
151: suffix: allatonce_merged_3D
152: nsize: 4
153: args: -M 10 -M 5 -M 10 -test_3D 1 -mat_product_algorithm allatonce_merged
154: output_file: output/ex89_3.out
156: test:
157: suffix: nonscalable
158: nsize: 4
159: args: -M 10 -N 10 -Z 10 -mat_product_algorithm nonscalable
160: output_file: output/ex89_5.out
162: TEST*/