Actual source code: ex31.c


  2: static char help[] = "Test partition. Reads a PETSc matrix and vector from a file and solves a linear system.\n\
  3: This   Input parameters include\n\
  4:   -f <input_file> : file to load \n\
  5:   -partition -mat_partitioning_view \n\\n";

  7: #include <petscksp.h>

  9: int main(int argc,char **args)
 10: {
 11:   KSP            ksp;             /* linear solver context */
 12:   Mat            A;               /* matrix */
 13:   Vec            x,b,u;           /* approx solution, RHS, exact solution */
 14:   PetscViewer    fd;              /* viewer */
 15:   char           file[PETSC_MAX_PATH_LEN];     /* input file name */
 16:   PetscBool      flg,partition=PETSC_FALSE,displayIS=PETSC_FALSE,displayMat=PETSC_FALSE;
 17:   PetscInt       its,m,n;
 18:   PetscReal      norm;
 19:   PetscMPIInt    size,rank;
 20:   PetscScalar    one = 1.0;

 22:   PetscInitialize(&argc,&args,(char*)0,help);
 23:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 24:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

 26:   PetscOptionsGetBool(NULL,NULL,"-partition",&partition,NULL);
 27:   PetscOptionsGetBool(NULL,NULL,"-displayIS",&displayIS,NULL);
 28:   PetscOptionsGetBool(NULL,NULL,"-displayMat",&displayMat,NULL);

 30:   /* Determine file from which we read the matrix.*/
 31:   PetscOptionsGetString(NULL,NULL,"-f",file,sizeof(file),&flg);

 34:   /* - - - - - - - - - - - - - - - - - - - - - - - -
 35:                            Load system
 36:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 37:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
 38:   MatCreate(PETSC_COMM_WORLD,&A);
 39:   MatLoad(A,fd);
 40:   PetscViewerDestroy(&fd);
 41:   MatGetLocalSize(A,&m,&n);

 44:   /* Create rhs vector of all ones */
 45:   VecCreate(PETSC_COMM_WORLD,&b);
 46:   VecSetSizes(b,m,PETSC_DECIDE);
 47:   VecSetFromOptions(b);
 48:   VecSet(b,one);

 50:   VecDuplicate(b,&x);
 51:   VecDuplicate(b,&u);
 52:   VecSet(x,0.0);

 54:   /* - - - - - - - - - - - - - - - - - - - - - - - -
 55:                       Test partition
 56:   - - - - - - - - - - - - - - - - - - - - - - - - - */
 57:   if (partition) {
 58:     MatPartitioning mpart;
 59:     IS              mis,nis,is;
 60:     PetscInt        *count;
 61:     Mat             BB;

 63:     if (displayMat) {
 64:       PetscPrintf(PETSC_COMM_WORLD,"Before partitioning/reordering, A:\n");
 65:       MatView(A,PETSC_VIEWER_DRAW_WORLD);
 66:     }

 68:     PetscMalloc1(size,&count);
 69:     MatPartitioningCreate(PETSC_COMM_WORLD, &mpart);
 70:     MatPartitioningSetAdjacency(mpart, A);
 71:     /* MatPartitioningSetVertexWeights(mpart, weight); */
 72:     MatPartitioningSetFromOptions(mpart);
 73:     MatPartitioningApply(mpart, &mis);
 74:     MatPartitioningDestroy(&mpart);
 75:     if (displayIS) {
 76:       PetscPrintf(PETSC_COMM_WORLD,"mis, new processor assignment:\n");
 77:       ISView(mis,PETSC_VIEWER_STDOUT_WORLD);
 78:     }

 80:     ISPartitioningToNumbering(mis,&nis);
 81:     if (displayIS) {
 82:       PetscPrintf(PETSC_COMM_WORLD,"nis:\n");
 83:       ISView(nis,PETSC_VIEWER_STDOUT_WORLD);
 84:     }

 86:     ISPartitioningCount(mis,size,count);
 87:     ISDestroy(&mis);
 88:     if (displayIS && rank == 0) {
 89:       PetscInt i;
 90:       PetscPrintf(PETSC_COMM_SELF,"[ %d ] count:\n",rank);
 91:       for (i=0; i<size; i++) PetscPrintf(PETSC_COMM_WORLD," %d",count[i]);
 92:       PetscPrintf(PETSC_COMM_WORLD,"\n");
 93:     }

 95:     ISInvertPermutation(nis, count[rank], &is);
 96:     PetscFree(count);
 97:     ISDestroy(&nis);
 98:     ISSort(is);
 99:     if (displayIS) {
100:       PetscPrintf(PETSC_COMM_WORLD,"inverse of nis - maps new local rows to old global rows:\n");
101:       ISView(is,PETSC_VIEWER_STDOUT_WORLD);
102:     }

104:     MatCreateSubMatrix(A,is,is,MAT_INITIAL_MATRIX,&BB);
105:     if (displayMat) {
106:       PetscPrintf(PETSC_COMM_WORLD,"After partitioning/reordering, A:\n");
107:       MatView(BB,PETSC_VIEWER_DRAW_WORLD);
108:     }

110:     /* need to move the vector also */
111:     ISDestroy(&is);
112:     MatDestroy(&A);
113:     A    = BB;
114:   }

116:   /* Create linear solver; set operators; set runtime options.*/
117:   KSPCreate(PETSC_COMM_WORLD,&ksp);
118:   KSPSetOperators(ksp,A,A);
119:   KSPSetFromOptions(ksp);

121:   /* - - - - - - - - - - - - - - - - - - - - - - - -
122:                            Solve system
123:         - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
124:   KSPSolve(ksp,b,x);
125:   KSPGetIterationNumber(ksp,&its);

127:   /* Check error */
128:   MatMult(A,x,u);
129:   VecAXPY(u,-1.0,b);
130:   VecNorm(u,NORM_2,&norm);
131:   PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
132:   PetscPrintf(PETSC_COMM_WORLD,"Residual norm %g\n",(double)norm);
133:   flg  = PETSC_FALSE;
134:   PetscOptionsGetBool(NULL,NULL, "-ksp_reason", &flg,NULL);
135:   if (flg) {
136:     KSPConvergedReason reason;
137:     KSPGetConvergedReason(ksp,&reason);
138:     PetscPrintf(PETSC_COMM_WORLD,"KSPConvergedReason: %D\n", reason);
139:   }

141:   /* Free work space.*/
142:   MatDestroy(&A)); PetscCall(VecDestroy(&b);
143:   VecDestroy(&u)); PetscCall(VecDestroy(&x);
144:   KSPDestroy(&ksp);

146:   PetscFinalize();
147:   return 0;
148: }

150: /*TEST

152:     test:
153:       args: -f ${DATAFILESPATH}/matrices/small -partition -mat_partitioning_type parmetis
154:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) parmetis
155:       output_file: output/ex31.out
156:       nsize: 3

158: TEST*/