1
2 static char help[] ="Tests sequential and parallel MatMatMatMult() and MatPtAP(). Modified from ex96.c \n\
3 -Mx <xg>, where <xg> = number of coarse grid points in the x-direction\n\
4 -My <yg>, where <yg> = number of coarse grid points in the y-direction\n\
5 -Mz <zg>, where <zg> = number of coarse grid points in the z-direction\n\
6 -Npx <npx>, where <npx> = number of processors in the x-direction\n\
7 -Npy <npy>, where <npy> = number of processors in the y-direction\n\
8 -Npz <npz>, where <npz> = number of processors in the z-direction\n\n";
9
10 /*
11 Example of usage: mpiexec -n 3 ./ex41 -Mx 10 -My 10 -Mz 10
12 */
13
14 #include <petscdm.h>
15 #include <petscdmda.h>
16
17 /* User-defined application contexts */
18 typedef struct {
19 PetscInt mx,my,mz; /* number grid points in x, y and z direction */
20 Vec localX,localF; /* local vectors with ghost region */
21 DM da;
22 Vec x,b,r; /* global vectors */
23 Mat J; /* Jacobian on grid */
24 } GridCtx;
25 typedef struct {
26 GridCtx fine;
27 GridCtx coarse;
28 PetscInt ratio;
29 Mat Ii; /* interpolation from coarse to fine */
30 } AppCtx;
31
32 #define COARSE_LEVEL 0
33 #define FINE_LEVEL 1
34
35 /*
36 Mm_ratio - ration of grid lines between fine and coarse grids.
37 */
main(int argc,char ** argv)38 int main(int argc,char **argv)
39 {
40 PetscErrorCode ierr;
41 AppCtx user;
42 PetscMPIInt size,rank;
43 PetscInt m,n,M,N,i,nrows;
44 PetscScalar one = 1.0;
45 PetscReal fill=2.0;
46 Mat A,P,R,C,PtAP,D;
47 PetscScalar *array;
48 PetscRandom rdm;
49 PetscBool Test_3D=PETSC_FALSE,flg;
50 const PetscInt *ia,*ja;
51
52 ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;
53 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
54 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
55
56 /* Get size of fine grids and coarse grids */
57 user.ratio = 2;
58 user.coarse.mx = 4; user.coarse.my = 4; user.coarse.mz = 4;
59
60 ierr = PetscOptionsGetInt(NULL,NULL,"-Mx",&user.coarse.mx,NULL);CHKERRQ(ierr);
61 ierr = PetscOptionsGetInt(NULL,NULL,"-My",&user.coarse.my,NULL);CHKERRQ(ierr);
62 ierr = PetscOptionsGetInt(NULL,NULL,"-Mz",&user.coarse.mz,NULL);CHKERRQ(ierr);
63 ierr = PetscOptionsGetInt(NULL,NULL,"-ratio",&user.ratio,NULL);CHKERRQ(ierr);
64 if (user.coarse.mz) Test_3D = PETSC_TRUE;
65
66 user.fine.mx = user.ratio*(user.coarse.mx-1)+1;
67 user.fine.my = user.ratio*(user.coarse.my-1)+1;
68 user.fine.mz = user.ratio*(user.coarse.mz-1)+1;
69
70 if (!rank) {
71 if (!Test_3D) {
72 ierr = PetscPrintf(PETSC_COMM_SELF,"coarse grids: %D %D; fine grids: %D %D\n",user.coarse.mx,user.coarse.my,user.fine.mx,user.fine.my);CHKERRQ(ierr);
73 } else {
74 ierr = PetscPrintf(PETSC_COMM_SELF,"coarse grids: %D %D %D; fine grids: %D %D %D\n",user.coarse.mx,user.coarse.my,user.coarse.mz,user.fine.mx,user.fine.my,user.fine.mz);CHKERRQ(ierr);
75 }
76 }
77
78 /* Set up distributed array for fine grid */
79 if (!Test_3D) {
80 ierr = DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,user.fine.mx,user.fine.my,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&user.fine.da);CHKERRQ(ierr);
81 } else {
82 ierr = DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,user.fine.mx,user.fine.my,user.fine.mz,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,
83 1,1,NULL,NULL,NULL,&user.fine.da);CHKERRQ(ierr);
84 }
85 ierr = DMSetFromOptions(user.fine.da);CHKERRQ(ierr);
86 ierr = DMSetUp(user.fine.da);CHKERRQ(ierr);
87
88 /* Create and set A at fine grids */
89 ierr = DMSetMatType(user.fine.da,MATAIJ);CHKERRQ(ierr);
90 ierr = DMCreateMatrix(user.fine.da,&A);CHKERRQ(ierr);
91 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr);
92 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr);
93
94 /* set val=one to A (replace with random values!) */
95 ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rdm);CHKERRQ(ierr);
96 ierr = PetscRandomSetFromOptions(rdm);CHKERRQ(ierr);
97 if (size == 1) {
98 ierr = MatGetRowIJ(A,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);CHKERRQ(ierr);
99 if (flg) {
100 ierr = MatSeqAIJGetArray(A,&array);CHKERRQ(ierr);
101 for (i=0; i<ia[nrows]; i++) array[i] = one;
102 ierr = MatSeqAIJRestoreArray(A,&array);CHKERRQ(ierr);
103 }
104 ierr = MatRestoreRowIJ(A,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);CHKERRQ(ierr);
105 } else {
106 Mat AA,AB;
107 ierr = MatMPIAIJGetSeqAIJ(A,&AA,&AB,NULL);CHKERRQ(ierr);
108 ierr = MatGetRowIJ(AA,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);CHKERRQ(ierr);
109 if (flg) {
110 ierr = MatSeqAIJGetArray(AA,&array);CHKERRQ(ierr);
111 for (i=0; i<ia[nrows]; i++) array[i] = one;
112 ierr = MatSeqAIJRestoreArray(AA,&array);CHKERRQ(ierr);
113 }
114 ierr = MatRestoreRowIJ(AA,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);CHKERRQ(ierr);
115 ierr = MatGetRowIJ(AB,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);CHKERRQ(ierr);
116 if (flg) {
117 ierr = MatSeqAIJGetArray(AB,&array);CHKERRQ(ierr);
118 for (i=0; i<ia[nrows]; i++) array[i] = one;
119 ierr = MatSeqAIJRestoreArray(AB,&array);CHKERRQ(ierr);
120 }
121 ierr = MatRestoreRowIJ(AB,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);CHKERRQ(ierr);
122 }
123 /* Set up distributed array for coarse grid */
124 if (!Test_3D) {
125 ierr = DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,user.coarse.mx,user.coarse.my,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&user.coarse.da);CHKERRQ(ierr);
126 } else {
127 ierr = DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,user.coarse.mx,user.coarse.my,user.coarse.mz,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,NULL,&user.coarse.da);CHKERRQ(ierr);
128 }
129 ierr = DMSetFromOptions(user.coarse.da);CHKERRQ(ierr);
130 ierr = DMSetUp(user.coarse.da);CHKERRQ(ierr);
131
132 /* Create interpolation between the fine and coarse grids */
133 ierr = DMCreateInterpolation(user.coarse.da,user.fine.da,&P,NULL);CHKERRQ(ierr);
134
135 /* Get R = P^T */
136 ierr = MatTranspose(P,MAT_INITIAL_MATRIX,&R);CHKERRQ(ierr);
137
138 /* C = R*A*P */
139 /* Developer's API */
140 ierr = MatProductCreate(R,A,P,&D);CHKERRQ(ierr);
141 ierr = MatProductSetType(D,MATPRODUCT_ABC);CHKERRQ(ierr);
142 ierr = MatProductSetFromOptions(D);CHKERRQ(ierr);
143 ierr = MatProductSymbolic(D);CHKERRQ(ierr);
144 ierr = MatProductNumeric(D);CHKERRQ(ierr);
145 ierr = MatProductNumeric(D);CHKERRQ(ierr); /* Test reuse symbolic D */
146
147 /* User's API */
148 { /* Test MatMatMatMult_Basic() */
149 Mat Adense,Cdense;
150 ierr = MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&Adense);CHKERRQ(ierr);
151 ierr = MatMatMatMult(R,Adense,P,MAT_INITIAL_MATRIX,fill,&Cdense);CHKERRQ(ierr);
152 ierr = MatMatMatMult(R,Adense,P,MAT_REUSE_MATRIX,fill,&Cdense);CHKERRQ(ierr);
153
154 ierr = MatMultEqual(D,Cdense,10,&flg);CHKERRQ(ierr);
155 if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"D*v != Cdense*v");
156 ierr = MatDestroy(&Adense);CHKERRQ(ierr);
157 ierr = MatDestroy(&Cdense);CHKERRQ(ierr);
158 }
159
160 ierr = MatMatMatMult(R,A,P,MAT_INITIAL_MATRIX,fill,&C);CHKERRQ(ierr);
161 ierr = MatMatMatMult(R,A,P,MAT_REUSE_MATRIX,fill,&C);CHKERRQ(ierr);
162 ierr = MatProductClear(C);CHKERRQ(ierr);
163
164 /* Test D == C */
165 ierr = MatEqual(D,C,&flg);CHKERRQ(ierr);
166 if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"D != C");
167
168 /* Test C == PtAP */
169 ierr = MatPtAP(A,P,MAT_INITIAL_MATRIX,fill,&PtAP);CHKERRQ(ierr);
170 ierr = MatPtAP(A,P,MAT_REUSE_MATRIX,fill,&PtAP);CHKERRQ(ierr);
171 ierr = MatEqual(C,PtAP,&flg);CHKERRQ(ierr);
172 if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"C != PtAP");
173 ierr = MatDestroy(&PtAP);CHKERRQ(ierr);
174
175 /* Clean up */
176 ierr = MatDestroy(&A);CHKERRQ(ierr);
177 ierr = PetscRandomDestroy(&rdm);CHKERRQ(ierr);
178 ierr = DMDestroy(&user.fine.da);CHKERRQ(ierr);
179 ierr = DMDestroy(&user.coarse.da);CHKERRQ(ierr);
180 ierr = MatDestroy(&P);CHKERRQ(ierr);
181 ierr = MatDestroy(&R);CHKERRQ(ierr);
182 ierr = MatDestroy(&C);CHKERRQ(ierr);
183 ierr = MatDestroy(&D);CHKERRQ(ierr);
184 ierr = PetscFinalize();
185 return ierr;
186 }
187
188 /*TEST
189
190 test:
191
192 test:
193 suffix: 2
194 nsize: 2
195 args: -matmatmatmult_via scalable
196
197 test:
198 suffix: 3
199 nsize: 2
200 args: -matmatmatmult_via nonscalable
201 output_file: output/ex111_1.out
202
203 TEST*/
204