Actual source code: ex27.c


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

138: /*TEST

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

145: TEST*/