1 2 static char help[] ="Tests sequential and parallel DMCreateMatrix(), MatMatMult() and MatPtAP()\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 This test is modified from ~src/ksp/tests/ex19.c. 12 Example of usage: mpiexec -n 3 ./ex96 -Mx 10 -My 10 -Mz 10 13 */ 14 15 #include <petscdm.h> 16 #include <petscdmda.h> 17 18 /* User-defined application contexts */ 19 typedef struct { 20 PetscInt mx,my,mz; /* number grid points in x, y and z direction */ 21 Vec localX,localF; /* local vectors with ghost region */ 22 DM da; 23 Vec x,b,r; /* global vectors */ 24 Mat J; /* Jacobian on grid */ 25 } GridCtx; 26 typedef struct { 27 GridCtx fine; 28 GridCtx coarse; 29 PetscInt ratio; 30 Mat Ii; /* interpolation from coarse to fine */ 31 } AppCtx; 32 33 #define COARSE_LEVEL 0 34 #define FINE_LEVEL 1 35 36 /* 37 Mm_ratio - ration of grid lines between fine and coarse grids. 38 */ 39 int main(int argc,char **argv) 40 { 41 PetscErrorCode ierr; 42 AppCtx user; 43 PetscInt Npx=PETSC_DECIDE,Npy=PETSC_DECIDE,Npz=PETSC_DECIDE; 44 PetscMPIInt size,rank; 45 PetscInt m,n,M,N,i,nrows; 46 PetscScalar one = 1.0; 47 PetscReal fill=2.0; 48 Mat A,A_tmp,P,C,C1,C2; 49 PetscScalar *array,none = -1.0,alpha; 50 Vec x,v1,v2,v3,v4; 51 PetscReal norm,norm_tmp,norm_tmp1,tol=100.*PETSC_MACHINE_EPSILON; 52 PetscRandom rdm; 53 PetscBool Test_MatMatMult=PETSC_TRUE,Test_MatPtAP=PETSC_TRUE,Test_3D=PETSC_TRUE,flg; 54 const PetscInt *ia,*ja; 55 56 ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr; 57 ierr = PetscOptionsGetReal(NULL,NULL,"-tol",&tol,NULL);CHKERRQ(ierr); 58 59 user.ratio = 2; 60 user.coarse.mx = 20; user.coarse.my = 20; user.coarse.mz = 20; 61 62 ierr = PetscOptionsGetInt(NULL,NULL,"-Mx",&user.coarse.mx,NULL);CHKERRQ(ierr); 63 ierr = PetscOptionsGetInt(NULL,NULL,"-My",&user.coarse.my,NULL);CHKERRQ(ierr); 64 ierr = PetscOptionsGetInt(NULL,NULL,"-Mz",&user.coarse.mz,NULL);CHKERRQ(ierr); 65 ierr = PetscOptionsGetInt(NULL,NULL,"-ratio",&user.ratio,NULL);CHKERRQ(ierr); 66 67 if (user.coarse.mz) Test_3D = PETSC_TRUE; 68 69 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); 70 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRMPI(ierr); 71 ierr = PetscOptionsGetInt(NULL,NULL,"-Npx",&Npx,NULL);CHKERRQ(ierr); 72 ierr = PetscOptionsGetInt(NULL,NULL,"-Npy",&Npy,NULL);CHKERRQ(ierr); 73 ierr = PetscOptionsGetInt(NULL,NULL,"-Npz",&Npz,NULL);CHKERRQ(ierr); 74 75 /* Set up distributed array for fine grid */ 76 if (!Test_3D) { 77 ierr = DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,user.coarse.mx,user.coarse.my,Npx,Npy,1,1,NULL,NULL,&user.coarse.da);CHKERRQ(ierr); 78 } else { 79 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,Npx,Npy,Npz,1,1,NULL,NULL,NULL,&user.coarse.da);CHKERRQ(ierr); 80 } 81 ierr = DMSetFromOptions(user.coarse.da);CHKERRQ(ierr); 82 ierr = DMSetUp(user.coarse.da);CHKERRQ(ierr); 83 84 /* This makes sure the coarse DMDA has the same partition as the fine DMDA */ 85 ierr = DMRefine(user.coarse.da,PetscObjectComm((PetscObject)user.coarse.da),&user.fine.da);CHKERRQ(ierr); 86 87 /* Test DMCreateMatrix() */ 88 /*------------------------------------------------------------*/ 89 ierr = DMSetMatType(user.fine.da,MATAIJ);CHKERRQ(ierr); 90 ierr = DMCreateMatrix(user.fine.da,&A);CHKERRQ(ierr); 91 ierr = DMSetMatType(user.fine.da,MATBAIJ);CHKERRQ(ierr); 92 ierr = DMCreateMatrix(user.fine.da,&C);CHKERRQ(ierr); 93 94 ierr = MatConvert(C,MATAIJ,MAT_INITIAL_MATRIX,&A_tmp);CHKERRQ(ierr); /* not work for mpisbaij matrix! */ 95 ierr = MatEqual(A,A_tmp,&flg);CHKERRQ(ierr); 96 if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"A != C"); 97 ierr = MatDestroy(&C);CHKERRQ(ierr); 98 ierr = MatDestroy(&A_tmp);CHKERRQ(ierr); 99 100 /*------------------------------------------------------------*/ 101 102 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr); 103 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr); 104 /* if (!rank) printf("A %d, %d\n",M,N); */ 105 106 /* set val=one to A */ 107 if (size == 1) { 108 ierr = MatGetRowIJ(A,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);CHKERRQ(ierr); 109 if (flg) { 110 ierr = MatSeqAIJGetArray(A,&array);CHKERRQ(ierr); 111 for (i=0; i<ia[nrows]; i++) array[i] = one; 112 ierr = MatSeqAIJRestoreArray(A,&array);CHKERRQ(ierr); 113 } 114 ierr = MatRestoreRowIJ(A,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);CHKERRQ(ierr); 115 } else { 116 Mat AA,AB; 117 ierr = MatMPIAIJGetSeqAIJ(A,&AA,&AB,NULL);CHKERRQ(ierr); 118 ierr = MatGetRowIJ(AA,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);CHKERRQ(ierr); 119 if (flg) { 120 ierr = MatSeqAIJGetArray(AA,&array);CHKERRQ(ierr); 121 for (i=0; i<ia[nrows]; i++) array[i] = one; 122 ierr = MatSeqAIJRestoreArray(AA,&array);CHKERRQ(ierr); 123 } 124 ierr = MatRestoreRowIJ(AA,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);CHKERRQ(ierr); 125 ierr = MatGetRowIJ(AB,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);CHKERRQ(ierr); 126 if (flg) { 127 ierr = MatSeqAIJGetArray(AB,&array);CHKERRQ(ierr); 128 for (i=0; i<ia[nrows]; i++) array[i] = one; 129 ierr = MatSeqAIJRestoreArray(AB,&array);CHKERRQ(ierr); 130 } 131 ierr = MatRestoreRowIJ(AB,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);CHKERRQ(ierr); 132 } 133 /* ierr = MatView(A, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */ 134 135 /* Create interpolation between the fine and coarse grids */ 136 ierr = DMCreateInterpolation(user.coarse.da,user.fine.da,&P,NULL);CHKERRQ(ierr); 137 ierr = MatGetLocalSize(P,&m,&n);CHKERRQ(ierr); 138 ierr = MatGetSize(P,&M,&N);CHKERRQ(ierr); 139 /* if (!rank) printf("P %d, %d\n",M,N); */ 140 141 /* Create vectors v1 and v2 that are compatible with A */ 142 ierr = VecCreate(PETSC_COMM_WORLD,&v1);CHKERRQ(ierr); 143 ierr = MatGetLocalSize(A,&m,NULL);CHKERRQ(ierr); 144 ierr = VecSetSizes(v1,m,PETSC_DECIDE);CHKERRQ(ierr); 145 ierr = VecSetFromOptions(v1);CHKERRQ(ierr); 146 ierr = VecDuplicate(v1,&v2);CHKERRQ(ierr); 147 ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rdm);CHKERRQ(ierr); 148 ierr = PetscRandomSetFromOptions(rdm);CHKERRQ(ierr); 149 150 /* Test MatMatMult(): C = A*P */ 151 /*----------------------------*/ 152 if (Test_MatMatMult) { 153 ierr = MatDuplicate(A,MAT_COPY_VALUES,&A_tmp);CHKERRQ(ierr); 154 ierr = MatMatMult(A_tmp,P,MAT_INITIAL_MATRIX,fill,&C);CHKERRQ(ierr); 155 156 /* Test MAT_REUSE_MATRIX - reuse symbolic C */ 157 alpha=1.0; 158 for (i=0; i<2; i++) { 159 alpha -= 0.1; 160 ierr = MatScale(A_tmp,alpha);CHKERRQ(ierr); 161 ierr = MatMatMult(A_tmp,P,MAT_REUSE_MATRIX,fill,&C);CHKERRQ(ierr); 162 } 163 /* Free intermediate data structures created for reuse of C=Pt*A*P */ 164 ierr = MatProductClear(C);CHKERRQ(ierr); 165 166 /* Test MatDuplicate() */ 167 /*----------------------------*/ 168 ierr = MatDuplicate(C,MAT_COPY_VALUES,&C1);CHKERRQ(ierr); 169 ierr = MatDuplicate(C1,MAT_COPY_VALUES,&C2);CHKERRQ(ierr); 170 ierr = MatDestroy(&C1);CHKERRQ(ierr); 171 ierr = MatDestroy(&C2);CHKERRQ(ierr); 172 173 /* Create vector x that is compatible with P */ 174 ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr); 175 ierr = MatGetLocalSize(P,NULL,&n);CHKERRQ(ierr); 176 ierr = VecSetSizes(x,n,PETSC_DECIDE);CHKERRQ(ierr); 177 ierr = VecSetFromOptions(x);CHKERRQ(ierr); 178 179 norm = 0.0; 180 for (i=0; i<10; i++) { 181 ierr = VecSetRandom(x,rdm);CHKERRQ(ierr); 182 ierr = MatMult(P,x,v1);CHKERRQ(ierr); 183 ierr = MatMult(A_tmp,v1,v2);CHKERRQ(ierr); /* v2 = A*P*x */ 184 ierr = MatMult(C,x,v1);CHKERRQ(ierr); /* v1 = C*x */ 185 ierr = VecAXPY(v1,none,v2);CHKERRQ(ierr); 186 ierr = VecNorm(v1,NORM_1,&norm_tmp);CHKERRQ(ierr); 187 ierr = VecNorm(v2,NORM_1,&norm_tmp1);CHKERRQ(ierr); 188 norm_tmp /= norm_tmp1; 189 if (norm_tmp > norm) norm = norm_tmp; 190 } 191 if (norm >= tol && !rank) { 192 ierr = PetscPrintf(PETSC_COMM_SELF,"Error: MatMatMult(), |v1 - v2|/|v2|: %g\n",(double)norm);CHKERRQ(ierr); 193 } 194 195 ierr = VecDestroy(&x);CHKERRQ(ierr); 196 ierr = MatDestroy(&C);CHKERRQ(ierr); 197 ierr = MatDestroy(&A_tmp);CHKERRQ(ierr); 198 } 199 200 /* Test P^T * A * P - MatPtAP() */ 201 /*------------------------------*/ 202 if (Test_MatPtAP) { 203 ierr = MatPtAP(A,P,MAT_INITIAL_MATRIX,fill,&C);CHKERRQ(ierr); 204 ierr = MatGetLocalSize(C,&m,&n);CHKERRQ(ierr); 205 206 /* Test MAT_REUSE_MATRIX - reuse symbolic C */ 207 alpha=1.0; 208 for (i=0; i<1; i++) { 209 alpha -= 0.1; 210 ierr = MatScale(A,alpha);CHKERRQ(ierr); 211 ierr = MatPtAP(A,P,MAT_REUSE_MATRIX,fill,&C);CHKERRQ(ierr); 212 } 213 214 /* Free intermediate data structures created for reuse of C=Pt*A*P */ 215 ierr = MatProductClear(C);CHKERRQ(ierr); 216 217 /* Test MatDuplicate() */ 218 /*----------------------------*/ 219 ierr = MatDuplicate(C,MAT_COPY_VALUES,&C1);CHKERRQ(ierr); 220 ierr = MatDuplicate(C1,MAT_COPY_VALUES,&C2);CHKERRQ(ierr); 221 ierr = MatDestroy(&C1);CHKERRQ(ierr); 222 ierr = MatDestroy(&C2);CHKERRQ(ierr); 223 224 /* Create vector x that is compatible with P */ 225 ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr); 226 ierr = MatGetLocalSize(P,&m,&n);CHKERRQ(ierr); 227 ierr = VecSetSizes(x,n,PETSC_DECIDE);CHKERRQ(ierr); 228 ierr = VecSetFromOptions(x);CHKERRQ(ierr); 229 230 ierr = VecCreate(PETSC_COMM_WORLD,&v3);CHKERRQ(ierr); 231 ierr = VecSetSizes(v3,n,PETSC_DECIDE);CHKERRQ(ierr); 232 ierr = VecSetFromOptions(v3);CHKERRQ(ierr); 233 ierr = VecDuplicate(v3,&v4);CHKERRQ(ierr); 234 235 norm = 0.0; 236 for (i=0; i<10; i++) { 237 ierr = VecSetRandom(x,rdm);CHKERRQ(ierr); 238 ierr = MatMult(P,x,v1);CHKERRQ(ierr); 239 ierr = MatMult(A,v1,v2);CHKERRQ(ierr); /* v2 = A*P*x */ 240 241 ierr = MatMultTranspose(P,v2,v3);CHKERRQ(ierr); /* v3 = Pt*A*P*x */ 242 ierr = MatMult(C,x,v4);CHKERRQ(ierr); /* v3 = C*x */ 243 ierr = VecAXPY(v4,none,v3);CHKERRQ(ierr); 244 ierr = VecNorm(v4,NORM_1,&norm_tmp);CHKERRQ(ierr); 245 ierr = VecNorm(v3,NORM_1,&norm_tmp1);CHKERRQ(ierr); 246 247 norm_tmp /= norm_tmp1; 248 if (norm_tmp > norm) norm = norm_tmp; 249 } 250 if (norm >= tol && !rank) { 251 ierr = PetscPrintf(PETSC_COMM_SELF,"Error: MatPtAP(), |v3 - v4|/|v3|: %g\n",(double)norm);CHKERRQ(ierr); 252 } 253 ierr = MatDestroy(&C);CHKERRQ(ierr); 254 ierr = VecDestroy(&v3);CHKERRQ(ierr); 255 ierr = VecDestroy(&v4);CHKERRQ(ierr); 256 ierr = VecDestroy(&x);CHKERRQ(ierr); 257 } 258 259 /* Clean up */ 260 ierr = MatDestroy(&A);CHKERRQ(ierr); 261 ierr = PetscRandomDestroy(&rdm);CHKERRQ(ierr); 262 ierr = VecDestroy(&v1);CHKERRQ(ierr); 263 ierr = VecDestroy(&v2);CHKERRQ(ierr); 264 ierr = DMDestroy(&user.fine.da);CHKERRQ(ierr); 265 ierr = DMDestroy(&user.coarse.da);CHKERRQ(ierr); 266 ierr = MatDestroy(&P);CHKERRQ(ierr); 267 ierr = PetscFinalize(); 268 return ierr; 269 } 270 271 /*TEST 272 273 test: 274 args: -Mx 10 -My 5 -Mz 10 275 output_file: output/ex96_1.out 276 277 test: 278 suffix: nonscalable 279 nsize: 3 280 args: -Mx 10 -My 5 -Mz 10 281 output_file: output/ex96_1.out 282 283 test: 284 suffix: scalable 285 nsize: 3 286 args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via scalable 287 output_file: output/ex96_1.out 288 289 test: 290 suffix: seq_scalable 291 nsize: 3 292 args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via scalable -inner_diag_matproduct_ab_via scalable -inner_offdiag_matproduct_ab_via scalable 293 output_file: output/ex96_1.out 294 295 test: 296 suffix: seq_sorted 297 nsize: 3 298 args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via scalable -inner_diag_matproduct_ab_via sorted -inner_offdiag_matproduct_ab_via sorted 299 output_file: output/ex96_1.out 300 301 test: 302 suffix: seq_scalable_fast 303 nsize: 3 304 args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via scalable -inner_diag_matproduct_ab_via scalable_fast -inner_offdiag_matproduct_ab_via scalable_fast 305 output_file: output/ex96_1.out 306 307 test: 308 suffix: seq_heap 309 nsize: 3 310 args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via scalable -inner_diag_matproduct_ab_via heap -inner_offdiag_matproduct_ab_via heap 311 output_file: output/ex96_1.out 312 313 test: 314 suffix: seq_btheap 315 nsize: 3 316 args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via scalable -inner_diag_matproduct_ab_via btheap -inner_offdiag_matproduct_ab_via btheap 317 output_file: output/ex96_1.out 318 319 test: 320 suffix: seq_llcondensed 321 nsize: 3 322 args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via scalable -inner_diag_matproduct_ab_via llcondensed -inner_offdiag_matproduct_ab_via llcondensed 323 output_file: output/ex96_1.out 324 325 test: 326 suffix: seq_rowmerge 327 nsize: 3 328 args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via scalable -inner_diag_matproduct_ab_via rowmerge -inner_offdiag_matproduct_ab_via rowmerge 329 output_file: output/ex96_1.out 330 331 test: 332 suffix: allatonce 333 nsize: 3 334 args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via allatonce 335 output_file: output/ex96_1.out 336 337 test: 338 suffix: allatonce_merged 339 nsize: 3 340 args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via allatonce_merged 341 output_file: output/ex96_1.out 342 343 TEST*/ 344