! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  SLEPc - Scalable Library for Eigenvalue Problem Computations
!  Copyright (c) 2002-2007, Universidad Politecnica de Valencia, Spain
!
!  This file is part of SLEPc. See the README file for conditions of use
!  and additional information.
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
!  This example can be found in src/examples
!
!  Program usage: mpirun -np n ex1f [-help] [-n <n>] [all SLEPc options] 
!
!  Description: Simple example that solves an eigensystem with the EPS object.
!  The standard symmetric eigenvalue problem to be solved corresponds to the 
!  Laplacian operator in 1 dimension. 
!
!  The command line options are:
!    -n <n>, where <n> = number of grid points = matrix size
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

program main
  implicit none

#include "finclude/petsc.h"
#include "finclude/petscvec.h"
#include "finclude/petscmat.h"
#include "finclude/slepc.h"
#include "finclude/slepceps.h"

  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
  !     Declarations
  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
  !     Variables:
  !     A     operator matrix
  !     eps   eigenproblem solver context

  Mat               :: A
  EPS               :: eps
  EPSType           :: type
  PetscReal         :: tol, error
  PetscScalar       :: kr, ki
  PetscTruth        :: flg
  PetscScalar       :: value(3)
  integer :: rank, n, nev, ierr, maxit, i, its, nconv
  integer :: col(3), Istart, Iend

  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  !     Beginning of program
  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
  
  call SlepcInitialize(PETSC_NULL_CHARACTER,ierr)
  call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
  n = 30
  call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-n',n,flg,ierr)

  if (rank .eq. 0) then
     write(*,100) n
  end if

100 format (/'1-D Laplacian Eigenproblem, n =',I3,' (Fortran)')

  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
  !     Compute the operator matrix that defines the eigensystem, Ax=kx
  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

  call MatCreate(PETSC_COMM_WORLD,A,ierr)
  call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr)
  call MatSetFromOptions(A,ierr)

  call MatGetOwnershipRange(A,Istart,Iend,ierr)

  if (Istart .eq. 0) then 
     i = 0
     col(1) = 0
     col(2) = 1
     value(1) =  2.0
     value(2) = -1.0
     call MatSetValues(A,1,i,2,col,value,INSERT_VALUES,ierr)
     Istart = Istart+1
  end if

  if (Iend .eq. n) then 
     i = n-1
     col(1) = n-2
     col(2) = n-1
     value(1) = -1.0
     value(2) =  2.0
     call MatSetValues(A,1,i,2,col,value,INSERT_VALUES,ierr)
     Iend = Iend-1
  end if

  value(1) = -1.0
  value(2) =  2.0
  value(3) = -1.0

  do i=Istart,Iend-1
     col(1) = i-1
     col(2) = i
     col(3) = i+1
     call MatSetValues(A,1,i,3,col,value,INSERT_VALUES,ierr)
  end do

  call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
  call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)

  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
  !     Create the eigensolver and display info
  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

  !     Create eigensolver context
  call EPSCreate(PETSC_COMM_WORLD,eps,ierr)

  !     Set operators. In this case, it is a standard eigenvalue problem
  call EPSSetOperators(eps,A,PETSC_NULL_OBJECT,ierr)
  call EPSSetProblemType(eps,EPS_HEP,ierr)

  !     Set solver parameters at runtime
  call EPSSetFromOptions(eps,ierr)

  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
  !     Solve the eigensystem
  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

  call EPSSolve(eps,ierr) 

  call EPSGetIterationNumber(eps,its,ierr)
  if (rank .eq. 0) then
     write(*,110) its
  end if
110 format (/' Number of iterations of the method:',I4)

  !     Optional: Get some information from the solver and display it

  call EPSGetType(eps,type,ierr)
  if (rank .eq. 0) then
     write(*,120) type
  end if
120 format (' Solution method: ',A)

  call EPSGetDimensions(eps,nev,PETSC_NULL_INTEGER,ierr)
  if (rank .eq. 0) then
     write(*,130) nev
  end if
130 format (' Number of requested eigenvalues:',I2)

  call EPSGetTolerances(eps,tol,maxit,ierr)
  if (rank .eq. 0) then
     write(*,140) tol, maxit
  end if
140 format (' Stopping condition: tol=',1P,E10.4,', maxit=',I4)

  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
  !     Display solution and clean up
  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

  !     Get number of converged eigenpairs
  call EPSGetConverged(eps,nconv,ierr)
  if (rank .eq. 0) then
     write(*,150) nconv
  end if
150 format (' Number of converged eigenpairs:',I2/)

  !     Display eigenvalues and relative errors
  if (nconv.gt.0) then

     if (rank .eq. 0) then
        write(*,*) '         k          ||Ax-kx||/||kx||'
        write(*,*) ' ----------------- ------------------'
     end if

     do i=0,nconv-1
        !     Get converged eigenpairs: i-th eigenvalue is stored in kr 
        !     (real part) and ki (imaginary part)
        call EPSGetEigenpair(eps,i,kr,ki,PETSC_NULL,PETSC_NULL,ierr)

        !     Compute the relative error associated to each eigenpair
        call EPSComputeRelativeError(eps,i,error,ierr)

        if (rank .eq. 0) then
           write(*,160) PetscRealPart(kr), error
        end if
160     format (1P,'   ',E12.4,'       ',E12.4)
     end do

     if (rank .eq. 0) then
        write(*,*)
     end if

  end if

  !     Free work space
  call EPSDestroy(eps,ierr)
  call MatDestroy(A,ierr)

  call SlepcFinalize(ierr)
end program main