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 AppCtx user; 42 PetscInt Npx=PETSC_DECIDE,Npy=PETSC_DECIDE,Npz=PETSC_DECIDE; 43 PetscMPIInt size,rank; 44 PetscInt m,n,M,N,i,nrows; 45 PetscScalar one = 1.0; 46 PetscReal fill=2.0; 47 Mat A,A_tmp,P,C,C1,C2; 48 PetscScalar *array,none = -1.0,alpha; 49 Vec x,v1,v2,v3,v4; 50 PetscReal norm,norm_tmp,norm_tmp1,tol=100.*PETSC_MACHINE_EPSILON; 51 PetscRandom rdm; 52 PetscBool Test_MatMatMult=PETSC_TRUE,Test_MatPtAP=PETSC_TRUE,Test_3D=PETSC_TRUE,flg; 53 const PetscInt *ia,*ja; 54 55 PetscFunctionBeginUser; 56 PetscCall(PetscInitialize(&argc,&argv,NULL,help)); 57 PetscCall(PetscOptionsGetReal(NULL,NULL,"-tol",&tol,NULL)); 58 59 user.ratio = 2; 60 user.coarse.mx = 20; user.coarse.my = 20; user.coarse.mz = 20; 61 62 PetscCall(PetscOptionsGetInt(NULL,NULL,"-Mx",&user.coarse.mx,NULL)); 63 PetscCall(PetscOptionsGetInt(NULL,NULL,"-My",&user.coarse.my,NULL)); 64 PetscCall(PetscOptionsGetInt(NULL,NULL,"-Mz",&user.coarse.mz,NULL)); 65 PetscCall(PetscOptionsGetInt(NULL,NULL,"-ratio",&user.ratio,NULL)); 66 67 if (user.coarse.mz) Test_3D = PETSC_TRUE; 68 69 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 70 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 71 PetscCall(PetscOptionsGetInt(NULL,NULL,"-Npx",&Npx,NULL)); 72 PetscCall(PetscOptionsGetInt(NULL,NULL,"-Npy",&Npy,NULL)); 73 PetscCall(PetscOptionsGetInt(NULL,NULL,"-Npz",&Npz,NULL)); 74 75 /* Set up distributed array for fine grid */ 76 if (!Test_3D) { 77 PetscCall(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)); 78 } else { 79 PetscCall(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)); 80 } 81 PetscCall(DMSetFromOptions(user.coarse.da)); 82 PetscCall(DMSetUp(user.coarse.da)); 83 84 /* This makes sure the coarse DMDA has the same partition as the fine DMDA */ 85 PetscCall(DMRefine(user.coarse.da,PetscObjectComm((PetscObject)user.coarse.da),&user.fine.da)); 86 87 /* Test DMCreateMatrix() */ 88 /*------------------------------------------------------------*/ 89 PetscCall(DMSetMatType(user.fine.da,MATAIJ)); 90 PetscCall(DMCreateMatrix(user.fine.da,&A)); 91 PetscCall(DMSetMatType(user.fine.da,MATBAIJ)); 92 PetscCall(DMCreateMatrix(user.fine.da,&C)); 93 94 PetscCall(MatConvert(C,MATAIJ,MAT_INITIAL_MATRIX,&A_tmp)); /* not work for mpisbaij matrix! */ 95 PetscCall(MatEqual(A,A_tmp,&flg)); 96 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"A != C"); 97 PetscCall(MatDestroy(&C)); 98 PetscCall(MatDestroy(&A_tmp)); 99 100 /*------------------------------------------------------------*/ 101 102 PetscCall(MatGetLocalSize(A,&m,&n)); 103 PetscCall(MatGetSize(A,&M,&N)); 104 /* if (rank == 0) printf("A %d, %d\n",M,N); */ 105 106 /* set val=one to A */ 107 if (size == 1) { 108 PetscCall(MatGetRowIJ(A,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg)); 109 if (flg) { 110 PetscCall(MatSeqAIJGetArray(A,&array)); 111 for (i=0; i<ia[nrows]; i++) array[i] = one; 112 PetscCall(MatSeqAIJRestoreArray(A,&array)); 113 } 114 PetscCall(MatRestoreRowIJ(A,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg)); 115 } else { 116 Mat AA,AB; 117 PetscCall(MatMPIAIJGetSeqAIJ(A,&AA,&AB,NULL)); 118 PetscCall(MatGetRowIJ(AA,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg)); 119 if (flg) { 120 PetscCall(MatSeqAIJGetArray(AA,&array)); 121 for (i=0; i<ia[nrows]; i++) array[i] = one; 122 PetscCall(MatSeqAIJRestoreArray(AA,&array)); 123 } 124 PetscCall(MatRestoreRowIJ(AA,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg)); 125 PetscCall(MatGetRowIJ(AB,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg)); 126 if (flg) { 127 PetscCall(MatSeqAIJGetArray(AB,&array)); 128 for (i=0; i<ia[nrows]; i++) array[i] = one; 129 PetscCall(MatSeqAIJRestoreArray(AB,&array)); 130 } 131 PetscCall(MatRestoreRowIJ(AB,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg)); 132 } 133 /* PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD)); */ 134 135 /* Create interpolation between the fine and coarse grids */ 136 PetscCall(DMCreateInterpolation(user.coarse.da,user.fine.da,&P,NULL)); 137 PetscCall(MatGetLocalSize(P,&m,&n)); 138 PetscCall(MatGetSize(P,&M,&N)); 139 /* if (rank == 0) printf("P %d, %d\n",M,N); */ 140 141 /* Create vectors v1 and v2 that are compatible with A */ 142 PetscCall(VecCreate(PETSC_COMM_WORLD,&v1)); 143 PetscCall(MatGetLocalSize(A,&m,NULL)); 144 PetscCall(VecSetSizes(v1,m,PETSC_DECIDE)); 145 PetscCall(VecSetFromOptions(v1)); 146 PetscCall(VecDuplicate(v1,&v2)); 147 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD,&rdm)); 148 PetscCall(PetscRandomSetFromOptions(rdm)); 149 150 /* Test MatMatMult(): C = A*P */ 151 /*----------------------------*/ 152 if (Test_MatMatMult) { 153 PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&A_tmp)); 154 PetscCall(MatMatMult(A_tmp,P,MAT_INITIAL_MATRIX,fill,&C)); 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 PetscCall(MatScale(A_tmp,alpha)); 161 PetscCall(MatMatMult(A_tmp,P,MAT_REUSE_MATRIX,fill,&C)); 162 } 163 /* Free intermediate data structures created for reuse of C=Pt*A*P */ 164 PetscCall(MatProductClear(C)); 165 166 /* Test MatDuplicate() */ 167 /*----------------------------*/ 168 PetscCall(MatDuplicate(C,MAT_COPY_VALUES,&C1)); 169 PetscCall(MatDuplicate(C1,MAT_COPY_VALUES,&C2)); 170 PetscCall(MatDestroy(&C1)); 171 PetscCall(MatDestroy(&C2)); 172 173 /* Create vector x that is compatible with P */ 174 PetscCall(VecCreate(PETSC_COMM_WORLD,&x)); 175 PetscCall(MatGetLocalSize(P,NULL,&n)); 176 PetscCall(VecSetSizes(x,n,PETSC_DECIDE)); 177 PetscCall(VecSetFromOptions(x)); 178 179 norm = 0.0; 180 for (i=0; i<10; i++) { 181 PetscCall(VecSetRandom(x,rdm)); 182 PetscCall(MatMult(P,x,v1)); 183 PetscCall(MatMult(A_tmp,v1,v2)); /* v2 = A*P*x */ 184 PetscCall(MatMult(C,x,v1)); /* v1 = C*x */ 185 PetscCall(VecAXPY(v1,none,v2)); 186 PetscCall(VecNorm(v1,NORM_1,&norm_tmp)); 187 PetscCall(VecNorm(v2,NORM_1,&norm_tmp1)); 188 norm_tmp /= norm_tmp1; 189 if (norm_tmp > norm) norm = norm_tmp; 190 } 191 if (norm >= tol && rank == 0) { 192 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error: MatMatMult(), |v1 - v2|/|v2|: %g\n",(double)norm)); 193 } 194 195 PetscCall(VecDestroy(&x)); 196 PetscCall(MatDestroy(&C)); 197 PetscCall(MatDestroy(&A_tmp)); 198 } 199 200 /* Test P^T * A * P - MatPtAP() */ 201 /*------------------------------*/ 202 if (Test_MatPtAP) { 203 PetscCall(MatPtAP(A,P,MAT_INITIAL_MATRIX,fill,&C)); 204 PetscCall(MatGetLocalSize(C,&m,&n)); 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 PetscCall(MatScale(A,alpha)); 211 PetscCall(MatPtAP(A,P,MAT_REUSE_MATRIX,fill,&C)); 212 } 213 214 /* Free intermediate data structures created for reuse of C=Pt*A*P */ 215 PetscCall(MatProductClear(C)); 216 217 /* Test MatDuplicate() */ 218 /*----------------------------*/ 219 PetscCall(MatDuplicate(C,MAT_COPY_VALUES,&C1)); 220 PetscCall(MatDuplicate(C1,MAT_COPY_VALUES,&C2)); 221 PetscCall(MatDestroy(&C1)); 222 PetscCall(MatDestroy(&C2)); 223 224 /* Create vector x that is compatible with P */ 225 PetscCall(VecCreate(PETSC_COMM_WORLD,&x)); 226 PetscCall(MatGetLocalSize(P,&m,&n)); 227 PetscCall(VecSetSizes(x,n,PETSC_DECIDE)); 228 PetscCall(VecSetFromOptions(x)); 229 230 PetscCall(VecCreate(PETSC_COMM_WORLD,&v3)); 231 PetscCall(VecSetSizes(v3,n,PETSC_DECIDE)); 232 PetscCall(VecSetFromOptions(v3)); 233 PetscCall(VecDuplicate(v3,&v4)); 234 235 norm = 0.0; 236 for (i=0; i<10; i++) { 237 PetscCall(VecSetRandom(x,rdm)); 238 PetscCall(MatMult(P,x,v1)); 239 PetscCall(MatMult(A,v1,v2)); /* v2 = A*P*x */ 240 241 PetscCall(MatMultTranspose(P,v2,v3)); /* v3 = Pt*A*P*x */ 242 PetscCall(MatMult(C,x,v4)); /* v3 = C*x */ 243 PetscCall(VecAXPY(v4,none,v3)); 244 PetscCall(VecNorm(v4,NORM_1,&norm_tmp)); 245 PetscCall(VecNorm(v3,NORM_1,&norm_tmp1)); 246 247 norm_tmp /= norm_tmp1; 248 if (norm_tmp > norm) norm = norm_tmp; 249 } 250 if (norm >= tol && rank == 0) { 251 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error: MatPtAP(), |v3 - v4|/|v3|: %g\n",(double)norm)); 252 } 253 PetscCall(MatDestroy(&C)); 254 PetscCall(VecDestroy(&v3)); 255 PetscCall(VecDestroy(&v4)); 256 PetscCall(VecDestroy(&x)); 257 } 258 259 /* Clean up */ 260 PetscCall(MatDestroy(&A)); 261 PetscCall(PetscRandomDestroy(&rdm)); 262 PetscCall(VecDestroy(&v1)); 263 PetscCall(VecDestroy(&v2)); 264 PetscCall(DMDestroy(&user.fine.da)); 265 PetscCall(DMDestroy(&user.coarse.da)); 266 PetscCall(MatDestroy(&P)); 267 PetscCall(PetscFinalize()); 268 return 0; 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_mat_product_algorithm scalable -inner_offdiag_mat_product_algorithm 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_mat_product_algorithm sorted -inner_offdiag_mat_product_algorithm 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_mat_product_algorithm scalable_fast -inner_offdiag_mat_product_algorithm 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_mat_product_algorithm heap -inner_offdiag_mat_product_algorithm 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_mat_product_algorithm btheap -inner_offdiag_mat_product_algorithm 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_mat_product_algorithm llcondensed -inner_offdiag_mat_product_algorithm 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_mat_product_algorithm rowmerge -inner_offdiag_mat_product_algorithm 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