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;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

149:   PetscFinalize();
150:   return 0;
151: }

153: /*TEST

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

161: TEST*/