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