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

  5: /*
  6:   Usage:
  7:      ex27 -f0 <mat_binaryfile>
  8: */

 10: #include <petscksp.h>
 11: extern PetscErrorCode PCShellApply_Matinv(PC, Vec, Vec);

 13: int main(int argc, char **args)
 14: {
 15:   KSP         ksp;
 16:   Mat         A, B, F, X;
 17:   Vec         x, b, u;                     /* approx solution, RHS, exact solution */
 18:   PetscViewer fd;                          /* viewer */
 19:   char        file[1][PETSC_MAX_PATH_LEN]; /* input file name */
 20:   PetscBool   flg;
 21:   PetscInt    M, N, i, its;
 22:   PetscReal   norm;
 23:   PetscScalar val = 1.0;
 24:   PetscMPIInt size;
 25:   PC          pc;

 27:   PetscFunctionBeginUser;
 28:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 29:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 30:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");

 32:   /* Read matrix and right-hand-side vector */
 33:   PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file[0], sizeof(file[0]), &flg));
 34:   PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate binary file with the -f option");

 36:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file[0], FILE_MODE_READ, &fd));
 37:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 38:   PetscCall(MatSetType(A, MATAIJ));
 39:   PetscCall(MatLoad(A, fd));
 40:   PetscCall(VecCreate(PETSC_COMM_WORLD, &b));
 41:   PetscCall(VecLoad(b, fd));
 42:   PetscCall(PetscViewerDestroy(&fd));

 44:   /*
 45:      If the loaded matrix is larger than the vector (due to being padded
 46:      to match the block size of the system), then create a new padded vector.
 47:   */
 48:   {
 49:     PetscInt     m, n, j, mvec, start, end, indx;
 50:     Vec          tmp;
 51:     PetscScalar *bold;

 53:     /* Create a new vector b by padding the old one */
 54:     PetscCall(MatGetLocalSize(A, &m, &n));
 55:     PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")", m, n);
 56:     PetscCall(VecCreate(PETSC_COMM_WORLD, &tmp));
 57:     PetscCall(VecSetSizes(tmp, m, PETSC_DECIDE));
 58:     PetscCall(VecSetFromOptions(tmp));
 59:     PetscCall(VecGetOwnershipRange(b, &start, &end));
 60:     PetscCall(VecGetLocalSize(b, &mvec));
 61:     PetscCall(VecGetArray(b, &bold));
 62:     for (j = 0; j < mvec; j++) {
 63:       indx = start + j;
 64:       PetscCall(VecSetValues(tmp, 1, &indx, bold + j, INSERT_VALUES));
 65:     }
 66:     PetscCall(VecRestoreArray(b, &bold));
 67:     PetscCall(VecDestroy(&b));
 68:     PetscCall(VecAssemblyBegin(tmp));
 69:     PetscCall(VecAssemblyEnd(tmp));
 70:     b = tmp;
 71:   }
 72:   PetscCall(VecDuplicate(b, &x));
 73:   PetscCall(VecDuplicate(b, &u));

 75:   /* Create dense matrices B and X. Set B as an identity matrix */
 76:   PetscCall(MatGetSize(A, &M, &N));
 77:   PetscCall(MatCreate(MPI_COMM_SELF, &B));
 78:   PetscCall(MatSetSizes(B, M, N, M, N));
 79:   PetscCall(MatSetType(B, MATSEQDENSE));
 80:   PetscCall(MatSeqDenseSetPreallocation(B, NULL));
 81:   for (i = 0; i < M; i++) PetscCall(MatSetValues(B, 1, &i, 1, &i, &val, INSERT_VALUES));
 82:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
 83:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));

 85:   PetscCall(MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &X));

 87:   /* Compute X=inv(A) by MatMatSolve() */
 88:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
 89:   PetscCall(KSPSetOperators(ksp, A, A));
 90:   PetscCall(KSPGetPC(ksp, &pc));
 91:   PetscCall(PCSetType(pc, PCLU));
 92:   PetscCall(KSPSetFromOptions(ksp));
 93:   PetscCall(KSPSetUp(ksp));
 94:   PetscCall(PCFactorGetMatrix(pc, &F));
 95:   PetscCall(MatMatSolve(F, B, X));
 96:   PetscCall(MatDestroy(&B));

 98:   /* Now, set X=inv(A) as a preconditioner */
 99:   PetscCall(PCSetType(pc, PCSHELL));
100:   PetscCall(PCShellSetContext(pc, X));
101:   PetscCall(PCShellSetApply(pc, PCShellApply_Matinv));
102:   PetscCall(KSPSetFromOptions(ksp));

104:   /* Solve preconditioned system A*x = b */
105:   PetscCall(KSPSolve(ksp, b, x));
106:   PetscCall(KSPGetIterationNumber(ksp, &its));

108:   /* Check error */
109:   PetscCall(MatMult(A, x, u));
110:   PetscCall(VecAXPY(u, -1.0, b));
111:   PetscCall(VecNorm(u, NORM_2, &norm));
112:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of iterations = %3" PetscInt_FMT "\n", its));
113:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Residual norm %g\n", (double)norm));

115:   /* Free work space.  */
116:   PetscCall(MatDestroy(&X));
117:   PetscCall(MatDestroy(&A));
118:   PetscCall(VecDestroy(&b));
119:   PetscCall(VecDestroy(&u));
120:   PetscCall(VecDestroy(&x));
121:   PetscCall(KSPDestroy(&ksp));
122:   PetscCall(PetscFinalize());
123:   return 0;
124: }

126: PetscErrorCode PCShellApply_Matinv(PC pc, Vec xin, Vec xout)
127: {
128:   Mat X;

130:   PetscFunctionBeginUser;
131:   PetscCall(PCShellGetContext(pc, &X));
132:   PetscCall(MatMult(X, xin, xout));
133:   PetscFunctionReturn(PETSC_SUCCESS);
134: }

136: /*TEST

138:     test:
139:       args: -f ${DATAFILESPATH}/matrices/small
140:       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
141:       output_file: output/ex27.out

143: TEST*/