1: static char help[] ="Tests MatPtAP() for MPIMAIJ and MPIAIJ \n ";

  3: #include <petscdmda.h>

  5: int main(int argc,char **argv)
  6: {
  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;

 20:   PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;
 21:   comm = PETSC_COMM_WORLD;
 22:   MPI_Comm_rank(comm,&rank);
 23:   MPI_Comm_size(comm,&size);
 24:   M = 10; N = 10; Z = 10;
 25:   dof  = 10;

 27:   PetscOptionsGetBool(NULL,NULL,"-test_3D",&Test_3D,NULL);
 28:   PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);
 29:   PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);
 30:   PetscOptionsGetInt(NULL,NULL,"-Z",&Z,NULL);
 31:   /* Set up distributed array for fine grid */
 32:   if (!Test_3D) {
 33:     DMDACreate2d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,M,N,PETSC_DECIDE,PETSC_DECIDE,dof,1,NULL,NULL,&coarsedm);
 34:   } else {
 35:     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);
 36:   }
 37:   DMSetFromOptions(coarsedm);
 38:   DMSetUp(coarsedm);

 40:   /* This makes sure the coarse DMDA has the same partition as the fine DMDA */
 41:   DMRefine(coarsedm,PetscObjectComm((PetscObject)coarsedm),&finedm);

 43:   /*------------------------------------------------------------*/
 44:   DMSetMatType(finedm,MATAIJ);
 45:   DMCreateMatrix(finedm,&A);

 47:   /* set val=one to A */
 48:   if (size == 1) {
 49:     MatGetRowIJ(A,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);
 50:     if (flg) {
 51:       MatSeqAIJGetArray(A,&array);
 52:       for (i=0; i<ia[nrows]; i++) array[i] = one;
 53:       MatSeqAIJRestoreArray(A,&array);
 54:     }
 55:     MatRestoreRowIJ(A,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);
 56:   } else {
 57:     Mat AA,AB;
 58:     MatMPIAIJGetSeqAIJ(A,&AA,&AB,NULL);
 59:     MatGetRowIJ(AA,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);
 60:     if (flg) {
 61:       MatSeqAIJGetArray(AA,&array);
 62:       for (i=0; i<ia[nrows]; i++) array[i] = one;
 63:       MatSeqAIJRestoreArray(AA,&array);
 64:     }
 65:     MatRestoreRowIJ(AA,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);
 66:     MatGetRowIJ(AB,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);
 67:     if (flg) {
 68:       MatSeqAIJGetArray(AB,&array);
 69:       for (i=0; i<ia[nrows]; i++) array[i] = one;
 70:       MatSeqAIJRestoreArray(AB,&array);
 71:     }
 72:     MatRestoreRowIJ(AB,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);
 73:   }
 74:   /* Create interpolation between the fine and coarse grids */
 75:   DMCreateInterpolation(coarsedm,finedm,&P,NULL);

 77:   /* Test P^T * A * P - MatPtAP() */
 78:   /*------------------------------*/
 79:   /* (1) Developer API */
 80:   MatProductCreate(A,P,NULL,&C);
 81:   MatProductSetType(C,MATPRODUCT_PtAP);
 82:   MatProductSetAlgorithm(C,"default");
 83:   MatProductSetFill(C,PETSC_DEFAULT);
 84:   MatProductSetFromOptions(C);
 85:   MatProductSymbolic(C);
 86:   MatProductNumeric(C);
 87:   MatProductNumeric(C); /* Test reuse of symbolic C */

 89:   MatPtAPMultEqual(A,P,C,10,&flg);
 90:   if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatProduct_PtAP");
 91:   MatDestroy(&C);

 93:   /* (2) User API */
 94:   MatPtAP(A,P,MAT_INITIAL_MATRIX,fill,&C);
 95:   /* Test MAT_REUSE_MATRIX - reuse symbolic C */
 96:   alpha=1.0;
 97:   for (i=0; i<1; i++) {
 98:     alpha -= 0.1;
 99:     MatScale(A,alpha);
100:     MatPtAP(A,P,MAT_REUSE_MATRIX,fill,&C);
101:   }

103:   /* Free intermediate data structures created for reuse of C=Pt*A*P */
104:   MatProductClear(C);

106:   MatPtAPMultEqual(A,P,C,10,&flg);
107:   if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatPtAP");

109:   MatDestroy(&C);
110:   MatDestroy(&A);
111:   MatDestroy(&P);
112:   DMDestroy(&finedm);
113:   DMDestroy(&coarsedm);
114:   PetscFinalize();
115:   return ierr;
116: }

118: /*TEST

120:    test:
121:       args: -M 10 -N 10 -Z 10
122:       output_file: output/ex89_1.out

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

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

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

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

148: TEST*/