program main

!*******************************************************************************
!
!! SSSIMP is a simple program to call ARPACK for a symmetric eigenproblem.
!
!  Discussion:
!
!    This example program is intended to illustrate the
!    simplest case of using ARPACK in considerable detail.
!    This code may be used to understand basic usage of ARPACK
!    and as a template for creating an interface to ARPACK.
!
!    This code shows how to use ARPACK to find a few eigenvalues
!    LAMBDA and corresponding eigenvectors X for the standard
!    eigenvalue problem:
!
!      A * X = LAMBDA * X
!
!    where A is an N by N real symmetric matrix.
!
!    The main points illustrated here are:
!
!    1) How to declare sufficient memory to find NEV
!       eigenvalues of largest magnitude.  Other options
!       are available.
!
!    2) Illustration of the reverse communication interface
!       needed to utilize the top level ARPACK routine SSAUPD
!       that computes the quantities needed to construct
!       the desired eigenvalues and eigenvectors (if requested).
!
!    3) How to extract the desired eigenvalues and eigenvectors
!       using the ARPACK routine SSEUPD.
!
!    The only things that must be supplied in order to use this
!    routine on your problem is:
!
!    * to change the array dimensions appropriately, 
!    * to specify WHICH eigenvalues you want to compute
!    * to supply a matrix-vector product
!      w <- A * v
!      in place of the call to AV( ) below.
!
!    Once usage of this routine is understood, you may wish to explore
!    the other available options to improve convergence, to solve generalized
!    problems, etc.  Look at the file ex-sym.doc in DOCUMENTS directory.
!
!
!    For this problem, we want to solve A*x = lambda*x in regular mode.
!
!    A is derived from the central difference discretization
!    of the 2-dimensional Laplacian on the unit square with
!    zero Dirichlet boundary condition.
!
!    OP = A  and  B = I.
!    Assume "call av (n,x,y)" computes y = A*x
!    Use mode 1 of SSAUPD.
!
!  Author:
!
!    Richard Lehoucq, Danny Sorensen, Chao Yang,
!    Department of Computational and Applied Mathematics,
!    Rice University,
!    Houston, Texas.
!
!  Storage:
! 
!    The maximum dimensions for all arrays are set here to accommodate 
!    a problem size of N <= MAXN
!
!    NEV is the number of eigenvalues requested.
!    See specifications for ARPACK usage below.
!
!    NCV is the largest number of basis vectors that will be used in 
!    the Implicitly Restarted Arnoldi Process.  Work per major iteration is
!    proportional to N*NCV*NCV.
!
!    You must set: 
! 
!    MAXN:   Maximum dimension of the A allowed. 
!    MAXNEV: Maximum NEV allowed. 
!    MAXNCV: Maximum NCV allowed. 
!
  implicit none

  integer, parameter :: maxn = 256
  integer, parameter :: maxnev = 10
  integer, parameter :: maxncv = 25

  integer, parameter :: ldv = maxn

  intrinsic abs
  real ax(maxn)
  character bmat  
  real d(maxncv,2)
  integer ido
  integer ierr
  integer info
  integer iparam(11)
  integer ipntr(11)
  integer ishfts
  integer j
  integer lworkl
  integer maxitr
  integer mode1
  integer n
  integer nconv
  integer ncv
  integer nev
  integer nx
  real resid(maxn)
  logical rvec
  external saxpy
  logical select(maxncv)
  real sigma
  real, external :: snrm2
  real tol
  real v(ldv,maxncv)
  character ( len = 2 ) which
  real workl(maxncv*(maxncv+8))
  real workd(3*maxn)
  real, parameter :: zero = 0.0E+00
!
!  The following include statement and assignments control trace output 
!  from the internal actions of ARPACK.  See debug.doc in the
!  DOCUMENTS directory for usage.  
!
!  Initially, the most useful information will be a breakdown of
!  time spent in the various stages of computation given by setting 
!  msaupd = 1.
!
  include 'debug.h'

  call timestamp ( )
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'SSSIMP:'
  write ( *, '(a)' ) '  FORTRAN90 version'
  write ( *, '(a)' ) '  A simple ARPACK calling program.'
  write ( *, '(a)' ) '  This program defines an eigenproblem for a'
  write ( *, '(a)' ) '  symmetric matrix.'

  ndigit = -3
  logfil = 6
  msgets = 0
  msaitr = 0
  msapps = 0
  msaupd = 1
  msaup2 = 0
  mseigt = 0
  mseupd = 0
!
!  Set dimensions for this problem.
!
  nx = 10
  n = nx * nx
!
!  Specifications for ARPACK usage are set below:
!
!  1) NEV = 4 asks for 4 eigenvalues to be computed.                            !
!  2) NCV = 20 sets the length of the Arnoldi factorization.
!
!  3) This is a standard problem(indicated by bmat  = 'I')
!
!  4) Ask for the NEV eigenvalues of largest magnitude
!     (indicated by which = 'LM')
!
!  See documentation in SSAUPD for the other options SM, LA, SA, LI, SI.
!
!  NEV and NCV must satisfy the following conditions:
!
!    NEV <= MAXNEV
!    NEV + 1 <= NCV <= MAXNCV
!
  nev = 4
  ncv = 20
  bmat = 'I'
  which = 'LM'

  if ( maxn < n ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'SSSIMP - Fatal error!'
    write ( *, '(a)' ) '  N is greater than MAXN '
    stop
  else if ( maxnev < nev ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'SSSIMP - Fatal error!'
    write ( *, '(a)' ) '  NEV is greater than MAXNEV '
    stop
  else if ( maxncv < ncv ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'SSSIMP - Fatal error!'
    write ( *, '(a)' ) '  NCV is greater than MAXNCV '
    stop
  end if
!
!  Specification of stopping rules and initial
!  conditions before calling SSAUPD
!
!  TOL determines the stopping criterion.  Expect
!    abs(lambdaC - lambdaT) < TOL*abs(lambdaC)
!  computed   true
!  If TOL <= 0, then TOL <- macheps (machine precision) is used.
!
!  IDO is the REVERSE COMMUNICATION parameter
!  used to specify actions to be taken on return
!  from SSAUPD. (See usage below.)
!  It MUST initially be set to 0 before the first
!  call to SSAUPD.
!
!  INFO on entry specifies starting vector information
!  and on return indicates error codes
!  Initially, setting INFO=0 indicates that a
!  random starting vector is requested to 
!  start the ARNOLDI iteration.  Setting INFO to
!  a nonzero value on the initial call is used 
!  if you want to specify your own starting 
!  vector. (This vector must be placed in RESID.)
!
!  The work array WORKL is used in SSAUPD as workspace.  Its dimension
!  LWORKL is set as illustrated below. 
!
  lworkl = ncv * ( ncv + 8 )
  tol = zero
  info = 0
  ido = 0
!
!  Specification of Algorithm Mode:
!
!  This program uses the exact shift strategy
!  (indicated by setting PARAM(1) = 1).
!
!  IPARAM(3) specifies the maximum number of Arnoldi iterations allowed.  
!
!  Mode 1 of SSAUPD is used (IPARAM(7) = 1). 
!
!  All these options can be changed by the user.  For details see the
!  documentation in SSAUPD.
!
  ishfts = 1
  maxitr = 300
  mode1 = 1

  iparam(1) = ishfts

  iparam(3) = maxitr

  iparam(7) = mode1
!
!  MAIN LOOP (Reverse communication loop)
!
!  Repeatedly call SSAUPD and take actions indicated by parameter 
!  IDO until convergence is indicated or MAXITR is exceeded.
!
  do

    call ssaupd ( ido, bmat, n, which, nev, tol, resid, &
      ncv, v, ldv, iparam, ipntr, workd, workl, &
      lworkl, info )

    if ( ido /= -1 .and. ido /= 1 ) then
      exit
    end if
!
!  Perform matrix vector multiplication
!
!    y <--- OP*x
!
!  The user supplies a matrix-vector multiplication routine that takes
!  workd(ipntr(1)) as the input, and return the result to workd(ipntr(2)).
!
    call av ( nx, workd(ipntr(1)), workd(ipntr(2)) )

   end do
!
!  Either we have convergence or there is an error.
!
  if ( info < 0 ) then
!
!  Error message. Check the documentation in SSAUPD.
!
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'SSSIMP - Fatal error!'
    write ( *, '(a,i6)' ) '  Error with SSAUPD, INFO = ', info
    write ( *, '(a)' ) '  Check documentation in SSAUPD.'

  else
!
!  No fatal errors occurred.
!  Post-Process using SSEUPD.
!
!  Computed eigenvalues may be extracted.
!
!  Eigenvectors may be also computed now if
!  desired.  (indicated by rvec = .true.)
!
!  The routine SSEUPD now called to do this
!  post processing (Other modes may require
!  more complicated post processing than mode1.)
!
    rvec = .true.

    call sseupd ( rvec, 'All', select, d, v, ldv, sigma, &
      bmat, n, which, nev, tol, resid, ncv, v, ldv, &
      iparam, ipntr, workd, workl, lworkl, ierr )
!
!  Eigenvalues are returned in the first column of the two dimensional 
!  array D and the corresponding eigenvectors are returned in the first 
!  NCONV (=IPARAM(5)) columns of the two dimensional array V if requested.
!
!  Otherwise, an orthogonal basis for the invariant subspace corresponding 
!  to the eigenvalues in D is returned in V.
!
    if ( ierr /= 0 ) then

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'SSSIMP - Fatal error!'
      write ( *, '(a,i6)' ) '  Error with SSEUPD, IERR = ', ierr
      write ( *, '(a)' ) '  Check the documentation of SSEUPD.'
!
!  Compute the residual norm
!
!    ||  A*x - lambda*x ||
! 
!  for the NCONV accurately computed eigenvalues and 
!  eigenvectors.  (iparam(5) indicates how many are 
!  accurate to the requested tolerance)
!
    else

      nconv =  iparam(5)

      do j = 1, nconv
        call av ( nx, v(1,j), ax )
        call saxpy ( n, -d(j,1), v(1,j), 1, ax, 1 )
        d(j,2) = snrm2 ( n, ax, 1)
        d(j,2) = d(j,2) / abs ( d(j,1) )
      end do
!
!  Display computed residuals.
!
      call smout ( 6, nconv, 2, d, maxncv, -6, &
        'Ritz values and relative residuals' )

    end if
!
!  Print additional convergence information.
!
    if ( info == 1) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  Maximum number of iterations reached.'
    else if ( info == 3) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  No shifts could be applied during implicit' &
        // ' Arnoldi update, try increasing NCV.'
    end if

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'SSSIMP:'
    write ( *, '(a)' ) '====== '
    write ( *, '(a)' ) ' '
    write ( *, '(a,i6)' ) '  Size of the matrix is ', n
    write ( *, '(a,i6)' ) '  The number of Ritz values requested is ', nev
    write ( *, '(a,i6)' ) &
      '  The number of Arnoldi vectors generated (NCV) is ', ncv
    write ( *, '(a)' ) '  What portion of the spectrum: ' // which
    write ( *, '(a,i6)' ) &
      '  The number of converged Ritz values is ', nconv
    write ( *, '(a,i6)' ) &
      '  The number of Implicit Arnoldi update iterations taken is ', iparam(3)
    write ( *, '(a,i6)' ) '  The number of OP*x is ', iparam(9)
    write ( *, '(a,g14.6)' ) '  The convergence criterion is ', tol

  end if

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'SSSIMP:'
  write ( *, '(a)' ) '  Normal end of execution.'

  write ( *, '(a)' ) ' '
  call timestamp ( )

  stop
end
subroutine av ( nx, v, w )

!*******************************************************************************
!
!! AV computes w <- A * V where A is a discretized Laplacian.
!
!  Discussion:
!
!    The matrix used is the 2 dimensional discrete Laplacian on unit
!    square with zero Dirichlet boundary condition.
!
!    Compute w <--- OP * v, where OP is the nx*nx by nx*nx block
!    tridiagonal matrix
!
!      | T -I           |
!      |-I  T -I        |
!      |   -I  T        |
!      |        ...  -I |
!      |           -I T |
!
!  Parameters:
!
!    Input, integer NX, the length of the vectors.
!
!    Input, real V(NX), the vector to be operated on by A.
!
!    Output, real W(NX), the result of A*V.
!
  implicit none

  integer nx

  real h2
  integer j
  integer lo
  integer n2
  real, parameter :: one = 1.0E+00
  real v(nx*nx)
  real w(nx*nx)

  call tv ( nx, v(1), w(1) )

  call saxpy ( nx, -one, v(nx+1), 1, w(1), 1 )

  do j = 2, nx-1
    lo = (j-1) * nx
    call tv ( nx, v(lo+1), w(lo+1) )
    call saxpy ( nx, -one, v(lo-nx+1), 1, w(lo+1), 1 )
    call saxpy ( nx, -one, v(lo+nx+1), 1, w(lo+1), 1 )
  end do

  lo = (nx-1) * nx
  call tv ( nx, v(lo+1), w(lo+1) )
  call saxpy ( nx, -one, v(lo-nx+1), 1, w(lo+1), 1 )
!
!  Scale the vector W by (1/H^2), where H is the mesh size.
!
  n2 = nx * nx
  h2 = one / real ( ( nx + 1 ) * ( nx + 1 ) )
  call sscal ( n2, one/h2, w, 1 )

  return
end
subroutine tv ( nx, x, y )

!*******************************************************************************
!
!! TV computes y <-- T*x.
!
!  Discussion:
!
!    The full matrix A is the 2 dimensional discrete Laplacian on unit
!    square with zero Dirichlet boundary condition.
!
!    A can be represented as the nx*nx by nx*nx block tridiagonal matrix
!
!      | T -I          |
!      |-I  T -I       |
!      |   -I  T       |
!      |        ...  -I|
!      |           -I T|
!
!    It is the job of this routine to compute products of the form T*x.
!
!  Parameters:
!
!    Input, integer NX, the length of the vectors.
!
!    Input, real X(NX), the vector to be multiplied by T.
!
!    Output, real Y(NX), the product T*X.
!
  implicit none

  integer nx

  integer j
  real dd
  real dl
  real du
  real x(nx)
  real y(nx)
!
!  Compute the matrix vector multiplication y<---T*x
!  where T is a nx by nx tridiagonal matrix with DD on the
!  diagonal, DL on the subdiagonal, and DU on the superdiagonal.
!
  dd  = 4.0E+00
  dl  = -1.0E+00
  du  = -1.0E+00

  y(1) =  dd * x(1) + du * x(2)

  do j = 2, nx-1
     y(j) = dl * x(j-1) + dd * x(j) + du * x(j+1)
  end do

  y(nx) = dl * x(nx-1) + dd * x(nx)

  return
end