Actual source code: ex12f.F90

  1: !
  2:       program main
  3: #include <petsc/finclude/petscksp.h>
  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      none
 17:       PetscReal        norm
 18:       Vec              x,b,u
 19:       Mat              A
 20:       character*(128)  f
 21:       PetscViewer      fd
 22:       MatInfo          info(MAT_INFO_SIZE)
 23:       KSP              ksp

 25:       none = -1.0
 26:       PetscCallA(PetscInitialize(ierr))

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

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

 36: ! Get information about matrix
 37:       PetscCallA(MatGetSize(A,m,n,ierr))
 38:       PetscCallA(MatGetLocalSize(A,mlocal,nlocal,ierr))
 39:       PetscCallA(MatGetInfo(A,MAT_GLOBAL_SUM,info,ierr))
 40:       write(*,100) m,                                                   &
 41:      &  n,                                                              &
 42:      &  mlocal,nlocal,                                                  &
 43:      &  info(MAT_INFO_BLOCK_SIZE),info(MAT_INFO_NZ_ALLOCATED),          &
 44:      &  info(MAT_INFO_NZ_USED),info(MAT_INFO_NZ_UNNEEDED),              &
 45:      &  info(MAT_INFO_MEMORY),info(MAT_INFO_ASSEMBLIES),                &
 46:      &  info(MAT_INFO_MALLOCS)

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

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

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

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

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

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

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