      subroutine subsp(b,v,g,h,d,dp,dtol,p,z,nf,nev,neq,
     &                 shift,tol,prt,its)

c     Purpose: Subspace iteration to extract lowest nf eigenpairs

c     Inputs:
c        v(*)     - Diagonal "mass" matrix
c        g(*)     - Subspace projection of stiffness matrix storage
c        h(*)     - Subspace projection of "mass" matrix storage
c        dp(*)    - Storage for previous eigenvalues
c        dtol(*)  - Storage for tolerance of computed eigenvalues
c        p(nev,*) - Storage for subspace eigenvectors
c        z(neq,*) - Storage for iteration vectors
c        nf       - Number of eigenvalues to converge
c        nev      - Size of subspace problem
c        neq      - Size of finite element equations
c        shift    - Value of shift 
c        tol      - Solution tolerance on eigenvalues
c        prt      - Flag, output each iteration result
c        its      - Maximum number of subspace iterations to perform

c     Outputs:
c        v(neq,*) - Eigenvectors
c        d(*)     - Eigenvalaues

      implicit none

      logical  conv,prt
      integer  nf,nev,neq,its
      integer  i,it,itt,itlim, j, k, n,nmas

      real*8   b(*),v(neq,*),g(*),h(*),d(*),dp(*)
      real*8   dtol(*),p(nev,*),z(neq,*)
      real*8   dm,shift,tol,tolmx,told

      integer         iodr,iodw,ipd,ipr,ipi
      common /iofild/ iodr,iodw,ipd,ipr,ipi

      integer         ior,iow
      common /iofile/ ior,iow

c     Compute initial iteration vectors

      call pconsd(v,nev*neq,0.0d0)
      nmas = 0

c     Count number of nonzero masses

      do 100 n = 1,neq
        if(b(n).ne.0.0d0) nmas = nmas + 1
100   continue

      nmas = nmas/nev
      i = 0
      j = 1
      do 110 n = 1,neq
        dm = b(n)
        if(dm.ne.0.0d0) then
          v(n,j) = dm
          i = i + 1
          if(mod(i,nmas).eq.0) j = j + 1
          j = min(j,nev)
        endif
110   continue
      do 120 i = 1,nev
        dp(i)   = 0.0d0
        dtol(i) = 1.0d0
        call scalev(v(1,i),neq)
120   continue

c     Compute new vectors and project 'a' onto 'g'

      told  = tol
      conv  = .false.
      itlim = its
      if(nev.eq.nf) itlim = 1
      do 300 it = 1,itlim
        itt = it

c       Project 'b' matrix to form 'h' and compute 'z' vectors

        call sprojb(b,v,z,h,neq,nev)

c       Project 'a' matrix to form 'g'

        call sproja(v,z,g,neq,nev)

c       Solve reduced eigenproblem 'g*p = h*p*d'

        call geig(g,h,d,p,v,nev,prt)

c       Check for convergence

        tolmx = 0.0d0
        do 200 n = 1,nev
          if(d(n).ne.0.0d0) dtol(n) = abs((d(n)-dp(n))/d(n))
          dp(n) = d(n)
          if(n.le.nf) tolmx = max(tolmx,dtol(n))
200     continue

        if(prt) then
          write(iow,2000) it,(d(n),n=1,nev)
          if(ior.lt.0) write(*,2000) it,(d(n),n=1,nev)
          if(itlim.gt.1) write(iow,2003) (dtol(n),n=1,nev)
          if(ior.lt.0.and.itlim.gt.1)
     &      write(*,2003) (dtol(n),n=1,nev)
        else
          if(ior.lt.0) write(*,2004) it,tolmx
        endif

c       Tolerance check

        do 210 n = 1,nf
          if(dtol(n).gt.told) go to 220
210     continue
        conv = .true.

c       Divide eigenvectors by eigenvalue to prevent overflows

220     do 235 i = 1,nev
          dm = d(i)
          if(p(i,i).lt.-0.00001d0) dm = -dm
          do 230 j = 1,nev
            p(j,i) = p(j,i)/dm
230       continue
235     continue

c       Compute new iteration vector 'u' from 'z'

        do 250 i = 1,neq
        do 250 j = 1,nev
          v(i,j) = 0.0d0
          do 240 k = 1,nev
            v(i,j) = v(i,j) + z(i,k)*p(k,j)
240       continue
250     continue
        if(conv) go to 305
300   continue

c     Scale vectors to have maximum element of 1.0

305   do 310 n = 1,nev
        d(n)  = 1.0/d(n) + shift
        dp(n) = sqrt(abs(d(n)))
        call scalev(v(1,n),neq)
310   continue

      write(iow,2001) itt,(d(n),n=1,nev)
      write(iow,2002)    (dp(n),n=1,nev)
      if(itt.gt.1) write(iow,2003) (dtol(n),n=1,nev)
      if(ior.lt.0) then
        write(*,2001) itt,(d(n),n=1,nev)
        write(*,2002)    (dp(n),n=1,nev)
        if(itt.gt.1) write(*,2003) (dtol(n),n=1,nev)
      endif

2000  format(/5x,'Current reciprocal shifted eigenvalues, iteration',
     &        i4/(4d20.8))
2001  format(/5x,'Solution  for  eigenvalues, iteration',i4/(4d20.8))
2002  format( 5x,'Square root of eigenvalues'/(4d20.8))
2003  format( 5x,'Current residuals'/(4d20.8))
2004  format( '+  Iteration',i3,' Max tol =',1p1e11.4)

      end

      subroutine chlbac(u,s,nn)

c     Purpose: Compute eigenvalues of general linear problem

c     Inputs:
c        u(*)     - Cholesky factor of "mass" subspace matrix
c        s(nn,*)  - Eigenvectors of standard problems
c        nn       - Size of subspace probleml

c     Outputs:
c        s(nn,*)  - Eigenvectors of generalized problem

      implicit none
      integer  i,j,jd,nn
      real*8 u(*),s(nn,nn)

c     Compute eigenvalues of general linear problem by backsubstitution

      j  = nn
      jd = nn*(nn+1)/2
      do 100 i = 1,nn
        s(nn,i) = s(nn,i)/u(jd)
100   continue

200   jd = jd - j
      j  = j - 1
      if(j.le.0) return
      do 300 i = 1,nn
        call colbac(u(jd+1),s(1,i),u(jd),j)
300   continue
      go to 200

      end

      subroutine chlfwd(u,g,s,nn)

c     Purpose: Forward solution using Choleski factors to project
c              onto a standard eigenproblem

c     Inputs:
c        u(*)     - Cholesky factor of "mass" subspace matrix
c        g(*)     - Subspace stiffness matrix
c        s(nn,*)  - Working storage

c     Outputs:
c        g(*)     - Standard eigenproblem matrix

      implicit none
      integer  i,id,im, j,jd, nn
      real*8   u(*),g(*),s(nn,nn),dot

      s(1,1) = g(1)/u(1)
      if(nn.eq.1) go to 300
      id = 1
      do 200 i = 2,nn
        s(1,i) = g(id+1)/u(1)
        im     = i - 1
        jd     = 0
        do 100 j = 1,im
          s(i,j)            = (g(id+j) -dot(u(id+1),s(1,j),im))/u(id+i)
          if(j.gt.1) s(j,i) = (g(id+j) -dot(u(jd+1),s(1,i),j-1))/u(jd+j)
          jd                = jd + j
100     continue
        id     = id + i
        s(i,i) = (g(id) - dot(u(id-im),s(1,i),im))/u(id)
200   continue

c     Complete projection

300   g(1) = s(1,1)/u(1)
      if(nn.eq.1) return
      jd = 2
      do 500 j = 2,nn
        g(jd) = s(j,1)/u(1)
        id    = 2
        do 400 i = 2,j
          im       = i - 1
          g(jd+im) = (s(j,i) - dot(u(id),g(jd),im))/u(id+im)
          id       = id + i
400     continue
        jd = jd + j
500   continue

      end

      subroutine chlfac(a,nn)

c     Purpose: Choleski factorization of a symmetric, positive definite matrix

c     Inputs
c        a(*)     - Matrix to factor
c        nn       - Size of matrix

c     Outputs
c        a(*)     - Factored matrix

      implicit none
      integer  i,id, j,jd,jm, nn
      real*8   a(*),dot

      a(1) = sqrt(a(1))
      if(nn.eq.1) return
      jd = 1
      do 200 j = 2,nn
        jm = j - 1
        id = 0
        do 100 i = 1,jm
          if(i-1.gt.0) a(jd+i) = a(jd+i) - dot(a(id+1),a(jd+1),i-1)
          id      = id + i
          a(jd+i) = a(jd+i)/a(id)
100     continue
        a(jd+j) = sqrt(a(jd+j) - dot(a(jd+1),a(jd+1),jm))
        jd      = jd + j
200   continue

      end

      subroutine colbac(u,s,d,jj)

c     Purpose: Projection to standard eigenproblem Backsubstitution macro

c     Inputs:
c        u(*)     - Upper triangular factor of matrix
c        d        - Diagonal to divide into solution
c        jj       - Length of column

c     Outputs:
c        s(*)     - Back substituion solution

      implicit none
      integer  j,jj
      real*8   u(*),s(*),d,dd

      dd = s(jj+1)
      do 100 j = 1,jj
        s(j) = s(j) - dd*u(j)
100   continue
      s(jj) = s(jj)/d

      end

      subroutine eisql(a,d,e,z,n,ierr)

c     Purpose: Solve standard eigen problem

c     Inputs:
c        a(*)     - Matrix of standard problem
c        e(*)     - Working vector
c        n        - Size of eigen problem

c     Inputs:
c        d(*)     - Eigenvalues
c        z(n,*)   - Eigenvalues of standard problem
c        ierr     - Error counter: zero if no errors

      implicit none
      integer  i,ii, j,jp1, k, l,l1, m,mml, n,n2, ierr
      real*8   b, c, f, g, h, hh, p, r, s
      real*8   scale,machep
      real*8   a(*),d(*),e(*),z(n,n)

      data machep/0.222044605d-15/
 
c     Eispac ql algorithm: Adapted from 'tred2' and 'tql2'

      n2 = 0
      do 100 i = 1,n
      do 100 j = 1,i
        n2     = n2 + 1
        z(i,j) = a(n2)
100   continue

      if(n.eq.1) go to 320
      n2 = n + 2
      do 300 ii = 2,n
        i     = n2 - ii
        l     = i - 1
        h     = 0.0d0
        scale = 0.0d0
        if(l.lt.2) go to 130
        do 120 k = 1,l
          scale = scale + abs(z(i,k))
120     continue
        if(scale.ne.0.0d0) go to 140
130     e(i) = z(i,l)
        go to 290
140     do 150 k = 1,l
          z(i,k) = z(i,k)/scale
          h      = h + z(i,k)*z(i,k)
150     continue
        f      = z(i,l)
        g      = -sign(sqrt(h),f)
        e(i)   = scale*g
        h      = h - f*g
        z(i,l) = f - g
        f      = 0.0d0
        do 240 j = 1,l
          z(j,i) = z(i,j)/h
          g      = 0.0d0
          do 180 k = 1,j
            g = g + z(j,k)*z(i,k)
180       continue
          jp1 = j + 1
          if(l.lt.jp1) go to 220
          do 200 k = jp1,l
            g = g + z(k,j)*z(i,k)
200       continue
220       e(j) = g/h
          f    = f + e(j)*z(i,j)
240     continue
        hh = f/(h+h)
        do 270 j = 1,l
          f    = z(i,j)
          g    = e(j) - hh*f
          e(j) = g
          do 260 k = 1,j
            z(j,k) = z(j,k) - f*e(k) - g*z(i,k)
260       continue
270     continue
290     d(i) = h
300   continue

c     Set transformation array for ql

320   d(1)   = z(1,1)
      z(1,1) = 1.0d0
      e(1)   = 0.0d0
      ierr   = 0
      if(n.eq.1) go to 950
      do 500 i = 2,n
        l = i - 1
        if(d(i).eq.0.0d0) go to 380
        do 370 j = 1,l
          g = 0.0d0
          do 340 k = 1,l
            g = g + z(i,k)*z(k,j)
340       continue
          do 360 k = 1,l
            z(k,j) = z(k,j) - g*z(k,i)
360       continue
370     continue
380     d(i)   = z(i,i)
        z(i,i) = 1.0d0
        do 400 j = 1,l
          z(i,j) = 0.0d0
          z(j,i) = 0.0d0
400     continue
500   continue

c     Begin 'ql' algorithm on tridagonal matrix now stored in 'd' and 'e

      do 600 i = 2,n
        e(i-1) = e(i)
600   continue
      f    = 0.0d0
      b    = 0.0d0
      e(n) = 0.0d0
      do 840 l = 1,n
        j = 0
        h = machep*(abs(d(l)) + abs(e(l)))
        if(b.lt.h) b = h
        do 710 m = l,n
          if(abs(e(m)).le.b) go to 720
710     continue
720     if(m.eq.l) go to 820
730     if(j.eq.30) go to 1000
        j    = j + 1
        l1   = l + 1
        g    = d(l)
        p    = (d(l1)-g)/(e(l)+e(l))
        r    = sqrt(p*p+1.0d0)
        d(l) = e(l)/(p+sign(r,p))
        h    = g - d(l)
        do 740 i = l1,n
          d(i) = d(i) - h
740     continue
        f   = f + h
        p   = d(m)
        c   = 1.0d0
        s   = 0.0d0
        mml = m - l
        do 800 ii = 1,mml
          i = m - ii
          g = c*e(i)
          h = c*p
          if(abs(p).ge.abs(e(i))) then
            c      = e(i)/p
            r      = sqrt(c*c+1.0d0)
            e(i+1) = s*p*r
            s      = c/r
            c      = 1.0d0/r
          else
            c      = p/e(i)
            r      = sqrt(c*c+1.0d0)
            e(i+1) = s*e(i)*r
            s      = 1.0d0/r
            c      = c*s
          endif
          p      = c*d(i) - s*g
          d(i+1) = h + s*(c*g + s*d(i))
          do 780 k = 1,n
            h        = z(k,i+1)
            z(k,i+1) = s*z(k,i) + c*h
            z(k,i  ) = c*z(k,i) - s*h
780       continue
800     continue
        e(l) = s*p
        d(l) = c*p
        if(abs(e(l)).gt.b) go to 730
820     d(l) = d(l) + f
840   continue

      do 900 ii = 2,n
        i = ii - 1
        k = i
        p = d(i)
        do 860 j = ii,n
        if(abs(d(j)).gt.abs(p)) then
          k = j
          p = d(j)
        endif
860     continue
        if(k.ne.i) then
          d(k) = d(i)
          d(i) = p
          do 880 j = 1,n
            p      = z(j,i)
            z(j,i) = z(j,k)
            z(j,k) = p
880       continue
        endif
900   continue
950   return

1000  ierr = l

      end

      subroutine geig(g,h,d,p,t,nev,prt)

c     Purpose: Solve general eigenproblem 'g*p = h*p*d'

c     Inputs:
c        g(*)     - Projected subspace array
c        h(*)     - Projected subspace array
c        t(*)     - Working vector
c        nev      - Size of subspace problem
c        prt      - Flag, output results of each iteration if true

c     Outputs:
c        d(*)     - Eigenvalues of subspace problem
c        p(nev,*) - Eigenvectors of subspace problem

      implicit none

      logical  prt
      integer  nev, ir

      real*8   g(*),h(*),d(*),p(nev,*),t(*)

c     Compute choleski factors of 'h'

      if(prt) call wprojm(g,nev,1)
      if(prt) call wprojm(h,nev,2)
      call chlfac(h,nev)

c     Compute standard eigenvalue problem matrix 'c'

      call chlfwd(h,g,p,nev)

c     Perform eignfunction decomposition of 'c'

      call eisql(g,d,t,p,nev,ir)

c     Compute vectors of original problem

      call chlbac(h,p,nev)

      end

      subroutine numass(b,neq,mq)

c     Purpose: Count number of non-zero entries in mass

c     Inputs
c        b(*)     - Diaonal mass
c        neq      - Length of b array

c     Outputs
c        nq       - Number of nonzero entries in b array

      implicit none
      integer  n,nn,neq,mq

      real*8   b(*)

      integer         ior,iow
      common /iofile/ ior,iow

      nn = 0
      do 10 n = 1,neq
        if(b(n).ne.0.0d0) nn = nn + 1
10    continue

      if(nn.lt.mq) write(iow,2000) nn
      if(ior.lt.0.and.nn.lt.mq) write(*,2000) nn

      mq = min(mq,nn)

2000  format(1x,'Subspace reduced to',i4,' by number of nonzero lumped',
     &          ' mass terms')

      end

      subroutine scalev(v,nn)

c     Purpose: Scale vector to have maximum element of 1.0

c     Inputs:
c        v(*)     - Vector to scale
c        nn       - Length of vector

c     Outputs:
c        v(*)     - Scaled vector

      implicit none
      integer n,nn

      real*8  v(*),vmax

      vmax = abs(v(1))
      do 100 n = 1,nn
        vmax = max(vmax,abs(v(n)))
100   continue

      do 110 n = 1,nn
        v(n) = v(n)/vmax
110   continue

      end

      subroutine sproja(v,z,g,neq,nev)

c     Compute subspace projection of 'aa' to form 'g'

c     Inputs:
c        aa(*)    - Tangent matrix to project
c        z(neq,*) - Iteration subspace (eigen) vectors
c        neq      - Number of equations in b(*)
c        nev      - Size of subspace 

c     Outputs:
c        g(*)     - Subspace projection of aa matrix

      implicit none

      integer  i, j, k, ma, neq,nev

      real*8   v(neq,*),z(neq,*),g(*)
      real*8   dot,engy

      integer         maxa
      parameter      (maxa=32000)
      real*8          aa
      common /adata/  aa(maxa)

      integer         maxf
      common /frdata/ maxf

      integer         iodr,iodw,ipd,ipr,ipi
      common /iofild/ iodr,iodw,ipd,ipr,ipi

      real*8          dimx,dimn
      integer                   nvb,npl
      common /nfrta/  dimx,dimn,nvb,npl

c     Forward reduce eigenvector estimates

      ma = maxf*nev + 1

c     Copy vectors 'v' into 'z'

      do 100 i = 1,nev
      do 100 j = 1,neq
        z(j,i) = v(j,i)
100   continue

c     Solve equations

      call rsolve(aa,z,aa(ma),ipd,ipr,maxf,nvb,neq,nev,engy,4)

c     Compute projection of stiffness

      k = 0
      do 200 j = 1,nev
      do 200 i = 1,j
        k    = k + 1
        g(k) = dot(v(1,i),z(1,j),neq)
200   continue

      end

      subroutine sprojb(b,v,t,h,neq,nev)

c     Purpose: Compute subspace projection of 'b' to form 'h'

c     Inputs:
c        b(*)     - Diagonal "mass" matrix
c        v(neq,*) - Iteration subspace (eigen) vectors
c        t(*)     - Working vector
c        neq      - Number of equations in b(*)
c        nev      - Size of subspace 

c     Outputs:
c        h(*)     - Subspace projection of b matrix

      implicit none
      integer  i, j, k, neq,nev
      real*8 b(1),v(neq,1),t(1),h(1),dot

c     Compute 'z' and 'b' projection to form 'h'

      do 130 j = 1,nev

c       Compute 'z' for a lumped mass

        do 100 i = 1,neq
          t(i) = v(i,j)*b(i)
100     continue

c       Project 'z' and 'v' vectors to form 'h'

        k = j*(j+1)/2
        do 110 i = j,nev
          h(k) = dot(t,v(1,i),neq)
          k    = k + i
110     continue
        do 120 i = 1,neq
          v(i,j) = t(i)
120     continue
130   continue

      end

      subroutine wprojm(a,nn,ia)

      implicit  none

      integer   i,ia, j, k, n,nn

      character ah(2)*1
      real*8    a(*)

      integer         ior,iow
      common /iofile/ ior,iow

      data ah(1),ah(2) /'g','h'/

      write(iow,2000) ah(ia)
      if(ior.lt.0) write(*,2000) ah(ia)

      i = 1
      do 100 n = 1,nn
        j = i + n - 1
        write(iow,2001) (a(k),k=i,j)
        if(ior.lt.0) write(*,2001) (a(k),k=i,j)
        i = i + n
100   continue

2000  format(1x,'Matrix ',a1)
2001  format(1p8d10.2)

      end
