1: !
  2: #include <petsc/finclude/petscksp.h>
  3: program main
  4:   use petscksp
  5:   implicit none

  7: !
  8: !  This example is a modified Fortran version of ex6.c.  It tests the use of
  9: !  options prefixes in PETSc. Two linear problems are solved in this program.
 10: !  The first problem is read from a file. The second problem is constructed
 11: !  from the first, by eliminating some of the entries of the linear matrix 'A'.

 13: !  Each solve is distinguished by a unique prefix - 'a' for the first, 'b'
 14: !  for the second.  With the prefix the user can distinguish between the various
 15: !  options (command line, from .petscrc file, etc.) for each of the solvers.
 16: !  Input arguments are:
 17: !     -f <input_file> : file to load.  For a 5X5 example of the 5-pt. stencil
 18: !                       use the file petsc/src/mat/examples/mat.ex.binary

 20:   PetscErrorCode ierr
 21:   PetscInt its
 22:   PetscBool flg
 23:   PetscScalar, parameter :: none = -1.0, five = 5.0
 24:   PetscReal norm
 25:   Vec x, b, u
 26:   Mat A
 27:   KSP ksp1, ksp2
 28:   character(len=PETSC_MAX_PATH_LEN) f
 29:   PetscViewer fd
 30:   IS isrow

 32:   PetscCallA(PetscInitialize(ierr))

 34: !     Read in matrix and RHS
 35:   PetscCallA(PetscOptionsGetString(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-f', f, flg, ierr))
 36:   PetscCallA(PetscViewerBinaryOpen(PETSC_COMM_WORLD, f, FILE_MODE_READ, fd, ierr))

 38:   PetscCallA(MatCreate(PETSC_COMM_WORLD, A, ierr))
 39:   PetscCallA(MatSetType(A, MATSEQAIJ, ierr))
 40:   PetscCallA(MatLoad(A, fd, ierr))

 42:   PetscCallA(VecCreate(PETSC_COMM_WORLD, b, ierr))
 43:   PetscCallA(VecLoad(b, fd, ierr))
 44:   PetscCallA(PetscViewerDestroy(fd, ierr))

 46: ! Set up solution
 47:   PetscCallA(VecDuplicate(b, x, ierr))
 48:   PetscCallA(VecDuplicate(b, u, ierr))

 50: ! Solve system-1
 51:   PetscCallA(KSPCreate(PETSC_COMM_WORLD, ksp1, ierr))
 52:   PetscCallA(KSPSetOptionsPrefix(ksp1, 'a', ierr))
 53:   PetscCallA(KSPAppendOptionsPrefix(ksp1, '_', ierr))
 54:   PetscCallA(KSPSetOperators(ksp1, A, A, ierr))
 55:   PetscCallA(KSPSetFromOptions(ksp1, ierr))
 56:   PetscCallA(KSPSolve(ksp1, b, x, ierr))

 58: ! Show result
 59:   PetscCallA(MatMult(A, x, u, ierr))
 60:   PetscCallA(VecAXPY(u, none, b, ierr))
 61:   PetscCallA(VecNorm(u, NORM_2, norm, ierr))
 62:   PetscCallA(KSPGetIterationNumber(ksp1, its, ierr))

 64:   write (6, 100) norm, its
 65: 100 format('Residual norm ', e11.4, ' iterations ', i5)

 67: ! Create system 2 by striping off some rows of the matrix
 68:   PetscCallA(ISCreateStride(PETSC_COMM_SELF, 5_PETSC_INT_KIND, 0_PETSC_INT_KIND, 1_PETSC_INT_KIND, isrow, ierr))
 69:   PetscCallA(MatZeroRowsIS(A, isrow, five, PETSC_NULL_VEC, PETSC_NULL_VEC, ierr))

 71: ! Solve system-2
 72:   PetscCallA(KSPCreate(PETSC_COMM_WORLD, ksp2, ierr))
 73:   PetscCallA(KSPSetOptionsPrefix(ksp2, 'b', ierr))
 74:   PetscCallA(KSPAppendOptionsPrefix(ksp2, '_', ierr))
 75:   PetscCallA(KSPSetOperators(ksp2, A, A, ierr))
 76:   PetscCallA(KSPSetFromOptions(ksp2, ierr))
 77:   PetscCallA(KSPSolve(ksp2, b, x, ierr))

 79: ! Show result
 80:   PetscCallA(MatMult(A, x, u, ierr))
 81:   PetscCallA(VecAXPY(u, none, b, ierr))
 82:   PetscCallA(VecNorm(u, NORM_2, norm, ierr))
 83:   PetscCallA(KSPGetIterationNumber(ksp2, its, ierr))
 84:   write (6, 100) norm, its

 86: ! Cleanup
 87:   PetscCallA(KSPDestroy(ksp1, ierr))
 88:   PetscCallA(KSPDestroy(ksp2, ierr))
 89:   PetscCallA(VecDestroy(b, ierr))
 90:   PetscCallA(VecDestroy(x, ierr))
 91:   PetscCallA(VecDestroy(u, ierr))
 92:   PetscCallA(MatDestroy(A, ierr))
 93:   PetscCallA(ISDestroy(isrow, ierr))

 95:   PetscCallA(PetscFinalize(ierr))
 96: end

 98: !/*TEST
 99: !
100: !    test:
101: !      args: -f ${DATAFILESPATH}/matrices/arco1 -options_left no
102: !      requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
103: !
104: !TEST*/