1 static char help[] ="Tests MatPtAP() for MPIMAIJ and MPIAIJ \n "; 2 3 #include <petscdmda.h> 4 5 int main(int argc,char **argv) 6 { 7 PetscErrorCode ierr; 8 DM coarsedm,finedm; 9 PetscMPIInt size,rank; 10 PetscInt M,N,Z,i,nrows; 11 PetscScalar one = 1.0; 12 PetscReal fill=2.0; 13 Mat A,P,C; 14 PetscScalar *array,alpha; 15 PetscBool Test_3D=PETSC_FALSE,flg; 16 const PetscInt *ia,*ja; 17 PetscInt dof; 18 MPI_Comm comm; 19 20 ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr; 21 comm = PETSC_COMM_WORLD; 22 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 23 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 24 M = 10; N = 10; Z = 10; 25 dof = 10; 26 27 ierr = PetscOptionsGetBool(NULL,NULL,"-test_3D",&Test_3D,NULL);CHKERRQ(ierr); 28 ierr = PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);CHKERRQ(ierr); 29 ierr = PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);CHKERRQ(ierr); 30 ierr = PetscOptionsGetInt(NULL,NULL,"-Z",&Z,NULL);CHKERRQ(ierr); 31 /* Set up distributed array for fine grid */ 32 if (!Test_3D) { 33 ierr = DMDACreate2d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,M,N,PETSC_DECIDE,PETSC_DECIDE,dof,1,NULL,NULL,&coarsedm);CHKERRQ(ierr); 34 } else { 35 ierr = DMDACreate3d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,M,N,Z,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,dof,1,NULL,NULL,NULL,&coarsedm);CHKERRQ(ierr); 36 } 37 ierr = DMSetFromOptions(coarsedm);CHKERRQ(ierr); 38 ierr = DMSetUp(coarsedm);CHKERRQ(ierr); 39 40 /* This makes sure the coarse DMDA has the same partition as the fine DMDA */ 41 ierr = DMRefine(coarsedm,PetscObjectComm((PetscObject)coarsedm),&finedm);CHKERRQ(ierr); 42 43 /*------------------------------------------------------------*/ 44 ierr = DMSetMatType(finedm,MATAIJ);CHKERRQ(ierr); 45 ierr = DMCreateMatrix(finedm,&A);CHKERRQ(ierr); 46 47 /* set val=one to A */ 48 if (size == 1) { 49 ierr = MatGetRowIJ(A,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);CHKERRQ(ierr); 50 if (flg) { 51 ierr = MatSeqAIJGetArray(A,&array);CHKERRQ(ierr); 52 for (i=0; i<ia[nrows]; i++) array[i] = one; 53 ierr = MatSeqAIJRestoreArray(A,&array);CHKERRQ(ierr); 54 } 55 ierr = MatRestoreRowIJ(A,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);CHKERRQ(ierr); 56 } else { 57 Mat AA,AB; 58 ierr = MatMPIAIJGetSeqAIJ(A,&AA,&AB,NULL);CHKERRQ(ierr); 59 ierr = MatGetRowIJ(AA,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);CHKERRQ(ierr); 60 if (flg) { 61 ierr = MatSeqAIJGetArray(AA,&array);CHKERRQ(ierr); 62 for (i=0; i<ia[nrows]; i++) array[i] = one; 63 ierr = MatSeqAIJRestoreArray(AA,&array);CHKERRQ(ierr); 64 } 65 ierr = MatRestoreRowIJ(AA,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);CHKERRQ(ierr); 66 ierr = MatGetRowIJ(AB,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);CHKERRQ(ierr); 67 if (flg) { 68 ierr = MatSeqAIJGetArray(AB,&array);CHKERRQ(ierr); 69 for (i=0; i<ia[nrows]; i++) array[i] = one; 70 ierr = MatSeqAIJRestoreArray(AB,&array);CHKERRQ(ierr); 71 } 72 ierr = MatRestoreRowIJ(AB,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);CHKERRQ(ierr); 73 } 74 /* Create interpolation between the fine and coarse grids */ 75 ierr = DMCreateInterpolation(coarsedm,finedm,&P,NULL);CHKERRQ(ierr); 76 77 /* Test P^T * A * P - MatPtAP() */ 78 /*------------------------------*/ 79 /* (1) Developer API */ 80 ierr = MatProductCreate(A,P,NULL,&C);CHKERRQ(ierr); 81 ierr = MatProductSetType(C,MATPRODUCT_PtAP);CHKERRQ(ierr); 82 ierr = MatProductSetAlgorithm(C,"default");CHKERRQ(ierr); 83 ierr = MatProductSetFill(C,PETSC_DEFAULT);CHKERRQ(ierr); 84 ierr = MatProductSetFromOptions(C);CHKERRQ(ierr); 85 ierr = MatProductSymbolic(C);CHKERRQ(ierr); 86 ierr = MatProductNumeric(C);CHKERRQ(ierr); 87 ierr = MatProductNumeric(C);CHKERRQ(ierr); /* Test reuse of symbolic C */ 88 89 ierr = MatPtAPMultEqual(A,P,C,10,&flg);CHKERRQ(ierr); 90 if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatProduct_PtAP"); 91 ierr = MatDestroy(&C);CHKERRQ(ierr); 92 93 /* (2) User API */ 94 ierr = MatPtAP(A,P,MAT_INITIAL_MATRIX,fill,&C);CHKERRQ(ierr); 95 /* Test MAT_REUSE_MATRIX - reuse symbolic C */ 96 alpha=1.0; 97 for (i=0; i<1; i++) { 98 alpha -= 0.1; 99 ierr = MatScale(A,alpha);CHKERRQ(ierr); 100 ierr = MatPtAP(A,P,MAT_REUSE_MATRIX,fill,&C);CHKERRQ(ierr); 101 } 102 103 /* Free intermediate data structures created for reuse of C=Pt*A*P */ 104 ierr = MatProductClear(C);CHKERRQ(ierr); 105 106 ierr = MatPtAPMultEqual(A,P,C,10,&flg);CHKERRQ(ierr); 107 if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatPtAP"); 108 109 ierr = MatDestroy(&C);CHKERRQ(ierr); 110 ierr = MatDestroy(&A);CHKERRQ(ierr); 111 ierr = MatDestroy(&P);CHKERRQ(ierr); 112 ierr = DMDestroy(&finedm);CHKERRQ(ierr); 113 ierr = DMDestroy(&coarsedm);CHKERRQ(ierr); 114 ierr = PetscFinalize(); 115 return ierr; 116 } 117 118 /*TEST 119 120 test: 121 args: -M 10 -N 10 -Z 10 122 output_file: output/ex89_1.out 123 124 test: 125 suffix: allatonce 126 nsize: 4 127 args: -M 10 -N 10 -Z 10 -matptap_via allatonce 128 output_file: output/ex89_1.out 129 130 test: 131 suffix: allatonce_merged 132 nsize: 4 133 args: -M 10 -M 5 -M 10 -matptap_via allatonce_merged 134 output_file: output/ex96_1.out 135 136 test: 137 suffix: allatonce_3D 138 nsize: 4 139 args: -M 10 -M 5 -M 10 -test_3D 1 -matptap_via allatonce 140 output_file: output/ex96_1.out 141 142 test: 143 suffix: allatonce_merged_3D 144 nsize: 4 145 args: -M 10 -M 5 -M 10 -test_3D 1 -matptap_via allatonce_merged 146 output_file: output/ex96_1.out 147 148 TEST*/ 149