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:   PetscFunctionBeginUser;
 23:   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
 24:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 25:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));

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

 31:   /* Determine file from which we read the matrix.*/
 32:   PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), &flg));
 33:   PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT, "Must indicate binary file with the -f option");

 35:   /* - - - - - - - - - - - - - - - - - - - - - - - -
 36:                            Load system
 37:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 38:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd));
 39:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 40:   PetscCall(MatLoad(A, fd));
 41:   PetscCall(PetscViewerDestroy(&fd));
 42:   PetscCall(MatGetLocalSize(A, &m, &n));
 43:   PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")", m, n);

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

 51:   PetscCall(VecDuplicate(b, &x));
 52:   PetscCall(VecDuplicate(b, &u));
 53:   PetscCall(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:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Before partitioning/reordering, A:\n"));
 66:       PetscCall(MatView(A, PETSC_VIEWER_DRAW_WORLD));
 67:     }

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

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

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

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

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

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

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

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

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

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

149:   PetscCall(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*/