text:
!---------------------------------------------------------------------!
!                                PHAML                                !
!                                                                     !
! The Parallel Hierarchical Adaptive MultiLevel code for solving      !
! linear elliptic partial differential equations of the form          !
! (PUx)x + (QUy)y + RU = F on 2D polygonal domains with mixed         !
! boundary conditions, and eigenvalue problems where F is lambda*U.   !
!                                                                     !
! PHAML is public domain software.  It was produced as part of work   !
! done by the U.S. Government, and is not subject to copyright in     !
! the United States.                                                  !
!                                                                     !
!     William F. Mitchell                                             !
!     Applied and Computational Mathematics Division                  !
!     National Institute of Standards and Technology                  !
!     william.mitchell@nist.gov                                       !
!     http://math.nist.gov/phaml                                      !
!                                                                     !
!---------------------------------------------------------------------!

module eigen

!----------------------------------------------------

!----------------------------------------------------
! This module contains routines for solve eigenproblems
!
! communication tags in this module are of the form 2xxx
!----------------------------------------------------

!----------------------------------------------------
! Other modules used:

use linsystype_mod
use global
use message_passing
use gridtype_mod
use petsc_interf
use mumps_interf
use hypre_interf
use superlu_interf
use linsys_util
use sort_mod
use krylov
use lapack_solve
use hbmg
use make_linsys
!----------------------------------------------------

implicit none
private
public eigen_arpack, eigen_blopex

contains

!          ------------
subroutine eigen_arpack(grid,procs,linsys,io_cntl,solver_cntl,still_sequential)
!          ------------

!----------------------------------------------------
! This routine solves the eigenvalue problem Ax = lambda*Mx where A is
! stiffness and M is mass.  It uses ARPACK to solve for the largest
! eigenvalue(s) and eigenvector(s) of
!
!                -1
! (A - lambda0*M)  M x = mu*x
!
! which gives lambda = 1/mu + lambda0 to give the eigenvalues closest to
! lambda0 (but they are computed by the Raleigh quotient rather than this
! formula).  If lambda0 == -huge(0), then rmin (the minimum value of the
! coefficient of u in the PDE, which should be smaller than the smallest
! eigenvalue) is used to get the smallest eigenvalues.
!
! This is the traditional "shift and invert" spectral transformation.  In
! ARPACK notation
!                      -1
!  OP = (A - lambda0*M)  M  and B = M
!
!----------------------------------------------------

!----------------------------------------------------
! Dummy arguments

type(grid_type), intent(inout) :: grid
type(proc_info), intent(in) :: procs
type(linsys_type), intent(inout), target :: linsys
type(io_options), intent(in) :: io_cntl
type(solver_options), intent(in) :: solver_cntl
logical, intent(in) :: still_sequential

!----------------------------------------------------
! Local variables:

! choose between using the symmetric solvers and nonsymmetric solvers

logical, parameter :: symm = .false.

real(my_real), allocatable :: arnoldi_vec(:,:), arp_workd(:), arp_resid(:), &
                 arp_workl(:), dr(:), di(:), evec(:,:), workev(:)
integer :: nev, ncv, lworkl, kev
integer :: arp_param(11), arp_pntr(14)
integer, allocatable :: iperm(:)
real(my_real) :: tol, maxsoln, minsoln
real(my_real) :: l2resid, l2rhs, solnMnorm, soln2norm, sums(2)
integer :: ido, info, iter, nproc, comm, neqown, allocstat, neq_total, &
           ndirch, neq, neq_full, objtype,brank,srank,objlid,loc_neq, ii, jj, &
           nresid_comp, lev
character(len=1) :: bmat, howmany
character(len=2) :: which
logical :: rvec, select(solver_cntl%arpack_cntl%ncv)
integer :: i, j, ss, eq
real(my_real) :: tempvec(linsys%neq_vert+linsys%neq_edge+linsys%neq_face), &
                 Ax(linsys%neq_vert+linsys%neq_edge+linsys%neq_face), &
                 Mx(linsys%neq_vert+linsys%neq_edge+linsys%neq_face)
real(my_real), pointer :: tmpptr(:), block(:,:)
integer, pointer :: ipiv(:)
real(my_real) :: my_resid(linsys%neq)
real(my_real) :: numer, denom, lambda0
logical :: i_made_lapack, i_made_petsc, i_made_mumps, i_made_superlu, &
           i_made_hypre
type(linsys_type) :: full_matrix
integer :: ni, nr, proc
integer, pointer :: irecv(:)
real(my_real), pointer :: rrecv(:)

!----------------------------------------------------
! Begin executable code

! in case we abort early

grid%arpack_iter = 0
grid%arpack_nconv = 0
grid%arpack_numop = 0
grid%arpack_numopb = 0
grid%arpack_numreo = 0
grid%arpack_info = 0

! master hangs around just to see if it needs to print a warning

if (my_proc(procs) == MASTER) then
   do
      call phaml_recv(procs,proc,irecv,ni,rrecv,nr,2201)
      select case (irecv(1))
      case (1) ! all done
         deallocate(irecv,stat=allocstat)
         exit
      case (2) ! ARPACK eigenvalue does not agree with Raleigh quotient
         call warning("Eigenvalue returned by ARPACK is not close to Raleigh quotient; may be wrong", &
                      "Eigenvalue number, eigenvalue and Raleigh quotient are:", &
                      (/irecv(2)/),rrecv)
         deallocate(irecv,rrecv,stat=allocstat)
      end select
   end do
   return
endif

! use linsys for solving the (statically condensed) shifted linear systems,
! and full_matrix for operations that need the whole matrix.  full_matrix
! points to entries in linsys with the following pointers/values different:
!  linsys%end_row => end_row_edge
!  linsys%matrix_val => condensed => shifted
!  linsys%rhs => rhs_cond
!  linsys%neq = neq_vert+neq_edge

full_matrix = linsys

full_matrix%end_row => linsys%end_row_face
full_matrix%matrix_val => linsys%stiffness
full_matrix%rhs => linsys%rhs_nocond
full_matrix%neq = linsys%neq_vert + linsys%neq_edge + linsys%neq_face

! convenience variables

nproc = num_proc(procs)
comm = slaves_communicator(procs)
nev = solver_cntl%num_eval
neq = linsys%neq
neq_full = full_matrix%neq
ss = grid%system_size
if (solver_cntl%lambda0 == -huge(0.0_my_real)) then
   lambda0 = linsys%rmin
else
   lambda0 = solver_cntl%lambda0
endif

! Make sure there is space in grid for the requested number of eigenvalues

if (size(grid%eigenvalue) /= nev) then
   deallocate(grid%eigenvalue,grid%eigenprob_l2_resid, &
              grid%eigenprob_variance,stat=allocstat)
   allocate(grid%eigenvalue(nev),grid%eigenprob_l2_resid(nev), &
            grid%eigenprob_variance(nev),stat=allocstat)
   if (allocstat /= 0) then
      ierr = ALLOC_FAILED
      call fatal("memory allocation failed in eigen_arpack",procs=procs)
      return
   endif
   grid%eigenvalue = 0.0_my_real
   grid%eigenprob_l2_resid = 0.0_my_real
   grid%eigenprob_variance = 0.0_my_real
endif

! copy the solutions from the grid to the linear system to return the
! old solution if we abort

do i=1,neq_full
   call eq_to_grid(full_matrix,full_matrix%gid(i),objtype,brank,srank, &
                   objlid,grid)
   select case (objtype)
   case (VERTEX_ID)
      full_matrix%solution(i) = grid%vertex_solution(objlid,srank,1)
      if (nev > 1) then
         full_matrix%evecs(i,:) = grid%vertex_solution(objlid,srank,2:nev)
      endif
   case (EDGE_ID)
      full_matrix%solution(i) = grid%edge(objlid)%solution(brank,srank,1)
      if (nev > 1) then
         do j=1,nev-1
            full_matrix%evecs(i,j) = grid%edge(objlid)%solution(brank,srank,j)
         end do
      endif
   case (ELEMENT_ID)
      full_matrix%solution(i) = grid%element(objlid)%solution(brank,srank,1)
      if (nev > 1) then
         do j=1,nev-1
            full_matrix%evecs(i,j) = grid%element(objlid)%solution(brank,srank,j)
         end do
      endif
   end select
end do

! all processors must own an equation for ARPACK to work

if (phaml_global_min(procs,grid%nvert_own,2000) == 0) then
   call warning("Some processor owns no equations; ARPACK would fail.", &
                "Returning old solution.")
   return
endif

! if there aren't enough equations, reduce number of Arnoldi vectors

ncv = solver_cntl%arpack_cntl%ncv

neqown = 0
ndirch = 0
do i=1,neq_full
   if (full_matrix%iown(i)) then
      neqown = neqown + 1
      if (full_matrix%equation_type(i) == DIRICHLET) ndirch = ndirch + 1
   endif
end do
if (still_sequential) then
   neq_total = neqown - ndirch
else
   neq_total = phaml_global_sum(procs,neqown-ndirch,2001)
endif
if (neq_total < ncv) ncv = max(1,neq_total-1)

! if ncv is too small, reduce number of eigenvalues

if (ncv - nev < 2) then
   nev = ncv - 2
   call warning( &
     "Not enough equations or Arnoldi vectors for requested number of eigenvalues.", &
     "Computing", intlist=(/nev/))
endif

! allocate memory for ARPACK

lworkl = 3*ncv*ncv + 6*ncv
allocate(arnoldi_vec(neqown,ncv), arp_workd(3*neqown), arp_resid(neqown), &
         arp_workl(lworkl), dr(nev+1), di(nev+1), evec(neqown,nev+1), &
         workev(3*ncv), iperm(nev), stat=allocstat)
if (allocstat /= 0) then
   ierr = ALLOC_FAILED
   call fatal("memory allocation failed in eigen_arpack",procs=procs)
   return
endif

! initalize convergence information

grid%eigen_linsys_max_l2_resid = 0.0_my_real
grid%eigen_linsys_ave_l2_resid = 0.0_my_real
nresid_comp = 0
grid%eigenprob_l2_resid = 0.0_my_real
grid%eigenprob_variance = 0.0_my_real

! create any needed alternate storage formats for solution methods

i_made_lapack = .false.
i_made_petsc = .false.
i_made_mumps = .false.
i_made_superlu = .false.
i_made_hypre = .false.

select case (solver_cntl%solver)

case (LAPACK_INDEFINITE_SOLVER)
   if (linsys%lapack_symm_band_exists) then
      call destroy_lapack_band(linsys%lapack_mat)
      linsys%lapack_symm_band_exists = .false.
   endif
   if (.not. linsys%lapack_gen_band_exists) then
      call make_lapack_gen_band(linsys%nlev+1,linsys,linsys%lapack_mat)
      linsys%lapack_gen_band_exists = .true.
      i_made_lapack = .true.
   endif

case (LAPACK_SPD_SOLVER)
   if (linsys%lapack_gen_band_exists) then
      call destroy_lapack_band(linsys%lapack_mat)
      linsys%lapack_gen_band_exists = .false.
   endif
   if (.not. linsys%lapack_symm_band_exists) then
      call make_lapack_symm_band(linsys%nlev+1,linsys,linsys%lapack_mat)
      if (ierr /= NO_ERROR) return
      linsys%lapack_symm_band_exists = .true.
      i_made_lapack = .true.
   endif

case (MUMPS_SPD_SOLVER, MUMPS_GEN_SOLVER, MUMPS_NONSYM_SOLVER)
   if (.not. linsys%mumps_matrix_exists) then
      call create_mumps_linear_system(linsys%mumps_matrix,linsys,procs, &
                                      solver_cntl,still_sequential)
      linsys%mumps_matrix_exists = .true.
      i_made_mumps = .true.
   endif

case (SUPERLU_SOLVER)
   if (.not. linsys%superlu_matrix_exists) then
      call create_superlu_linear_system(linsys%superlu_matrix,linsys, &
                                        procs,still_sequential)
      linsys%superlu_matrix_exists = .true.
      i_made_superlu = .true.
   endif

case (PETSC_RICHARDSON_SOLVER, PETSC_CHEBYCHEV_SOLVER, PETSC_CG_SOLVER, &
      PETSC_GMRES_SOLVER,      PETSC_TCQMR_SOLVER,     PETSC_BCGS_SOLVER,&
      PETSC_CGS_SOLVER,        PETSC_TFQMR_SOLVER,     PETSC_CR_SOLVER, &
      PETSC_LSQR_SOLVER,       PETSC_BICG_SOLVER)
   if (.not. linsys%petsc_matrix_exists) then
      if (solver_cntl%petsc_matrix_free) then
         call create_petsc_linear_system_mf(linsys,linsys%petsc_matrix, &
                                         still_sequential)
      else
         call create_petsc_linear_system(linsys,linsys%petsc_matrix, &
                                         still_sequential,procs)
      endif
      linsys%petsc_matrix_exists = .true.
      i_made_petsc = .true.
   endif

case (HYPRE_BOOMERAMG_SOLVER, HYPRE_PCG_SOLVER, HYPRE_GMRES_SOLVER)
   if (.not. linsys%hypre_matrix_exists) then
      call create_hypre_linear_system(linsys%hypre_matrix,linsys,procs, &
                                      still_sequential)
      linsys%hypre_matrix_exists = .true.
      i_made_hypre = .true.
   endif

end select

! chose between closest to lambda0 (LM), larger than lambda0 (LR),
! and smaller than lambda0 (SR) (assumes symm=.false., different if true)

select case (solver_cntl%lambda0_side)
case (EIGEN_BOTH)
   which = "LM"                         ! compute largest magnitude eigenvalues
case (EIGEN_RIGHT)
   which = "LR"                         ! compute largest real part eigenvalues
case (EIGEN_LEFT)
   which = "SR"                         ! compute smallest real part e-values
case default
   ierr = USER_INPUT_ERROR
   call fatal("bad value for lambda0_side in eigen_arpack", &
              intlist=(/solver_cntl%lambda0_side/),procs=procs)
   stop
end select

! set parameters for ARPACK

ido = 0                                 ! indicates first call
bmat = "G"                              ! generalized eigenvalue problem
tol = solver_cntl%arpack_cntl%tol       ! stopping criteria
arp_param = 0                           ! set unused entries to 0
arp_param(1) = 1                        ! shifts are not provided by the user
arp_param(3) = solver_cntl%arpack_cntl%maxit ! maximum number of iterations
arp_param(7) = 3                        ! shift-invert mode

! set initial guess for eigenvector
! TEMP would it be better to use some combination of existing eigenvectors
! instead of the first one?

call all_to_mine(full_matrix,full_matrix%solution(1:),arp_resid)
if (maxval(abs(arp_resid)) < 10*epsilon(0.0_my_real)) then
   info = 0                                 ! random initial residual vector
else
   info = 1                                 ! nonrandom initial residual vector
endif

! main loop

iter = 0
do
   iter = iter + 1

! call ARPACK routine [p| ][s|d][s|n]aupd for
! [parallel|sequential][single|double] precision [symmetric|nonsymmetric]
! eigensolver

   if (my_real == kind(1.0e0)) then
      if (nproc > 1 .and. .not. still_sequential) then
       if (symm) then
         call pssaupd(comm,ido,bmat,neqown,which,nev,tol,arp_resid,ncv, &
                     arnoldi_vec,neqown,arp_param,arp_pntr,arp_workd, &
                     arp_workl,lworkl,info)
       else
         call psnaupd(comm,ido,bmat,neqown,which,nev,tol,arp_resid,ncv, &
                     arnoldi_vec,neqown,arp_param,arp_pntr,arp_workd, &
                     arp_workl,lworkl,info)
       endif
      else
       if (symm) then
         call ssaupd(ido,bmat,neqown,which,nev,tol,arp_resid,ncv, &
                     arnoldi_vec,neqown,arp_param,arp_pntr,arp_workd, &
                     arp_workl,lworkl,info)
       else
         call snaupd(ido,bmat,neqown,which,nev,tol,arp_resid,ncv, &
                     arnoldi_vec,neqown,arp_param,arp_pntr,arp_workd, &
                     arp_workl,lworkl,info)
       endif
      endif
   elseif (my_real == kind(1.0d0)) then
      if (nproc > 1 .and. .not. still_sequential) then
       if (symm) then
         call pdsaupd(comm,ido,bmat,neqown,which,nev,tol,arp_resid,ncv, &
                     arnoldi_vec,neqown,arp_param,arp_pntr,arp_workd, &
                     arp_workl,lworkl,info)
       else
         call pdnaupd(comm,ido,bmat,neqown,which,nev,tol,arp_resid,ncv, &
                     arnoldi_vec,neqown,arp_param,arp_pntr,arp_workd, &
                     arp_workl,lworkl,info)
       endif
      else
       if (symm) then
         call dsaupd(ido,bmat,neqown,which,nev,tol,arp_resid,ncv, &
                     arnoldi_vec,neqown,arp_param,arp_pntr,arp_workd, &
                     arp_workl,lworkl,info)
       else
         call dnaupd(ido,bmat,neqown,which,nev,tol,arp_resid,ncv, &
                     arnoldi_vec,neqown,arp_param,arp_pntr,arp_workd, &
                     arp_workl,lworkl,info)
       endif
      endif
   else
      ierr = PHAML_INTERNAL_ERROR
      call fatal("my_real is neither single nor double precision. Can't call ARPACK",procs=procs)
      return
   endif

   if (info /= 0) then
      call warning("info from ARPACK *aupd is nonzero",intlist=(/info/))
   endif

! Determine which task to perform

   select case (ido)

   case (-1,1) ! multiply by (A-lambda0*M)^-1*M
               ! If ido==1, then leave off the multiplication by M

! Copy the vector from ARPACK into rhs

      if (ido == -1) then
         call mine_to_all(full_matrix,arp_workd(arp_pntr(1):arp_pntr(1)+neqown-1), &
                          full_matrix%rhs,0.0_my_real)
      else
         call mine_to_all(full_matrix,arp_workd(arp_pntr(3):arp_pntr(3)+neqown-1), &
                          full_matrix%rhs,0.0_my_real)
      endif

! For FuDoP multigrid, get rhs entries from other processors, making sure they
! match my level of the corresponding vertex

      if (solver_cntl%solver == MG_SOLVER .and. .not. still_sequential .and. &
          nproc > 1) then
         do lev=full_matrix%nlev+1,2,-1
            call basis_change(lev,TO_HIER,full_matrix,skip_matrix=.true.)
         end do
         call sum_fudop_vect(full_matrix%rhs,procs,full_matrix, &
                             2310+iter+1,2320+iter+1,2325+iter+1)
         do lev=2,full_matrix%nlev+1
            call basis_change(lev,TO_NODAL,full_matrix,skip_matrix=.true.)
         end do
! second mine_to_all to fix any changes at my points due to basis changes
         if (ido == -1) then
            call mine_to_all(full_matrix,arp_workd(arp_pntr(1):arp_pntr(1)+neqown-1), &
                             full_matrix%rhs)
         else
            call mine_to_all(full_matrix,arp_workd(arp_pntr(3):arp_pntr(3)+neqown-1), &
                             full_matrix%rhs)
         endif
      endif

! for ido==-1, multiply by M putting the result in rhs

      if (ido == -1) then
         tempvec = full_matrix%rhs
         full_matrix%matrix_val => full_matrix%mass
         call matrix_times_vector(tempvec,full_matrix%rhs,full_matrix,procs, &
                                  still_sequential,2100+iter,2110+iter, &
                                  2115+iter,2120+iter,2130+iter,2135+iter, &
                                  nocomm2=.true.)
! For FuDoP multigrid, get the product for equations owned by other processors,
! making sure it is converted to my level for that equation
         if (solver_cntl%solver == MG_SOLVER .and. &
             .not. still_sequential .and. nproc > 1) then
            do i=1,neq
               if (.not. full_matrix%iown(i)) &
                  full_matrix%rhs(i) = 0.0_my_real
            end do
            do lev=full_matrix%nlev+1,2,-1
               call basis_change(lev,TO_HIER,full_matrix,skip_matrix=.true.)
            end do
            call sum_fudop_vect(full_matrix%rhs,procs,full_matrix, &
                                2330+2*iter,2340+2*iter,2345+2*iter)
            do lev=2,full_matrix%nlev+1
               call basis_change(lev,TO_NODAL,full_matrix,skip_matrix=.true.)
            end do
         endif
      endif

! multiply by (A-lambda0*M)^-1, i.e., call a solver

! perform static condensation on the rhs of full_matrix, putting the
! result in the rhs of linsys; b1 - A21 A2^-1 b2 where 2 is face bases

      if (full_matrix%begin_level(full_matrix%nlev+2) /= &
          full_matrix%begin_level(full_matrix%nlev+3)) then

         linsys%rhs(1:neq) = full_matrix%rhs(1:neq)
         full_matrix%matrix_val => full_matrix%shifted
         eq = full_matrix%begin_level(full_matrix%nlev+2)
         do
            if (eq >= full_matrix%begin_level(full_matrix%nlev+3)) exit
            call eq_to_grid(full_matrix,full_matrix%gid(eq),objtype,brank, &
                            srank,objlid,grid)
            block => full_matrix%elem_block(objlid)%matrix
            ipiv => full_matrix%elem_block(objlid)%ipiv
            loc_neq = full_matrix%elem_block(objlid)%neq
            tempvec(1:loc_neq) = full_matrix%rhs(eq:eq+loc_neq-1)
            if (my_real == kind(0.0)) then
              call sgetrs("N",loc_neq,1,block,loc_neq,ipiv,tempvec,loc_neq,info)
            else
              call dgetrs("N",loc_neq,1,block,loc_neq,ipiv,tempvec,loc_neq,info)
            endif
            if (info /= 0) then
               ierr = PHAML_INTERNAL_ERROR
               call fatal("lapack sgetrs solution failed", &
                          intlist=(/info/),procs=procs)
               stop
            endif
            do i=1,loc_neq
               do j=full_matrix%begin_row(eq+i-1),full_matrix%end_row(eq+i-1)
                  if (full_matrix%column_index(j) == NO_ENTRY) cycle
                  if (full_matrix%column_index(j) >= eq .and. &
                      full_matrix%column_index(j) < eq+loc_neq) cycle
                  linsys%rhs(full_matrix%column_index(j)) = &
                                linsys%rhs(full_matrix%column_index(j)) - &
                                get_matval(full_matrix,full_matrix%matrix_val, &
                                full_matrix%column_index(j),eq+i-1)*tempvec(i)
               end do
            end do
            eq = eq + loc_neq
         end do

         if (.not. still_sequential .and. nproc > 1) then
            call sum_fudop_vect(linsys%rhs(1:neq),procs,linsys,2510+iter, &
                                2520+iter,2530+iter)
         endif

      endif

! solve the linear system by the chosen solver

      select case (solver_cntl%solver)
      case (MG_SOLVER)
         call multigrid(grid,procs,linsys,io_cntl,solver_cntl,still_sequential,&
                        no_master=.true.)
      case (CG_SOLVER)
         call phaml_cg(linsys,procs,grid,io_cntl,solver_cntl,still_sequential)
      case (GMRES_SOLVER)
         call phaml_gmres(linsys,procs,grid,io_cntl,solver_cntl, &
                          still_sequential)
      case (LAPACK_INDEFINITE_SOLVER)
         call lapack_indef(linsys%nlev+1,linsys,linsys%lapack_mat)
      case (LAPACK_SPD_SOLVER)
         call lapack_spd(linsys%nlev+1,linsys,linsys%lapack_mat)
      case (MUMPS_SPD_SOLVER, MUMPS_GEN_SOLVER, MUMPS_NONSYM_SOLVER)
         call change_mumps_rhs(linsys%mumps_matrix,linsys,linsys%rhs, &
                               procs,still_sequential)
         call mumps_solve(linsys%mumps_matrix,linsys,procs, &
                          still_sequential,noshadow=.true.)
      case (SUPERLU_SOLVER)
         call change_superlu_rhs(linsys%superlu_matrix,linsys,linsys%rhs, &
                                 procs,still_sequential)
         call superlu_solve(linsys%superlu_matrix,linsys,procs, &
                            still_sequential)
      case (PETSC_RICHARDSON_SOLVER, PETSC_CHEBYCHEV_SOLVER, PETSC_CG_SOLVER, &
            PETSC_GMRES_SOLVER,      PETSC_TCQMR_SOLVER,     PETSC_BCGS_SOLVER,&
            PETSC_CGS_SOLVER,        PETSC_TFQMR_SOLVER,     PETSC_CR_SOLVER, &
            PETSC_LSQR_SOLVER,       PETSC_BICG_SOLVER)
         where (linsys%equation_type == DIRICHLET) linsys%rhs = linsys%solution(1:)
         call change_petsc_rhs(linsys%rhs,linsys%petsc_matrix)
         call petsc_solve(linsys,linsys%petsc_matrix,solver_cntl,io_cntl, &
                          still_sequential,grid,procs)
      case (HYPRE_BOOMERAMG_SOLVER, HYPRE_PCG_SOLVER, HYPRE_GMRES_SOLVER)
         call change_hypre_rhs(linsys%hypre_matrix,linsys,linsys%rhs,procs, &
                               still_sequential)
         call zero_hypre_solution(linsys%hypre_matrix,procs,still_sequential)
         call hypre_solve(linsys%hypre_matrix,linsys,procs,solver_cntl, &
                          still_sequential)
      case default
         ierr = USER_INPUT_ERROR
         call fatal("illegal selection for solver",procs=procs)
         return
      end select

! get solution for unowned equations that neigbor an owned equation

      if (.not. still_sequential) then
         call exchange_neigh_vect(linsys%solution(1:neq),procs,linsys, &
                                  2360+iter,2370+iter,2380+iter)
      endif

! solve for the face bases; x2 = A2^-1 (b2 - A12 x1)
! note full_matrix%solution and linsys%solution point to the same place

      full_matrix%matrix_val => full_matrix%shifted
      eq = full_matrix%begin_level(full_matrix%nlev+2)
      do
         if (eq >= full_matrix%begin_level(full_matrix%nlev+3)) exit
         call eq_to_grid(full_matrix,full_matrix%gid(eq),objtype,brank, &
                         srank,objlid,grid)
         block => full_matrix%elem_block(objlid)%matrix
         ipiv => full_matrix%elem_block(objlid)%ipiv
         loc_neq = full_matrix%elem_block(objlid)%neq
         do i=1,loc_neq
            full_matrix%solution(eq+i-1) = full_matrix%rhs(eq+i-1)
            do j=full_matrix%begin_row(eq+i-1),full_matrix%end_row(eq+i-1)
               if (full_matrix%column_index(j) == NO_ENTRY) cycle
               if (full_matrix%column_index(j) >= eq .and. &
                   full_matrix%column_index(j) < eq+loc_neq) cycle
               full_matrix%solution(eq+i-1) = full_matrix%solution(eq+i-1)- &
                  full_matrix%matrix_val(j)*full_matrix%solution(full_matrix%column_index(j))
            end do
         end do
         if (my_real == kind(0.0)) then
            call sgetrs("N",loc_neq,1,block,loc_neq,ipiv, &
                        full_matrix%solution(eq:eq+loc_neq-1),loc_neq,info)
         else
            call dgetrs("N",loc_neq,1,block,loc_neq,ipiv, &
                        full_matrix%solution(eq:eq+loc_neq-1),loc_neq,info)
         endif
         if (info /= 0) then
            ierr = PHAML_INTERNAL_ERROR
            call fatal("lapack spotrs solution failed for solving for face bases", &
                       intlist=(/info/),procs=procs)
            stop
         endif
         eq = eq + loc_neq
      end do

! compute the residual from the linear system solver (does not include
! face bases)

      my_resid = 0.0_my_real
      call matrix_times_vector(linsys%solution(1:neq),my_resid,linsys, &
                               procs,still_sequential,2410+iter,2420+iter, &
                               2425+iter,2430+iter,2440+iter,2445+iter, &
                               nocomm2=.true.)
      my_resid = my_resid-linsys%rhs(1:neq)
      do i=1,neq
         if (.not. linsys%iown(i)) my_resid(i) = 0.0_my_real
         if (linsys%equation_type(i)==DIRICHLET) my_resid(i) = 0.0_my_real
      end do
      l2resid = sum(my_resid**2)
      l2rhs = sum(linsys%rhs(1:neq)**2, &
                  mask=linsys%equation_type(1:neq)/=DIRICHLET)
      sums = phaml_global_sum(procs,(/l2resid,l2rhs/),2450+iter)
      l2resid = sqrt(sums(1))
      l2rhs = sqrt(sums(2))
      if (l2rhs /= 0.0_my_real) l2resid = l2resid/l2rhs
      grid%eigen_linsys_max_l2_resid = max(grid%eigen_linsys_max_l2_resid, &
                                           l2resid)
      grid%eigen_linsys_ave_l2_resid = grid%eigen_linsys_max_l2_resid + l2resid
      nresid_comp = nresid_comp + 1

! copy the result into ARPACK's workspace

      call all_to_mine(full_matrix,full_matrix%solution(1:neq_full), &
                       arp_workd(arp_pntr(2):arp_pntr(2)+neqown-1))

   case (2) ! multiply by M

! Copy the vector from ARPACK into rhs

      call mine_to_all(full_matrix,arp_workd(arp_pntr(1):arp_pntr(1)+neqown-1), &
                       full_matrix%rhs,0.0_my_real)

! For FuDoP multigrid, get rhs entries from other processors, making sure they
! match my level of the corresponding vertex

      if (solver_cntl%solver == MG_SOLVER .and. .not. still_sequential .and. &
          nproc > 1) then
         do lev=full_matrix%nlev+1,2,-1
            call basis_change(lev,TO_HIER,full_matrix,skip_matrix=.true.)
         end do
         call sum_fudop_vect(full_matrix%rhs,procs,full_matrix, &
                             2350+iter+1,2360+iter+1,2365+iter+1)
         do lev=2,full_matrix%nlev+1
            call basis_change(lev,TO_NODAL,full_matrix,skip_matrix=.true.)
         end do
! second mine_to_all to fix any changes at my points due to basis changes
         call mine_to_all(full_matrix,arp_workd(arp_pntr(1):arp_pntr(1)+neqown-1), &
                          full_matrix%rhs)
      endif

! multiply by M putting the result in rhs

      tempvec = full_matrix%rhs
      full_matrix%matrix_val => full_matrix%mass
      call matrix_times_vector(tempvec,full_matrix%rhs,full_matrix,procs, &
                               still_sequential,2140+iter,2150+iter, &
                               2155+iter,2160+iter,2170+iter,2175+iter, &
                               nocomm2=.true.)

! copy the result into ARPACK's workspace

      call all_to_mine(full_matrix,full_matrix%rhs, &
                       arp_workd(arp_pntr(2):arp_pntr(2)+neqown-1))

   case default ! ARPACK done
      exit

   end select

end do

! finish computing the average residual

if (nresid_comp /= 0) then
   grid%eigen_linsys_ave_l2_resid = grid%eigen_linsys_ave_l2_resid/nresid_comp
endif

! gather some statistics from ARPACK

grid%arpack_iter = arp_param(3)
grid%arpack_nconv = arp_param(5)
grid%arpack_numop = arp_param(9)
grid%arpack_numopb = arp_param(10)
grid%arpack_numreo = arp_param(11)
grid%arpack_info = 0

! check for failure

if (info /= 0) then
   call warning("ARPACK *aupd failed. info and nconv are",&
                intlist=(/info,grid%arpack_nconv/))
   grid%arpack_info = info
endif

! extract eigenvalue and compute eigenvector

rvec = .true. ! compute eigenvector
howmany = "A" ! compute all eigenvectors

kev = nev
if (my_real == kind(1.0e0)) then
   if (nproc > 1 .and. .not. still_sequential) then
    if (symm) then
      call psseupd(comm,rvec,howmany,select,dr,evec,neqown,lambda0, &
                  bmat,neqown,which,kev,tol,arp_resid, &
                  ncv,arnoldi_vec,neqown,arp_param,arp_pntr,arp_workd, &
                  arp_workl,lworkl,info)
    else
      call psneupd(comm,rvec,howmany,select,dr,di,evec,neqown,lambda0, &
                  0.0_my_real,workev,bmat,neqown,which,kev,tol,arp_resid, &
                  ncv,arnoldi_vec,neqown,arp_param,arp_pntr,arp_workd, &
                  arp_workl,lworkl,info)
    endif
   else
    if (symm) then
      call sseupd(rvec,howmany,select,dr,evec,neqown,lambda0, &
                  bmat,neqown,which,kev,tol,arp_resid, &
                  ncv,arnoldi_vec,neqown,arp_param,arp_pntr,arp_workd, &
                  arp_workl,lworkl,info)
    else
      call sneupd(rvec,howmany,select,dr,di,evec,neqown,lambda0, &
                  0.0_my_real,workev,bmat,neqown,which,kev,tol,arp_resid, &
                  ncv,arnoldi_vec,neqown,arp_param,arp_pntr,arp_workd, &
                  arp_workl,lworkl,info)
    endif
   endif
elseif (my_real == kind(1.0d0)) then
   if (nproc > 1 .and. .not. still_sequential) then
    if (symm) then
      call pdseupd(comm,rvec,howmany,select,dr,evec,neqown,lambda0, &
                  bmat,neqown,which,kev,tol,arp_resid, &
                  ncv,arnoldi_vec,neqown,arp_param,arp_pntr,arp_workd, &
                  arp_workl,lworkl,info)
    else
      call pdneupd(comm,rvec,howmany,select,dr,di,evec,neqown,lambda0, &
                  0.0_my_real,workev,bmat,neqown,which,kev,tol,arp_resid, &
                  ncv,arnoldi_vec,neqown,arp_param,arp_pntr,arp_workd, &
                  arp_workl,lworkl,info)
    endif
   else
    if (symm) then
      call dseupd(rvec,howmany,select,dr,evec,neqown,lambda0, &
                  bmat,neqown,which,kev,tol,arp_resid, &
                  ncv,arnoldi_vec,neqown,arp_param,arp_pntr,arp_workd, &
                  arp_workl,lworkl,info)
    else
      call dneupd(rvec,howmany,select,dr,di,evec,neqown,lambda0, &
                  0.0_my_real,workev,bmat,neqown,which,kev,tol,arp_resid, &
                  ncv,arnoldi_vec,neqown,arp_param,arp_pntr,arp_workd, &
                  arp_workl,lworkl,info)
    endif
   endif
else
   ierr = PHAML_INTERNAL_ERROR
   call fatal("my_real is neither single nor double precision. Can't call ARPACK",procs=procs)
   return
endif

! check for failure

if (info /= 0) then
   call warning("ARPACK *eupd failed.",&
                intlist=(/info/))
   grid%arpack_info = info
endif

! sort the eigenvalues

call sort(dr,nev,iperm,1,info)

! For each eigenvalue/eigenvector

do jj=1,nev
   ii = iperm(jj)

! use tmpptr to point to the solution vector in the linear system data struct

   if (jj == 1) then
      tmpptr => full_matrix%solution(1:neq_full)
   else
      tmpptr => full_matrix%evecs(:,jj-1)
   endif

! copy the eigenvectors into the solutions

   call mine_to_all(full_matrix,evec(:,ii),tmpptr,0.0_my_real)
   if (.not. still_sequential) then
      call exchange_fudop_vect(tmpptr,procs,full_matrix,2060+ii,2070+ii, &
                               2075+ii)
   endif

! find the min and max solution values

   maxsoln = -huge(0.0_my_real)
   minsoln = huge(0.0_my_real)
   do i=1,neq_full
      if (full_matrix%equation_type(i) == DIRICHLET) cycle
      if (.not. full_matrix%iown(i)) cycle
      maxsoln = max(maxsoln,tmpptr(i))
      minsoln = min(minsoln,tmpptr(i))
   end do
   if (.not. still_sequential) then
      maxsoln = phaml_global_max(procs,maxsoln,2920+ii)
      minsoln = phaml_global_min(procs,minsoln,2930+ii)
   endif

! Scale the eigenvectors so the selected norm is 1.0 and so the maximum
! absolution value is in the positive direction

   select case (solver_cntl%scale_evec)

   case (SCALE_LINF)

! Scale eigenvectors to have L infinity norm 1.0

      if (abs(minsoln) > abs(maxsoln)) then
         if (minsoln /= 0.0_my_real) tmpptr = tmpptr/minsoln
      else
         if (maxsoln /= 0.0_my_real) tmpptr = tmpptr/maxsoln
      endif

   case (SCALE_L2)

! Scale eigenvectors to have L2 norm 1.0

      soln2norm = 0.0_my_real
      do i=1,neq_full
         if (full_matrix%equation_type(i) == DIRICHLET) cycle
         if (.not. full_matrix%iown(i)) cycle
         soln2norm = soln2norm + tmpptr(i)**2
      end do
      if (.not. still_sequential) then
         soln2norm = phaml_global_sum(procs,soln2norm,2810+ii)
      endif
      soln2norm = sqrt(soln2norm)
      if (abs(minsoln) > abs(maxsoln)) then
         soln2norm = -soln2norm
      endif
      if (soln2norm /= 0.0_my_real) tmpptr = tmpptr/soln2norm

   case (SCALE_M)

! Scale eigenvectors to have M-norm 1.0

      full_matrix%matrix_val => full_matrix%mass
      call matrix_times_vector(tmpptr,Mx,full_matrix,procs,still_sequential, &
                               2850+ii,2860+ii,2865+ii,2870+ii,2880+ii, &
                               2885+ii,nocomm2=.true.)
      solnMnorm = 0.0_my_real
      do i=1,neq_full
         if (full_matrix%equation_type(i) == DIRICHLET) cycle
         if (.not. full_matrix%iown(i)) cycle
         solnMnorm = solnMnorm + tmpptr(i)*Mx(i)
      end do
      if (.not. still_sequential) then
         solnMnorm = phaml_global_sum(procs,solnMnorm,2810+ii)
      endif
      solnMnorm = sqrt(solnMnorm)
      if (abs(minsoln) > abs(maxsoln)) then
         solnMnorm = -solnMnorm
      endif
      if (solnMnorm /= 0.0_my_real) tmpptr = tmpptr/solnMnorm

   end select

! Compute Ax and Mx

   Ax = 0.0_my_real
   Mx = 0.0_my_real
   full_matrix%matrix_val => full_matrix%stiffness
   call matrix_times_vector(tmpptr,Ax,full_matrix,procs,still_sequential, &
                            2210+ii,2220+ii,2225+ii,2230+ii,2240+ii, &
                            2245+ii,nocomm2=.true.)
   full_matrix%matrix_val => full_matrix%mass
   call matrix_times_vector(tmpptr,Mx,full_matrix,procs,still_sequential, &
                            2250+ii,2260+ii,2265+ii,2270+ii,2280+ii,2285+ii, &
                            nocomm2=.true.)

! Compute the Raleigh quotient x^T A x / x^T M x to use as eigenvalues

   numer = 0.0_my_real
   denom = 0.0_my_real
   do i=1,neq_full
      if (full_matrix%equation_type(i) == DIRICHLET) cycle
      if (.not. full_matrix%iown(i)) cycle
      numer = tmpptr(i)*Ax(i) + numer
      denom = tmpptr(i)*Mx(i) + denom
   end do
   if (.not. still_sequential) then
      numer = phaml_global_sum(procs,numer,2020+ii)
      denom = phaml_global_sum(procs,denom,2030+ii)
   endif
   grid%eigenvalue(jj) = numer/denom

! If the Raleigh quotient is not close to the eigenvalue that ARPACK returned,
! it may be wrong.  Tell the master to print a warning.

   if (dr(ii) /= 0.0_my_real) then
      if (abs((dr(ii)-grid%eigenvalue(jj))/dr(ii)) > 1.0_my_real) then
         if (my_proc(procs) == 1) then
            call phaml_send(procs,MASTER,(/2,jj/),2, &
                            (/dr(ii),grid%eigenvalue(jj)/),2,2201)
         endif
      endif
   endif

! Compute the l2 norm of the residual || Ax - lambda Mx || / || lambda Mx ||

   l2resid = 0.0_my_real
   l2rhs = 0.0_my_real
   do i=1,neq_full
      if (full_matrix%equation_type(i) == DIRICHLET) cycle
      if (.not. full_matrix%iown(i)) cycle
      l2resid = l2resid + (Ax(i)-grid%eigenvalue(jj)*Mx(i))**2
      l2rhs = l2rhs + Mx(i)**2
   end do
   if (.not. still_sequential) then
      l2resid = phaml_global_sum(procs,l2resid,2040+ii)
      l2rhs = phaml_global_sum(procs,l2rhs,2050+ii)
   endif
   l2resid = sqrt(l2resid)
   l2rhs = abs(grid%eigenvalue(jj))*sqrt(l2rhs)
   if (l2rhs /= 0.0_my_real) l2resid = l2resid/l2rhs
   grid%eigenprob_l2_resid(jj) = l2resid

! the variance x^T A^T M^-1 A x / x^T M x

   numer = 0.0_my_real
   denom = 0.0_my_real
   do i=1,neq_full
      if (full_matrix%equation_type(i) == DIRICHLET) cycle
      if (.not. full_matrix%iown(i)) cycle
      numer = numer + Ax(i)**2
      denom = denom + Mx(i)**2
   end do
   if (.not. still_sequential) then
      numer = phaml_global_sum(procs,numer,2080+ii)
      denom = phaml_global_sum(procs,denom,2090+ii)
   endif
   if (denom /= 0.0_my_real) then
      grid%eigenprob_variance(jj) = numer/denom - grid%eigenvalue(jj)**2
   else
      grid%eigenprob_variance(jj) = 0.0_my_real
   endif

end do ! jj

! free workspace

deallocate(arnoldi_vec, arp_workd, arp_resid, arp_workl, dr, di, evec, &
           workev, iperm, stat=allocstat)
if (allocstat /= 0) then
   call warning("deallocation failed in eigen_arpack")
endif
if (i_made_lapack) then
   call destroy_lapack_band(linsys%lapack_mat)
   linsys%lapack_symm_band_exists = .false.
   linsys%lapack_gen_band_exists = .false.
endif
if (i_made_petsc) then
   call destroy_petsc_linear_system(linsys%petsc_matrix)
   linsys%petsc_matrix_exists = .false.
endif
if (i_made_mumps) then
   call destroy_mumps_linear_system(linsys%mumps_matrix,procs)
   linsys%mumps_matrix_exists = .false.
endif
if (i_made_superlu) then
   call destroy_superlu_linear_system(linsys%superlu_matrix,procs)
   linsys%superlu_matrix_exists = .false.
endif
if (i_made_hypre) then
   call destroy_hypre_linear_system(linsys%hypre_matrix,procs)
   linsys%hypre_matrix_exists = .false.
endif

! tell the master we're done

if (my_proc(procs) == 1) then
   call phaml_send(procs,MASTER,(/1/),1,(/0.0_my_real/),0,2201)
endif

end subroutine eigen_arpack

!        ----------
function get_matval(linear_system,matval,row,col)
!        ----------

!----------------------------------------------------
! This routine get the (row,col) matrix value in matval
!----------------------------------------------------

!----------------------------------------------------
! Dummy arguments

type(linsys_type), intent(in) :: linear_system
real(my_real), intent(in) :: matval(:)
integer, intent(in) :: row,col
real(my_real) :: get_matval
!----------------------------------------------------
! Local variables:

integer :: i

!----------------------------------------------------
! Begin executable code

! TEMP this is used to build the transpose block in static condensation.
!      If that's all it is used for, this can be made more efficient when
!      the matrix is symmetric, and maybe other times, too

! seach row for column col to get the value; return 0 if col is not found

get_matval = 0.0_my_real
do i=linear_system%begin_row(row), linear_system%end_row(row)
   if (linear_system%column_index(i) == col) then
      get_matval = matval(i)
      exit
   endif
end do

end function get_matval

!          ------------
subroutine eigen_blopex(grid,procs,linsys,io_cntl,solver_cntl,still_sequential)
!          ------------

!----------------------------------------------------
! This routine solves the eigenvalue problem Ax = lambda*Mx where A is
! stiffness and M is mass.  It uses BLOPEX.
!----------------------------------------------------

!----------------------------------------------------
! Dummy arguments

type(grid_type), intent(inout) :: grid
type(proc_info), intent(in) :: procs
type(linsys_type), intent(inout), target :: linsys
type(io_options), intent(in) :: io_cntl
type(solver_options), intent(in) :: solver_cntl
logical, intent(in) :: still_sequential

!----------------------------------------------------
! Local variables:

logical :: i_made_hypre, i_made_petsc
integer :: i, j, ii, nproc, nev, neq, neq_full, ss, comm, objtype, brank, &
           srank, objlid, allocstat, info, indx, istrt, indx0
real(my_real) :: maxsoln, minsoln, lambda0, soln2norm, solnMnorm, l2resid, &
                 l2rhs, numer, denom, temp, temp2, temp3
type(linsys_type) :: full_matrix
integer :: ni, nr, proc
integer, pointer :: irecv(:)
real(my_real), pointer :: rrecv(:)
real(my_real) :: Ax(linsys%neq_vert+linsys%neq_edge+linsys%neq_face), &
                 Mx(linsys%neq_vert+linsys%neq_edge+linsys%neq_face)
real(my_real), pointer :: tmpptr(:)
integer :: iperm(solver_cntl%num_eval)
!----------------------------------------------------
! Begin executable code

! master hangs around just to see if it needs to print a warning

if (my_proc(procs) == MASTER) then
   do
      call phaml_recv(procs,proc,irecv,ni,rrecv,nr,2201)
      select case (irecv(1))
      case (1) ! all done
         deallocate(irecv,stat=allocstat)
         exit
      case (2) ! lobpcg_solve returned an error code
         call warning("BLOPEX returned error code",intlist=(/irecv(2),irecv(3)/))
         deallocate(irecv,stat=allocstat)
         if (associated(rrecv)) deallocate(rrecv,stat=allocstat)
      end select
   end do
   return
endif

! TEMP101124 zero out Dirichlet columns -- only sequential and degree=1

do i=1,linsys%neq
   do j=linsys%begin_row(i),linsys%end_row(i)
      if (linsys%column_index(j) == NO_ENTRY) cycle
      if (linsys%equation_type(linsys%column_index(j)) == DIRICHLET .and. &
          linsys%column_index(j) /= i) then
         linsys%matrix_val(j) = 0
         linsys%stiffness(j) = 0
         linsys%mass(j) = 0
         linsys%shifted(j) = 0
         linsys%condensed(j) = 0
      endif
   end do
end do
! end TEMP101124

! use linsys for solving the (statically condensed) shifted linear systems,
! and full_matrix for operations that need the whole matrix.  full_matrix
! points to entries in linsys with the following pointers/values different:
!  linsys%end_row => end_row_edge
!  linsys%matrix_val => condensed => shifted
!  linsys%rhs => rhs_cond
!  linsys%neq = neq_vert+neq_edge

full_matrix = linsys

full_matrix%end_row => linsys%end_row_face
full_matrix%matrix_val => linsys%stiffness
full_matrix%rhs => linsys%rhs_nocond
full_matrix%neq = linsys%neq_vert + linsys%neq_edge + linsys%neq_face

! convenience variables

nproc = num_proc(procs)
comm = slaves_communicator(procs)
nev = solver_cntl%num_eval
neq = linsys%neq
neq_full = full_matrix%neq
ss = grid%system_size
if (solver_cntl%lambda0 == -huge(0.0_my_real)) then
   lambda0 = solver_cntl%lambda0
elseif (solver_cntl%lambda0_side == EIGEN_RIGHT) then
   lambda0 = -solver_cntl%lambda0
else
   lambda0 = solver_cntl%lambda0
endif

! Make sure there is space in grid for the requested number of eigenvalues

if (size(grid%eigenvalue) /= nev) then
   deallocate(grid%eigenvalue,grid%eigenprob_l2_resid, &
              grid%eigenprob_variance,stat=allocstat)
   allocate(grid%eigenvalue(nev),grid%eigenprob_l2_resid(nev), &
            grid%eigenprob_variance(nev),stat=allocstat)
   if (allocstat /= 0) then
      ierr = ALLOC_FAILED
      call fatal("memory allocation failed in eigen_blopex",procs=procs)
      return
   endif
   grid%eigenvalue = 0.0_my_real
   grid%eigenprob_l2_resid = 0.0_my_real
   grid%eigenprob_variance = 0.0_my_real
endif

! copy the solutions from the grid to the linear system to return the
! old solution if we abort

do i=1,neq_full
   call eq_to_grid(full_matrix,full_matrix%gid(i),objtype,brank,srank, &
                   objlid,grid)
   select case (objtype)
   case (VERTEX_ID)
      full_matrix%solution(i) = grid%vertex_solution(objlid,srank,1)
      if (nev > 1) then
         full_matrix%evecs(i,:) = grid%vertex_solution(objlid,srank,2:nev)
      endif
   case (EDGE_ID)
      full_matrix%solution(i) = grid%edge(objlid)%solution(brank,srank,1)
      if (nev > 1) then
         do j=1,nev-1
            full_matrix%evecs(i,j) = grid%edge(objlid)%solution(brank,srank,j)
         end do
      endif
   case (ELEMENT_ID)
      full_matrix%solution(i) = grid%element(objlid)%solution(brank,srank,1)
      if (nev > 1) then
         do j=1,nev-1
            full_matrix%evecs(i,j) = grid%element(objlid)%solution(brank,srank,j)
         end do
      endif
   end select
end do

! create any needed alternate storage formats for solution methods

i_made_hypre = .false.
i_made_petsc = .false.

select case (solver_cntl%solver)

! TEMP for later expansion
!case (HYPRE_BOOMERAMG_SOLVER, HYPRE_PCG_SOLVER, HYPRE_GMRES_SOLVER)
!   if (.not. linsys%hypre_matrix_exists) then
!      call create_hypre_linear_system(linsys%hypre_matrix,linsys,procs, &
!                                      still_sequential)
!      linsys%hypre_matrix_exists = .true.
!      i_made_hypre = .true.
!   endif

case (PETSC_RICHARDSON_SOLVER, PETSC_CHEBYCHEV_SOLVER, PETSC_CG_SOLVER, &
      PETSC_GMRES_SOLVER,      PETSC_TCQMR_SOLVER,     PETSC_BCGS_SOLVER,&
      PETSC_CGS_SOLVER,        PETSC_TFQMR_SOLVER,     PETSC_CR_SOLVER, &
      PETSC_LSQR_SOLVER,       PETSC_BICG_SOLVER)
   if (my_proc(procs) /= MASTER) then
      where (linsys%equation_type == DIRICHLET) linsys%rhs = linsys%solution(1:)
      if (.not. linsys%petsc_matrix_exists) then
         if (solver_cntl%petsc_matrix_free) then
            call create_petsc_linear_system_mf(linsys,linsys%petsc_matrix, &
                                               still_sequential)
         else
            call create_petsc_linear_system(linsys,linsys%petsc_matrix, &
                                            still_sequential,procs)
         endif
         linsys%petsc_matrix_exists = .true.
         i_made_petsc = .true.
      endif
   endif

case default
   call fatal("BLOPEX currently requires a petsc solver",procs=procs)
   stop

end select

! solve the eigenvalue problem

select case (solver_cntl%solver)

! TEMP for later expansion
!case (HYPRE_BOOMERAMG_SOLVER, HYPRE_PCG_SOLVER, HYPRE_GMRES_SOLVER)
!
!   if (my_proc(procs) /= MASTER) then
!      call hypre_lobpcg_solve_f(linsys,full_matrix,linsys%hypre_matrix, &
!                                solver_cntl,io_cntl,still_sequential,grid,procs)
!   endif

case (PETSC_RICHARDSON_SOLVER, PETSC_CHEBYCHEV_SOLVER, PETSC_CG_SOLVER, &
      PETSC_GMRES_SOLVER,      PETSC_TCQMR_SOLVER,     PETSC_BCGS_SOLVER,&
      PETSC_CGS_SOLVER,        PETSC_TFQMR_SOLVER,     PETSC_CR_SOLVER, &
      PETSC_LSQR_SOLVER,       PETSC_BICG_SOLVER)
   if (my_proc(procs) /= MASTER) then
      call petsc_lobpcg_solve_f(linsys,full_matrix,linsys%petsc_matrix, &
                                solver_cntl,io_cntl,still_sequential,grid,procs)
   endif

case default
   call fatal("BLOPEX currently requires a petsc solver",procs=procs)
   stop

end select

! For each eigenvalue/eigenvector

do ii=1,nev

! use tmpptr to point to the solution vector in the linear system data struct

   if (ii == 1) then
      tmpptr => full_matrix%solution(1:neq_full)
   else
      tmpptr => full_matrix%evecs(:,ii-1)
   endif

! find the min and max solution values

   maxsoln = -huge(0.0_my_real)
   minsoln = huge(0.0_my_real)
   do i=1,neq_full
      if (full_matrix%equation_type(i) == DIRICHLET) cycle
      if (.not. full_matrix%iown(i)) cycle
      maxsoln = max(maxsoln,tmpptr(i))
      minsoln = min(minsoln,tmpptr(i))
   end do
   if (.not. still_sequential) then
      maxsoln = phaml_global_max(procs,maxsoln,2920+ii)
      minsoln = phaml_global_min(procs,minsoln,2930+ii)
   endif

! Scale the eigenvectors so the selected norm is 1.0 and so the maximum
! absolution value is in the positive direction

   select case (solver_cntl%scale_evec)

   case (SCALE_LINF)

! Scale eigenvectors to have L infinity norm 1.0

      if (abs(minsoln) > abs(maxsoln)) then
         if (minsoln /= 0.0_my_real) tmpptr = tmpptr/minsoln
      else
         if (maxsoln /= 0.0_my_real) tmpptr = tmpptr/maxsoln
      endif

   case (SCALE_L2)

! Scale eigenvectors to have L2 norm 1.0

      soln2norm = 0.0_my_real
      do i=1,neq_full
         if (full_matrix%equation_type(i) == DIRICHLET) cycle
         if (.not. full_matrix%iown(i)) cycle
         soln2norm = soln2norm + tmpptr(i)**2
      end do
      if (.not. still_sequential) then
         soln2norm = phaml_global_sum(procs,soln2norm,2810+ii)
      endif
      soln2norm = sqrt(soln2norm)
      if (abs(minsoln) > abs(maxsoln)) then
         soln2norm = -soln2norm
      endif
      if (soln2norm /= 0.0_my_real) tmpptr = tmpptr/soln2norm

   case (SCALE_M)

! Scale eigenvectors to have M-norm 1.0

      full_matrix%matrix_val => full_matrix%mass
      call matrix_times_vector(tmpptr,Mx,full_matrix,procs,still_sequential, &
                               2850+ii,2860+ii,2865+ii,2870+ii,2880+ii, &
                               2885+ii,nocomm2=.true.)
      solnMnorm = 0.0_my_real
      do i=1,neq_full
         if (full_matrix%equation_type(i) == DIRICHLET) cycle
         if (.not. full_matrix%iown(i)) cycle
         solnMnorm = solnMnorm + tmpptr(i)*Mx(i)
      end do
      if (.not. still_sequential) then
         solnMnorm = phaml_global_sum(procs,solnMnorm,2810+ii)
      endif
      solnMnorm = sqrt(solnMnorm)
      if (abs(minsoln) > abs(maxsoln)) then
         solnMnorm = -solnMnorm
      endif
      if (solnMnorm /= 0.0_my_real) tmpptr = tmpptr/solnMnorm

   end select

! Compute Ax and Mx

   Ax = 0.0_my_real
   Mx = 0.0_my_real
   full_matrix%matrix_val => full_matrix%stiffness
   call matrix_times_vector(tmpptr,Ax,full_matrix,procs,still_sequential, &
                            2210+ii,2220+ii,2225+ii,2230+ii,2240+ii, &
                            2245+ii,nodirch=.true.,nocomm2=.true.)
   if (lambda0 /= -huge(0.0_my_real) .and. &
       solver_cntl%lambda0_side == EIGEN_RIGHT) then
      Ax = -Ax
   endif
   full_matrix%matrix_val => full_matrix%mass
   call matrix_times_vector(tmpptr,Mx,full_matrix,procs,still_sequential, &
                            2250+ii,2260+ii,2265+ii,2270+ii,2280+ii,2285+ii, &
                            nodirch=.true.,nocomm2=.true.)

! Compute the Raleigh quotient x^T A x / x^T M x to use as eigenvalues
! This is necessary for SHIFT_SQUARE because we don't know which square root
! to use, but we could use the transformed eigenvalue for SHIFT_INVERT.

   numer = 0.0_my_real
   denom = 0.0_my_real
   do i=1,neq_full
      if (full_matrix%equation_type(i) == DIRICHLET) cycle
      if (.not. full_matrix%iown(i)) cycle
      numer = tmpptr(i)*Ax(i) + numer
      denom = tmpptr(i)*Mx(i) + denom
   end do
   if (.not. still_sequential) then
      numer = phaml_global_sum(procs,numer,2020+ii)
      denom = phaml_global_sum(procs,denom,2030+ii)
   endif
   grid%eigenvalue(ii) = numer/denom

! Compute the l2 norm of the residual || Ax - lambda Mx || / || lambda Mx ||

   l2resid = 0.0_my_real
   l2rhs = 0.0_my_real
   do i=1,neq_full
      if (full_matrix%equation_type(i) == DIRICHLET) cycle
      if (.not. full_matrix%iown(i)) cycle
      l2resid = l2resid + (Ax(i)-grid%eigenvalue(ii)*Mx(i))**2
      call eq_to_grid(full_matrix,full_matrix%gid(i),objtype,brank,srank, &
                      objlid,grid)
      l2rhs = l2rhs + Mx(i)**2
   end do
   if (.not. still_sequential) then
      l2resid = phaml_global_sum(procs,l2resid,2040+ii)
      l2rhs = phaml_global_sum(procs,l2rhs,2050+ii)
   endif
   l2resid = sqrt(l2resid)
   l2rhs = abs(grid%eigenvalue(ii))*sqrt(l2rhs)
   if (l2rhs /= 0.0_my_real) l2resid = l2resid/l2rhs
   grid%eigenprob_l2_resid(ii) = l2resid

! the variance x^T A^T M^-1 A x / x^T M x

   numer = 0.0_my_real
   denom = 0.0_my_real
   do i=1,neq_full
      if (full_matrix%equation_type(i) == DIRICHLET) cycle
      if (.not. full_matrix%iown(i)) cycle
      numer = numer + Ax(i)**2
      denom = denom + Mx(i)**2
   end do
   if (.not. still_sequential) then
      numer = phaml_global_sum(procs,numer,2080+ii)
      denom = phaml_global_sum(procs,denom,2090+ii)
   endif
   if (denom /= 0.0_my_real) then
      grid%eigenprob_variance(ii) = numer/denom - grid%eigenvalue(ii)**2
   else
      grid%eigenprob_variance(ii) = 0.0_my_real
   endif

end do ! ii

! sort the eigenvalues, if they are on the left or both sides of lambda0

if (nev > 1 .and. lambda0 /= -huge(0.0_my_real)) then
   allocate(tmpptr(neq_full),stat=allocstat)
   if (allocstat /= 0) then
      ierr = ALLOC_FAILED
      call fatal("memory allocation failed in eigen_blopex",procs=procs)
      return
   endif

   call sort(grid%eigenvalue,nev,iperm,1,info)

   do i=1,nev
      indx=abs(iperm(i))
      if (iperm(indx) > 0) iperm(indx) = -iperm(indx)
   end do
   do istrt = 1,nev
      if (iperm(istrt) > 0) cycle
      indx = istrt
      indx0 = indx
      temp = grid%eigenvalue(istrt)
      temp2 = grid%eigenprob_l2_resid(istrt)
      temp3 = grid%eigenprob_variance(istrt)
      if (istrt == 1) then
         tmpptr = full_matrix%solution(1:)
      else
         tmpptr = full_matrix%evecs(:,istrt-1)
      endif
      do while (iperm(indx) < 0)
         grid%eigenvalue(indx) = grid%eigenvalue(-iperm(indx))
         grid%eigenprob_l2_resid(indx) = grid%eigenprob_l2_resid(-iperm(indx))
         grid%eigenprob_variance(indx) = grid%eigenprob_variance(-iperm(indx))
         if (indx == 1) then
            full_matrix%solution(1:) = full_matrix%evecs(:,-iperm(indx)-1)
         elseif (-iperm(indx) == 1) then
            full_matrix%evecs(:,indx-1) = full_matrix%solution(1:)
         else
            full_matrix%evecs(:,indx-1) = full_matrix%evecs(:,-iperm(indx)-1)
         endif
         indx0 = indx
         iperm(indx) = -iperm(indx)
         indx = iperm(indx)
      end do
      grid%eigenvalue(indx0) = temp
      grid%eigenprob_l2_resid(indx0) = temp2
      grid%eigenprob_variance(indx0) = temp3
      if (indx0 == 1) then
         full_matrix%solution(1:) = tmpptr
      else
         full_matrix%evecs(:,indx0-1) = tmpptr
      endif
   end do

   deallocate(tmpptr)
endif

! get rid of any extra versions of the linear system

if (i_made_hypre) then
   call destroy_hypre_linear_system(linsys%hypre_matrix,procs)
   linsys%hypre_matrix_exists = .false.
endif
if (i_made_petsc) then
   call destroy_petsc_linear_system(linsys%petsc_matrix)
   linsys%petsc_matrix_exists = .false.
endif

! tell the master we're done

if (my_proc(procs) == 1) then
   call phaml_send(procs,MASTER,(/1/),1,(/0.0_my_real/),0,2201)
endif

end subroutine eigen_blopex

end module eigen

max_stars_repo_name: qsnake/phaml

max_issues_repo_name: qsnake/phaml

max_forks_repo_name: qsnake/phaml


text:
C$Procedure GETFNM ( Get a filename from standard input )

      SUBROUTINE GETFNM ( PRMPT, FSTAT, FNAME, VALID, MESSG )
      IMPLICIT NONE

C$ Abstract
C
C     This routine prompts the user for a valid filename.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     None.
C
C$ Keywords
C
C     UTILITY
C
C$ Declarations

      EXTERNAL              PROMPT

      CHARACTER*(*)         PRMPT
      CHARACTER*(*)         FSTAT
      CHARACTER*(*)         FNAME
      LOGICAL               VALID
      CHARACTER*(*)         MESSG

C$ Brief_I/O
C
C     Variable  I/O  Description
C     --------  ---  --------------------------------------------------
C     PRMPT      I   The prompt to use when asking for the filename.
C     FSTAT      I   Status of the file: 'OLD' or 'NEW'.
C     FNAME      O   A valid filename typed in by the user.
C     VALID      O   A logical flag indicating a valid filename.
C     MESSG      O   A descriptive message for an invalid filename.
C
C$ Detailed_Input
C
C     PRMPT    is a character string that will be displayed from the
C              active position of the cursor to the end of string
C              that lets a user know that input is expected.
C
C     FSTAT    This is the status of the filename entered. It should
C              be 'OLD' when prompting for the filename of a file which
C              already exists, and 'NEW' when prompting for the filename
C              of a file which does not already exist or is to be over
C              written.
C
C$ Detailed_Output
C
C     FNAME    is a character string that contains a valid filename
C              typed in by the user. A valid filename is defined
C              simply to be a nonblank character string with no
C              embedded blanks, nonprinting characters, or characters
C              having decimal values > 126.
C
C     VALID    A logical flag which indicates whether or not the
C              filename entered is valid, i.e., a nonblank character
C              string with no leading or embedded blanks, which
C              satisfies the constraints for validity imposed.
C
C     MESSG    A brief descriptive message which describes why a
C              particular filename was not valid. Blank if a valid
C              filename is entered.
C
C$ Parameters
C
C     None.
C
C$ Exceptions
C
C     None.
C
C$ Files
C
C     None.
C
C$ Particulars
C
C     This is a utility that allows you to "easily" request a valid,
C     filename from a program user.  At a high level, it frees you
C     from the peculiarities of a particular FORTRAN's implementation
C     of cursor control.
C
C     A valid filename is defined as a nonblank character string with
C     no embedded blanks, nonprinting characters, or characters with
C     decimal values > 126. Leading blanks are removed, and trailing
C     blanks are ignored.
C
C     If an invalid filename is entered, this routine provides a
C     descriptive error message and halts the execution of the
C     process which called it by using a Fortran STOP.
C
C$ Examples
C
C     EXAMPLE 1:
C
C        FNAME = ' '
C        PRMPT = 'Filename? '
C        FSTAT = 'OLD'
C
C        CALL GETFNM( PRMPT, FSTAT, FNAME, VALID, MESSG )
C
C     The user sees the following displayed on his screen:
C
C        Filename? _
C
C     where the underbar, '_', represents the cursor position.
C
C$ Restrictions
C
C     None.
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     K.R. Gehringer (JPL)
C
C$ Version
C
C-    SPICELIB Version 5.17.0, 10-MAR-2014 (BVS)
C
C        Updated for SUN-SOLARIS-64BIT-INTEL.
C
C-    SPICELIB Version 5.16.0, 10-MAR-2014 (BVS)
C
C        Updated for PC-LINUX-64BIT-IFORT.
C
C-    SPICELIB Version 5.15.0, 10-MAR-2014 (BVS)
C
C        Updated for PC-CYGWIN-GFORTRAN.
C
C-    SPICELIB Version 5.14.0, 10-MAR-2014 (BVS)
C
C        Updated for PC-CYGWIN-64BIT-GFORTRAN.
C
C-    SPICELIB Version 5.13.0, 14-DEC-2010 (EDW)
C
C        Declared PROMPT as EXTERNAL.
C
C        Unfied Version and Revision sections, eliminated Revision
C        section. Corrected error in 09-DEC-1999 Version entry.
C        Version ID changed to 5.0.9 from 7.0.0.
C
C-    Beta Version 5.12.0, 13-MAY-2010 (BVS)
C
C        Updated for SUN-SOLARIS-INTEL.
C
C-    Beta Version 5.11.0, 13-MAY-2010 (BVS)
C
C        Updated for SUN-SOLARIS-64BIT-NATIVE_C.
C
C-    Beta Version 5.10.0, 13-MAY-2010 (BVS)
C
C        Updated for PC-LINUX-64BIT-GFORTRAN.
C
C-    Beta Version 5.9.0, 13-MAY-2010 (BVS)
C
C        Updated for MAC-OSX-64BIT-IFORT.
C
C-    Beta Version 5.8.0, 13-MAY-2010 (BVS)
C
C        Updated for MAC-OSX-64BIT-GFORTRAN.
C
C-    Beta Version 5.7.0, 18-MAR-2009 (BVS)
C
C        Updated for PC-LINUX-GFORTRAN.
C
C-    Beta Version 5.6.0, 18-MAR-2009 (BVS)
C
C        Updated for MAC-OSX-GFORTRAN.
C
C-    Beta Version 5.5.0, 19-FEB-2008 (BVS)
C
C        Updated for PC-LINUX-IFORT.
C
C-    Beta Version 5.4.0, 14-NOV-2006 (BVS)
C
C        Updated for MAC-OSX-IFORT.
C
C-    Beta Version 5.3.0, 26-OCT-2005 (BVS)
C
C        Updated for SUN-SOLARIS-64BIT-GCC_C.
C
C-    Beta Version 5.2.0, 03-JAN-2005 (BVS)
C
C        Updated for PC-CYGWIN.
C
C-    Beta Version 5.1.1, 17-JUL-2002 (BVS)
C
C        Added MAC-OSX environments.
C
C-    Beta Version 5.1.0, 16-AUG-2000 (WLT)
C
C        Added the PC-LINUX environment
C
C-    Beta Version 5.0.9, 09-DEC-1999 (WLT)
C
C        This routine now calls EXPFNM_2 only in UNIX environments
C
C-    Beta Version 5.0.0, 20-JAN-1998 (NJB)
C
C        Now calls EXPFNM_2 to attempt to expand environment variables.
C
C        Fixed a typo or two at various places in the header.
C
C-    Beta Version 4.0.1, 25-APR-1994 (KRG)
C
C        Removed some incorrect comments from the $ Particulars section
C        of the header. Something about a looping structure that is not
C        a part of the code now, if it ever was.
C
C        Fixed a typo or two at various places in the header.
C
C-    Beta Version 4.0.0, 29-SEP-1993 (KRG)
C
C        Added the character reperesnted by decimal 127 to the BADCHR.
C        It should have been there, but it wasn't.
C
C-    Beta Version 3.0.0, 10-SEP-1993 (KRG)
C
C        Made the file status variable FSTAT case insensitive.
C
C        Added code to the  file status .EQ. 'NEW' case to set the
C        valid flag to .FALSE. and set an appropriate error message
C        about the file already existing.
C
C-    Beta Version 2.0.0, 02-APR-1993 (KRG)
C
C        The variable BADCHR was not saved which caused problems on some
C        computers. This variable is now saved.
C
C-    Beta Version 1.0.0, 01-JUN-1992 (KRG)
C
C-&

C$ Index_Entries
C
C      prompt for a filename with error handling
C
C-&

C
C     SPICELIB Functions
C
      INTEGER               CPOS
      INTEGER               LASTNB
      INTEGER               RTRIM

      LOGICAL               EXISTS
      LOGICAL               RETURN
      LOGICAL               FAILED

      INTEGER               MAXNAM
      PARAMETER           ( MAXNAM = 1000 )

      CHARACTER*(MAXNAM)    NAMBUF

C
C     Local Parameters
C
      CHARACTER*(*)         DPRMPT
      PARAMETER           ( DPRMPT = 'Filename? ' )

      CHARACTER*(*)         OSTAT
      PARAMETER           ( OSTAT  = 'OLD' )

      CHARACTER*(*)         NSTAT
      PARAMETER           ( NSTAT  = 'NEW' )

C
C     Local Variables
C
      CHARACTER*(162)       BADCHR
      CHARACTER*(3)         STATUS

      INTEGER               I
      INTEGER               LENGTH

      LOGICAL               FIRST
C
C     Saved Variables
C
      SAVE BADCHR
      SAVE FIRST
C
C     Initial Values
C
      DATA                  FIRST   / .TRUE. /
C
C     Standard SPICE error handling.
C
      IF ( RETURN () ) THEN
         RETURN
      ELSE
         CALL CHKIN ( 'GETFNM' )
      END IF
C
C     If this is the first time this routine has been called, initialize
C     the ``bad character'' string.
C
      IF ( FIRST ) THEN

         FIRST = .FALSE.

         DO I = 0, 32
            BADCHR(I+1:I+1) = CHAR( I )
         END DO

         DO I = 1, 129
            BADCHR(33+I:33+I) = CHAR( 126+I )
         END DO

      END IF
C
C     Set the value of the valid flag to .TRUE.. We might as well assume
C     that the name entered will be a valid one.
C
      VALID = .TRUE.
C
C     Left justify and convert the file status to upper case for
C     comparisons.
C
      CALL LJUST ( FSTAT,  STATUS )
      CALL UCASE ( STATUS, STATUS )
C
C     Check to see if we have a valid status for the filename.
C
      IF ( ( STATUS .NE. OSTAT ) .AND.
     .     ( STATUS .NE. NSTAT ) ) THEN

         VALID = .FALSE.
         MESSG = 'The status '''//STATUS//''' was not recognized.'
         CALL CHKOUT ( 'GETFNM' )
         RETURN

      END IF
C
C     Read in a potential filename, and test it for validity.
C
      IF ( PRMPT .EQ. ' ' ) THEN

         CALL PROMPT( DPRMPT, FNAME )

      ELSE

         CALL PROMPT ( PRMPT, FNAME )

      END IF

C
C     The string we just obtained could be an environment variable.
C     If so, expand it.
C
      CALL EXPFNM_2 ( FNAME, NAMBUF )

      IF ( FAILED() ) THEN
         CALL CHKOUT ( 'GETFNM' )
         RETURN
      END IF

      FNAME = NAMBUF

      IF ( FNAME .EQ. ' ' ) THEN

         VALID = .FALSE.
         MESSG = 'A blank filename is not valid.'
         CALL CHKOUT ( 'GETFNM' )
         RETURN

      END IF

C
C     Left justify the filename.
C
      CALL LJUST ( FNAME, FNAME )

C
C     Check for bad characters in the filename.
C
      LENGTH = LASTNB( FNAME )
      I      = CPOS ( FNAME(:LENGTH), BADCHR, 1 )

      IF ( I .GT. 0 ) THEN

         VALID = .FALSE.
         MESSG = 'Invalid filename. Illegal character encountered:' //
     .           ' decimal value: #'
         CALL REPMI  ( MESSG, '#', ICHAR(FNAME(I:I)), MESSG )
         CALL CHKOUT ( 'GETFNM' )
         RETURN

      END IF

C
C     We know that the filename that was entered was nonblank and had
C     no bad characters. So, now we take care of the status question.
C
      IF ( STATUS .EQ. OSTAT ) THEN

         IF ( .NOT. EXISTS( FNAME(1:RTRIM(FNAME))  ) ) THEN

            VALID = .FALSE.
            MESSG = 'The file does not exist.'
            CALL CHKOUT ( 'GETFNM' )
            RETURN

         END IF

      ELSE IF ( STATUS .EQ. NSTAT ) THEN

         IF ( EXISTS( FNAME(1:RTRIM(FNAME))  ) ) THEN

            VALID = .FALSE.
            MESSG = 'The file already exists.'
            CALL CHKOUT ( 'GETFNM' )
            RETURN

         END IF

      END IF

C
C     At this point, we have done the best we can. If the status
C     was new, we might still have an invalid filename, but the
C     exact reasons for its invalidity are system dependent.
C
      CALL CHKOUT ( 'GETFNM' )

      RETURN
      END

ext: f

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
C
C
*+ plot_getuv2

       subroutine plot_getuv2 (minirt, prompt, default, iuv, status)
C      -------------------------------------------------------------
C
C  Prompts for a map U,V coordinate window.
C
C  Given:
C    mini redtape of map
       integer          minirt(8)
C    prompt
       character*(*)    prompt
C    default
       character*(*)    default
C Returned:
C    UV-coordinate window
       integer          iuv(4)
C
C Updated:
C   error status
       integer          status
C
C  Prompts on the input device for a map area, specified as four integers,
C  (IU1,IU2,IV1,IV2).  If only one integer is entered then the square area
C  (-IU1,IU1,IU1,-IU1) is returned.  Note that the ordering convention
C   IU1<IU2 and IV1>IV2  is enforced if necessary.  If the redtape common
C  blocks contain parameters for a current map, then checks are also made
C  that the area lies within the map bounds.  If the input string is null,
C  the DEFAULT string is used as input.  If the default string is null or
C  '*', the current values of IUV are used as default.
C
C  The user may optionally specify "CURSOR" in response to the prompt and
C  if then a map is plotted on the current plot device (and open) then
C  the UV-range may be specified by positioning the cursor on the graphics
C  screen in the usual manner.
C
C  The STATUS value should be zero on entry:  the returned status value
C  is zero if successful.  Other possible values are:
C
C      - invalid u,v window (ILL_UVWIND)
C      - u,v window outside map (UV_OUTMAP)
C      - unexpected system I/O error, or user break
C
C  (DJT, 14 October 86)
C  (PA,  22 January 88)
C  (PA,  17 December 91)
C  (PA,  17 January 93)
*-
       character*20   string
       integer        iuv_n(4), iu1, iu2, iv1, iv2
       integer        iu, iv, ls, errdev
       real*4         uv_c1(2), uv_c2(2)
       logical        batch, chr_cmatch, verify_value, cursor_available,
     *                cmd_dblev

       include '/mrao/include/iolib_errors.inc'
       include '/mrao/include/maplib_errors.inc'
       include '../include/plt_basic_defn.inc'
       include '../include/plt_image_defn.inc'

        equivalence (iuv_n(1),iu1),(iuv_n(2),iu2),(iuv_n(3),iv1),
     *              (iuv_n(4),iv2)

    1  if (status.eq.0) then

         write(string,'(3(i4,'',''),i4)')iuv

C    Read U,V window from command line

         call io_getwrd(prompt,default,string,ls,status)
         if (status.ne.0) then
           call cmd_err(status,'PLOT_GETUV',' ')
           return
         end if

C Look for "CURSOR" response
         call chr_chucas(string)
         if (chr_cmatch(string(1:ls),'CURSOR')) then

           call pgqcur(cursor_available)
           if (frame_init .and. cursor_available) then

C ... look for cursor input
             verify_value = .true.
             if (cmd_dblev(1)) then
               call io_wrout(
     *         '.. mark the UV-range with the graphics cursor')
               call io_wrout(
     *         '   mark opposite corners of the window')
               call io_wrout(
     *         '   type any character to enter or a Q to quit')
             endif
             call plot_cursor_get( uv_c1, status )
             call plot_cursor_get( uv_c2, status )
             iuv_n(1) = min(nint(uv_c1(1)),nint(uv_c2(1)))
             iuv_n(2) = max(nint(uv_c1(1)),nint(uv_c2(1)))
             iuv_n(3) = max(nint(uv_c1(2)),nint(uv_c2(2)))
             iuv_n(4) = min(nint(uv_c1(2)),nint(uv_c2(2)))

          else

C ... no plot on output device therefore do not use cursor input
             call cmd_wrerr('UV-RANGE',
     *                       'NO plot frame or NO cursor available')
             status = ill_uvwind

           end if

         else if (ls.gt.0) then

C .. normal terminal prompt
           verify_value = .true.
           call chr_chctoi(string(1:ls),iu1,status)
           call io_nxtwrd(string,ls,status)
           if (ls.gt.0) then
             call chr_chctoi(string(1:ls),iu2,status)
             call io_nxtwrd(string,ls,status)
             call chr_chctoi(string(1:ls),iv1,status)
             call io_nxtwrd(string,ls,status)
             call chr_chctoi(string(1:ls),iv2,status)
           else
             iu1=-iu1
             iu2=-iu1
             iv1=iu2
             iv2=iu1
             ls=1
           endif
         end if
C
C    Check for validity
C
         if (verify_value) then
           if (ls.eq.0 .or. status.eq.bad_syntax) then
             call io_wrout('*** bad syntax,  U1,U2,V1,V2  wanted')
             status=ill_uvwind
           else
C
C    Enforce IU1 < IU2 and IV1 > IV2
C
             if (iu1.gt.iu2) then
               iu=iu1
               iu1=iu2
               iu2=iu
             endif
             if (iv1.lt.iv2) then
               iv=iv1
               iv1=iv2
               iv2=iv
             endif
C
C    Check against map bounds
               if (iu1.lt.minirt(1) .or. iu2.gt.minirt(2) .or.
     :             iv1.gt.minirt(3) .or. iv2.lt.minirt(4) ) then
                 call io_enqerr( errdev )
                 write(errdev,'(x,a,4i6)')
     :             '*** area outside map, restrict to :',
     :             minirt(1),minirt(2),minirt(3),minirt(4)
                 status=uv_outmap
               endif
           endif

           if (status.eq.ill_uvwind .or. status.eq.uv_outmap) then
             call io_enqbch( batch )
             if (.not.batch) status=0
             call io_setcli(' ')
             goto 1
           endif

           iuv(1)=iu1
           iuv(2)=iu2
           iuv(3)=iv1
           iuv(4)=iv2

         endif
       endif

       end

max_stars_repo_name: CavendishAstrophysics/anmap

max_issues_repo_name: CavendishAstrophysics/anmap

max_forks_repo_name: CavendishAstrophysics/anmap