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