Actual source code: ex31.c

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

  6: #include <petscksp.h>

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

148:   PetscCall(PetscFinalize());
149:   return 0;
150: }

152: /*TEST

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

160: TEST*/