Actual source code: ex31.c

petsc-3.12.2 2019-11-22
Report Typos and Errors

  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: /*T
  8:    Concepts: KSP^solving a linear system
  9:    Processors: n
 10: T*/


 13:  #include <petscksp.h>

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

 29:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
 30:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 31:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

 33:   PetscOptionsGetBool(NULL,NULL,"-partition",&partition,NULL);
 34:   PetscOptionsGetBool(NULL,NULL,"-displayIS",&displayIS,NULL);
 35:   PetscOptionsGetBool(NULL,NULL,"-displayMat",&displayMat,NULL);

 37:   /* Determine file from which we read the matrix.*/
 38:   PetscOptionsGetString(NULL,NULL,"-f",file,PETSC_MAX_PATH_LEN,&flg);
 39:   if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate binary file with the -f option");

 41:   /* - - - - - - - - - - - - - - - - - - - - - - - -
 42:                            Load system
 43:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 44:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
 45:   MatCreate(PETSC_COMM_WORLD,&A);
 46:   MatLoad(A,fd);
 47:   PetscViewerDestroy(&fd);
 48:   MatGetLocalSize(A,&m,&n);
 49:   if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%D, %D)", m, n);

 51:   /* Create rhs vector of all ones */
 52:   VecCreate(PETSC_COMM_WORLD,&b);
 53:   VecSetSizes(b,m,PETSC_DECIDE);
 54:   VecSetFromOptions(b);
 55:   VecSet(b,one);

 57:   VecDuplicate(b,&x);
 58:   VecDuplicate(b,&u);
 59:   VecSet(x,0.0);

 61:   /* - - - - - - - - - - - - - - - - - - - - - - - -
 62:                       Test partition
 63:   - - - - - - - - - - - - - - - - - - - - - - - - - */
 64:   if (partition) {
 65:     MatPartitioning mpart;
 66:     IS              mis,nis,is;
 67:     PetscInt        *count;
 68:     Mat             BB;

 70:     if (displayMat) {
 71:       PetscPrintf(PETSC_COMM_WORLD,"Before partitioning/reordering, A:\n");
 72:       MatView(A,PETSC_VIEWER_DRAW_WORLD);
 73:     }

 75:     PetscMalloc1(size,&count);
 76:     MatPartitioningCreate(PETSC_COMM_WORLD, &mpart);
 77:     MatPartitioningSetAdjacency(mpart, A);
 78:     /* MatPartitioningSetVertexWeights(mpart, weight); */
 79:     MatPartitioningSetFromOptions(mpart);
 80:     MatPartitioningApply(mpart, &mis);
 81:     MatPartitioningDestroy(&mpart);
 82:     if (displayIS) {
 83:       PetscPrintf(PETSC_COMM_WORLD,"mis, new processor assignment:\n");
 84:       ISView(mis,PETSC_VIEWER_STDOUT_WORLD);
 85:     }

 87:     ISPartitioningToNumbering(mis,&nis);
 88:     if (displayIS) {
 89:       PetscPrintf(PETSC_COMM_WORLD,"nis:\n");
 90:       ISView(nis,PETSC_VIEWER_STDOUT_WORLD);
 91:     }

 93:     ISPartitioningCount(mis,size,count);
 94:     ISDestroy(&mis);
 95:     if (displayIS && !rank) {
 96:       PetscInt i;
 97:       PetscPrintf(PETSC_COMM_SELF,"[ %d ] count:\n",rank);
 98:       for (i=0; i<size; i++) {PetscPrintf(PETSC_COMM_WORLD," %d",count[i]);}
 99:       PetscPrintf(PETSC_COMM_WORLD,"\n");
100:     }

102:     ISInvertPermutation(nis, count[rank], &is);
103:     PetscFree(count);
104:     ISDestroy(&nis);
105:     ISSort(is);
106:     if (displayIS) {
107:       PetscPrintf(PETSC_COMM_WORLD,"inverse of nis - maps new local rows to old global rows:\n");
108:       ISView(is,PETSC_VIEWER_STDOUT_WORLD);
109:     }

111:     MatCreateSubMatrix(A,is,is,MAT_INITIAL_MATRIX,&BB);
112:     if (displayMat) {
113:       PetscPrintf(PETSC_COMM_WORLD,"After partitioning/reordering, A:\n");
114:       MatView(BB,PETSC_VIEWER_DRAW_WORLD);
115:     }

117:     /* need to move the vector also */
118:     ISDestroy(&is);
119:     MatDestroy(&A);
120:     A    = BB;
121:   }

123:   /* Create linear solver; set operators; set runtime options.*/
124:   KSPCreate(PETSC_COMM_WORLD,&ksp);
125:   KSPSetOperators(ksp,A,A);
126:   KSPSetFromOptions(ksp);

128:   /* - - - - - - - - - - - - - - - - - - - - - - - -
129:                            Solve system
130:         - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
131:   KSPSolve(ksp,b,x);
132:   KSPGetIterationNumber(ksp,&its);

134:   /* Check error */
135:   MatMult(A,x,u);
136:   VecAXPY(u,-1.0,b);
137:   VecNorm(u,NORM_2,&norm);
138:   PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
139:   PetscPrintf(PETSC_COMM_WORLD,"Residual norm %g\n",(double)norm);
140:   flg  = PETSC_FALSE;
141:   PetscOptionsGetBool(NULL,NULL, "-ksp_reason", &flg,NULL);
142:   if (flg) {
143:     KSPConvergedReason reason;
144:     KSPGetConvergedReason(ksp,&reason);
145:     PetscPrintf(PETSC_COMM_WORLD,"KSPConvergedReason: %D\n", reason);
146:   }

148:   /* Free work space.*/
149:   MatDestroy(&A); VecDestroy(&b);
150:   VecDestroy(&u); VecDestroy(&x);
151:   KSPDestroy(&ksp);

153:   PetscFinalize();
154:   return ierr;
155: }

157: /*TEST

159:     test:
160:       args: -f ${DATAFILESPATH}/matrices/small -partition -mat_partitioning_type parmetis 
161:       requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES) parmetis
162:       output_file: output/ex31.out
163:       nsize: 3
164:  
165: TEST*/