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

  7: !
  8: !  This example is the Fortran version of ex6.c.  The program reads a PETSc matrix
  9: !  and vector from a file and solves a linear system.  Input arguments are:
 10: !        -f <input_file> : file to load.  For example see $PETSC_DIR/share/petsc/datafiles/matrices
 11: !

 13:   PetscErrorCode ierr
 14:   PetscInt its, m, n, mlocal, nlocal
 15:   PetscBool flg
 16:   PetscScalar, parameter :: none = -1.0
 17:   PetscReal norm
 18:   Vec x, b, u
 19:   Mat A
 20:   character(len=128) f
 21:   PetscViewer fd
 22:   MatInfo info
 23:   KSP ksp

 25:   PetscCallA(PetscInitialize(ierr))

 27: ! Read in matrix and RHS
 28:   PetscCallA(PetscOptionsGetString(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-f', f, flg, ierr))
 29:   PetscCallA(PetscViewerBinaryOpen(PETSC_COMM_WORLD, f, FILE_MODE_READ, fd, ierr))

 31:   PetscCallA(MatCreate(PETSC_COMM_WORLD, A, ierr))
 32:   PetscCallA(MatSetType(A, MATSEQAIJ, ierr))
 33:   PetscCallA(MatLoad(A, fd, ierr))

 35: ! Get information about matrix
 36:   PetscCallA(MatGetSize(A, m, n, ierr))
 37:   PetscCallA(MatGetLocalSize(A, mlocal, nlocal, ierr))
 38:   PetscCallA(MatGetInfo(A, MAT_GLOBAL_SUM, info, ierr))
 39:   write (*, 100) m, &
 40:     n, &
 41:     mlocal, nlocal, &
 42:     info%BLOCK_SIZE, info%NZ_ALLOCATED, &
 43:     info%NZ_USED, info%NZ_UNNEEDED, &
 44:     info%MEMORY, info%ASSEMBLIES, &
 45:     info%MALLOCS

 47: 100 format(4(i4, 1x), 7(1pe9.2, 1x))
 48:   PetscCallA(VecCreate(PETSC_COMM_WORLD, b, ierr))
 49:   PetscCallA(VecLoad(b, fd, ierr))
 50:   PetscCallA(PetscViewerDestroy(fd, ierr))

 52: ! Set up solution
 53:   PetscCallA(VecDuplicate(b, x, ierr))
 54:   PetscCallA(VecDuplicate(b, u, ierr))

 56: ! Solve system
 57:   PetscCallA(KSPCreate(PETSC_COMM_WORLD, ksp, ierr))
 58:   PetscCallA(KSPSetOperators(ksp, A, A, ierr))
 59:   PetscCallA(KSPSetFromOptions(ksp, ierr))
 60:   PetscCallA(KSPSolve(ksp, b, x, ierr))

 62: ! Show result
 63:   PetscCallA(MatMult(A, x, u, ierr))
 64:   PetscCallA(VecAXPY(u, none, b, ierr))
 65:   PetscCallA(VecNorm(u, NORM_2, norm, ierr))
 66:   PetscCallA(KSPGetIterationNumber(ksp, its, ierr))
 67:   write (6, 101) norm, its
 68: 101 format('Residual norm ', 1pe9.2, ' iterations ', i5)

 70: ! Cleanup
 71:   PetscCallA(KSPDestroy(ksp, ierr))
 72:   PetscCallA(VecDestroy(b, ierr))
 73:   PetscCallA(VecDestroy(x, ierr))
 74:   PetscCallA(VecDestroy(u, ierr))
 75:   PetscCallA(MatDestroy(A, ierr))

 77:   PetscCallA(PetscFinalize(ierr))
 78: end

 80: !/*TEST
 81: !
 82: !    test:
 83: !      args: -f ${DATAFILESPATH}/matrices/arco1 -options_left no
 84: !      requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
 85: !
 86: !TEST*/