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 */ 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; 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 ierr = MatMatMatMult(R,A,P,MAT_INITIAL_MATRIX,fill,&C);CHKERRQ(ierr); 140 ierr = MatMatMatMult(R,A,P,MAT_REUSE_MATRIX,fill,&C);CHKERRQ(ierr); 141 ierr = MatFreeIntermediateDataStructures(C);CHKERRQ(ierr); 142 143 /* Test C == PtAP */ 144 ierr = MatPtAP(A,P,MAT_INITIAL_MATRIX,fill,&PtAP);CHKERRQ(ierr); 145 ierr = MatPtAP(A,P,MAT_REUSE_MATRIX,fill,&PtAP);CHKERRQ(ierr); 146 ierr = MatEqual(C,PtAP,&flg);CHKERRQ(ierr); 147 if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Matrices are not equal"); 148 ierr = MatDestroy(&PtAP);CHKERRQ(ierr); 149 150 /* Clean up */ 151 ierr = MatDestroy(&A);CHKERRQ(ierr); 152 ierr = PetscRandomDestroy(&rdm);CHKERRQ(ierr); 153 ierr = DMDestroy(&user.fine.da);CHKERRQ(ierr); 154 ierr = DMDestroy(&user.coarse.da);CHKERRQ(ierr); 155 ierr = MatDestroy(&P);CHKERRQ(ierr); 156 ierr = MatDestroy(&R);CHKERRQ(ierr); 157 ierr = MatDestroy(&C);CHKERRQ(ierr); 158 ierr = PetscFinalize(); 159 return ierr; 160 } 161 162 /*TEST 163 164 test: 165 166 test: 167 suffix: 2 168 nsize: 2 169 args: -matmatmatmult_via scalable 170 171 test: 172 suffix: 3 173 nsize: 2 174 args: -matmatmatmult_via nonscalable 175 output_file: output/ex111_1.out 176 177 TEST*/ 178