Actual source code: ex27.c

  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, (char *)0, 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));
 74:   PetscCall(VecSet(x, 0.0));

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

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

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

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

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

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

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

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

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

137: /*TEST

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

144: TEST*/