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,none = -1.0,alpha; 15 Vec x,v1,v2,v3,v4; 16 PetscReal norm,norm_tmp,norm_tmp1,tol=100.*PETSC_MACHINE_EPSILON; 17 PetscRandom rdm; 18 PetscBool Test_3D=PETSC_FALSE,flg; 19 const PetscInt *ia,*ja; 20 PetscInt dof; 21 MPI_Comm comm; 22 23 ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr; 24 comm = PETSC_COMM_WORLD; 25 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 26 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 27 M = 10; N = 10; Z = 10; 28 dof = 10; 29 30 ierr = PetscOptionsGetBool(NULL,NULL,"-test_3D",&Test_3D,NULL);CHKERRQ(ierr); 31 ierr = PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);CHKERRQ(ierr); 32 ierr = PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);CHKERRQ(ierr); 33 ierr = PetscOptionsGetInt(NULL,NULL,"-Z",&Z,NULL);CHKERRQ(ierr); 34 /* Set up distributed array for fine grid */ 35 if (!Test_3D) { 36 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); 37 } else { 38 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); 39 } 40 ierr = DMSetFromOptions(coarsedm);CHKERRQ(ierr); 41 ierr = DMSetUp(coarsedm);CHKERRQ(ierr); 42 43 /* This makes sure the coarse DMDA has the same partition as the fine DMDA */ 44 ierr = DMRefine(coarsedm,PetscObjectComm((PetscObject)coarsedm),&finedm);CHKERRQ(ierr); 45 46 /*------------------------------------------------------------*/ 47 ierr = DMSetMatType(finedm,MATAIJ);CHKERRQ(ierr); 48 ierr = DMCreateMatrix(finedm,&A);CHKERRQ(ierr); 49 50 /* set val=one to A */ 51 if (size == 1) { 52 ierr = MatGetRowIJ(A,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);CHKERRQ(ierr); 53 if (flg) { 54 ierr = MatSeqAIJGetArray(A,&array);CHKERRQ(ierr); 55 for (i=0; i<ia[nrows]; i++) array[i] = one; 56 ierr = MatSeqAIJRestoreArray(A,&array);CHKERRQ(ierr); 57 } 58 ierr = MatRestoreRowIJ(A,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);CHKERRQ(ierr); 59 } else { 60 Mat AA,AB; 61 ierr = MatMPIAIJGetSeqAIJ(A,&AA,&AB,NULL);CHKERRQ(ierr); 62 ierr = MatGetRowIJ(AA,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);CHKERRQ(ierr); 63 if (flg) { 64 ierr = MatSeqAIJGetArray(AA,&array);CHKERRQ(ierr); 65 for (i=0; i<ia[nrows]; i++) array[i] = one; 66 ierr = MatSeqAIJRestoreArray(AA,&array);CHKERRQ(ierr); 67 } 68 ierr = MatRestoreRowIJ(AA,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);CHKERRQ(ierr); 69 ierr = MatGetRowIJ(AB,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);CHKERRQ(ierr); 70 if (flg) { 71 ierr = MatSeqAIJGetArray(AB,&array);CHKERRQ(ierr); 72 for (i=0; i<ia[nrows]; i++) array[i] = one; 73 ierr = MatSeqAIJRestoreArray(AB,&array);CHKERRQ(ierr); 74 } 75 ierr = MatRestoreRowIJ(AB,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);CHKERRQ(ierr); 76 } 77 /* Create interpolation between the fine and coarse grids */ 78 ierr = DMCreateInterpolation(coarsedm,finedm,&P,NULL);CHKERRQ(ierr); 79 /* Create vectors v1 and v2 that are compatible with A */ 80 ierr = MatCreateVecs(A,&v1,NULL);CHKERRQ(ierr); 81 ierr = VecDuplicate(v1,&v2);CHKERRQ(ierr); 82 ierr = PetscRandomCreate(comm,&rdm);CHKERRQ(ierr); 83 ierr = PetscRandomSetFromOptions(rdm);CHKERRQ(ierr); 84 85 /* Test P^T * A * P - MatPtAP() */ 86 /*------------------------------*/ 87 ierr = MatPtAP(A,P,MAT_INITIAL_MATRIX,fill,&C);CHKERRQ(ierr); 88 /* Test MAT_REUSE_MATRIX - reuse symbolic C */ 89 alpha=1.0; 90 for (i=0; i<1; i++) { 91 alpha -= 0.1; 92 ierr = MatScale(A,alpha);CHKERRQ(ierr); 93 ierr = MatPtAP(A,P,MAT_REUSE_MATRIX,fill,&C);CHKERRQ(ierr); 94 } 95 96 /* Free intermediate data structures created for reuse of C=Pt*A*P */ 97 ierr = MatFreeIntermediateDataStructures(C);CHKERRQ(ierr); 98 99 /* Create vector x that is compatible with P */ 100 ierr = MatCreateVecs(P,&x,NULL);CHKERRQ(ierr); 101 ierr = VecDuplicate(x,&v3);CHKERRQ(ierr); 102 ierr = VecDuplicate(x,&v4);CHKERRQ(ierr); 103 104 norm = 0.0; 105 for (i=0; i<10; i++) { 106 ierr = VecSetRandom(x,rdm);CHKERRQ(ierr); 107 ierr = MatMult(P,x,v1);CHKERRQ(ierr); 108 ierr = MatMult(A,v1,v2);CHKERRQ(ierr); /* v2 = A*P*x */ 109 110 ierr = MatMultTranspose(P,v2,v3);CHKERRQ(ierr); /* v3 = Pt*A*P*x */ 111 ierr = MatMult(C,x,v4);CHKERRQ(ierr); /* v3 = C*x */ 112 ierr = VecAXPY(v4,none,v3);CHKERRQ(ierr); 113 ierr = VecNorm(v4,NORM_1,&norm_tmp);CHKERRQ(ierr); 114 ierr = VecNorm(v3,NORM_1,&norm_tmp1);CHKERRQ(ierr); 115 116 norm_tmp /= norm_tmp1; 117 if (norm_tmp > norm) norm = norm_tmp; 118 } 119 if (norm >= tol && !rank) { 120 ierr = PetscPrintf(PETSC_COMM_SELF,"Error: MatPtAP(), |v3 - v4|/|v3|: %g\n",(double)norm);CHKERRQ(ierr); 121 } 122 123 ierr = MatDestroy(&C);CHKERRQ(ierr); 124 ierr = MatDestroy(&A);CHKERRQ(ierr); 125 ierr = MatDestroy(&P);CHKERRQ(ierr); 126 ierr = VecDestroy(&v3);CHKERRQ(ierr); 127 ierr = VecDestroy(&v4);CHKERRQ(ierr); 128 ierr = VecDestroy(&x);CHKERRQ(ierr); 129 ierr = VecDestroy(&v1);CHKERRQ(ierr); 130 ierr = VecDestroy(&v2);CHKERRQ(ierr); 131 ierr = DMDestroy(&finedm);CHKERRQ(ierr); 132 ierr = DMDestroy(&coarsedm);CHKERRQ(ierr); 133 ierr = PetscRandomDestroy(&rdm);CHKERRQ(ierr); 134 ierr = PetscFinalize(); 135 return ierr; 136 } 137 138 /*TEST 139 140 test: 141 args: -M 10 -N 10 -Z 10 142 output_file: output/ex89_1.out 143 144 test: 145 suffix: allatonce 146 nsize: 4 147 args: -M 10 -N 10 -Z 10 -matmaijptap_via allatonce 148 output_file: output/ex89_1.out 149 150 test: 151 suffix: allatonce_merged 152 nsize: 4 153 args: -M 10 -M 5 -M 10 -matmaijptap_via allatonce_merged 154 output_file: output/ex96_1.out 155 156 test: 157 suffix: allatonce_3D 158 nsize: 4 159 args: -M 10 -M 5 -M 10 -test_3D 1 -matmaijptap_via allatonce 160 output_file: output/ex96_1.out 161 162 test: 163 suffix: allatonce_merged_3D 164 nsize: 4 165 args: -M 10 -M 5 -M 10 -test_3D 1 -matmaijptap_via allatonce_merged 166 output_file: output/ex96_1.out 167 168 TEST*/ 169