Actual source code: ex89.c

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

  3: #include <petscdmda.h>

  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;

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

 26:   PetscOptionsGetBool(NULL,NULL,"-test_3D",&Test_3D,NULL);
 27:   PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);
 28:   PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);
 29:   PetscOptionsGetInt(NULL,NULL,"-Z",&Z,NULL);
 30:   /* Set up distributed array for fine grid */
 31:   if (!Test_3D) {
 32:     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:     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:   DMSetFromOptions(coarsedm);
 37:   DMSetUp(coarsedm);

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

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

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

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

 88:   { /* Test MatProductView() */
 89:     PetscViewer viewer;
 90:     PetscViewerASCIIOpen(comm,NULL, &viewer);
 91:     PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
 92:     MatProductView(C,viewer);
 93:     PetscViewerPopFormat(viewer);
 94:     PetscViewerDestroy(&viewer);
 95:   }

 97:   MatPtAPMultEqual(A,P,C,10,&flg);
 99:   MatDestroy(&C);

101:   /* (2) User API */
102:   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:     MatScale(A,alpha);
108:     MatPtAP(A,P,MAT_REUSE_MATRIX,fill,&C);
109:   }

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

114:   MatPtAPMultEqual(A,P,C,10,&flg);

117:   MatDestroy(&C);
118:   MatDestroy(&A);
119:   MatDestroy(&P);
120:   DMDestroy(&finedm);
121:   DMDestroy(&coarsedm);
122:   PetscFinalize();
123:   return 0;
124: }

126: /*TEST

128:    test:
129:       args: -M 10 -N 10 -Z 10
130:       output_file: output/ex89_1.out

132:    test:
133:       suffix: allatonce
134:       nsize: 4
135:       args: -M 10 -N 10 -Z 10
136:       output_file: output/ex89_2.out

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

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

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

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

162: TEST*/