program fnct_mat
c **************** Program FNCT_MAT ***************
c This program is designed to read in a real square
c matrix, perform a function on it, and return this
c new array.  Possible functions, using X as the
c input matrix, are: 
c 
c (1) X^(-1/2), NOTE: X must be real symmetric, 
c                     and positive definite. 
c 
c (2) X^(+1/2), NOTE: X must be real symmetric, 
c                     and positive definite. 
c 
c (3) X^(-1), NOTE: X must be real symmetric, 
c                   and have non-zero eigenvalues.
c 
c (4) a power series expansion of a matrix 
c     to find the transformation matrix: 
c     U = exp(X) = 1 + X + X**2/2! + X**3/3! + 
c                            ... + X**N/N! + ...
c
c Copyright (C) 1992 Jeff Nichols and Jack Simons
c
c **************************************************
      implicit real (a-h,o-z)
      include "limits.h"
c
c     Core needed for matrix exp(+  or - 1/2):
c         1 symmetric array (input) of order mmax
c         1 square array (eigenvectors) of order mmax
c         1 vector (eigenvalues) of length mmax
c         2 scratch vectors of length mmax
c         1 scratch array of order mmax
c
c     Core needed for e exp(matrix) power series expansion:
c         4 square arrays of order mmax
c     total = 4*mmax*mmax
c
c     Dimension PCORE for largest possible function case.
c
      dimension pcore(4*mmax*mmax),ipntr(10)
c
c     Calculate and show memory use.
c
      itotal = 4*mmax**2 + 10
      call modinfo('funct_mat   ')
      call shwmem('funct_mat   ',itotal)
      call goon
   10 continue
      write(*,*)
      write(*,*)' Input the order of the your input array. '
      read(*,*)msize
      if(msize.gt.mmax)then
        write(*,1010)mmax
        goto 10
      endif
   20 continue
      write(*,*)'  ' 
      write(*,*)' Input the function you would like to perform '
      write(*,*)' on your input matrix, X. '
      write(*,*)'  '
      write(*,*)' (1) X^(-1/2), NOTE: X must be real symmetric,  '
      write(*,*)'                     and positive definite.  '
      write(*,*)'  '
      write(*,*)' (2) X^(+1/2), NOTE: X must be real symmetric,  '
      write(*,*)'                     and positive definite.  '
      write(*,*)'  '
      write(*,*)' (3) X^(-1), NOTE: X must be real symmetric,  '
      write(*,*)'                   and have non-zero eigenvalues. '
      write(*,*)'  '
      write(*,*)' (4) a power series expansion of a matrix  '
      write(*,*)'     to find the transformation matrix:  '
      write(*,*)'     U = exp(X) = 1 + X + X**2/2! + X**3/3! +  '
      write(*,*)'                            ... + X**N/N! + ... '
      write(*,*)'  '
      write(*,*)' (5) QUIT '
      read(*,*)iftype
      if(iftype.eq.1.or.iftype.eq.2.or.iftype.eq.3)then
        ipntr(1) = 1
        ipntr(2) = ipntr(1) + (msize*(msize+1))/2
        ipntr(3) = ipntr(2) + msize*msize
        ipntr(4) = ipntr(3) + msize
        ipntr(5) = ipntr(4) + msize
        ipntr(6) = ipntr(5) + msize
        ipntr(7) = ipntr(6) + msize*msize
        call funct1(pcore(ipntr(1)),pcore(ipntr(2)),pcore(ipntr(3)),
     &              pcore(ipntr(4)),pcore(ipntr(5)),pcore(ipntr(6)),
     &              msize,iftype)
      elseif(iftype.eq.4)then
        ipntr(1) = 1
        ipntr(2) = ipntr(1) + msize*msize
        ipntr(3) = ipntr(2) + msize*msize
        ipntr(4) = ipntr(3) + msize*msize
        ipntr(5) = ipntr(4) + msize*msize
        call funct2(pcore(ipntr(1)),pcore(ipntr(2)),pcore(ipntr(3)),
     &              pcore(ipntr(4)),msize)
      elseif(iftype.eq.5)then
        call qmexit
      else
        write(*,*)' Please input one of the appropriate reponses. '
        goto 20
      endif
      call qmexit
      stop
 1010 format(' Maximum order possible is less than ',I2,/,
     &       ' Please try again. ')
      end
      subroutine funct1(asym,avects,avals,scr1,scr2,scr3,
     &                  msize,iftype)
      implicit real (a-h,o-z)
      logical yesno
      character*1 ans
      character*30 dfname, fname, dscrpt
      dimension asym((msize*(msize+1))/2),
     &          avects(msize,msize),
     &          avals(msize),
     &          scr1(msize),
     &          scr2(msize),
     &          scr3(msize,msize)
   10 continue
      dscrpt = 'input matrix'
      ids = 12
      dfname = 'funct_mat.in'
      idf = 12
      fname = '                              '
      ifn = 30
      iflnum = 32
      call rfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno)
      if(yesno)then
        call rdmatd(fname,ifn,dscrpt,ids,scr3,msize,msize,
     &              iflnum)
c
c       Close input matrix file.
c
        close(unit=32)
      else
        call rdmatu(dscrpt,ids,scr3,msize,msize)
      endif
c
c     Put the input array into a symmetric array.
c      
      ij = 0
      do 40 i = 1,msize
        do 30 j = 1,i
          ij = ij + 1
          asym(ij) = scr3(i,j)
   30   continue
   40 continue
c
c     Write input array back to user.
c
      write(*,*)' The array to perform the function on is: '
      call tprntd(asym,msize)
      write(*,*)' is this correct? <y or n> '
      read(*,'(a)')ans
      if(ans.eq.'Y'.or.ans.eq.'y')then
c
c       Diagonalize the array.
c
        call diaghe(asym,msize,avals,msize,avects,msize,scr1)
c
c       Write the eigenvalues to user.
c
        write(*,*)' The eigenvalues of this matrix are: '
        call matout(avals,msize,1,msize)
c
c       write the eigenvectors to user.
c
        write(*,*)' their corresponding eigenvectors are: '
        call matout(avects,msize,msize,msize)
      else
        write(*,*)' Please try again. '
        goto 10
      endif
c
c     Construct ASYM exp(-1/2 or +1/2 or -1).
c
      do 70 i = 1,msize
        if(iftype.eq.1)then
          avals(i) = sqrt(avals(i))
          avals(i) = 1.0e0/avals(i)
        elseif(iftype.eq.2)then
          avals(i) = sqrt(avals(i))
        elseif(iftype.eq.3)then
          avals(i) = 1.0e0/avals(i)
        endif
   70 continue
      do 90 i = 1,msize
        do 80 j = 1,msize
          scr3(i,j) = avals(i)*avects(j,i)
   80   continue
   90 continue
      ij = 0
      do 120 i = 1,msize
        do 110 j = 1,i
          ij = ij + 1
          asym(ij) = 0.0e0
          do 100 k = 1,msize
            asym(ij) = asym(ij) + avects(i,k)*scr3(k,j)
  100     continue
  110   continue
  120 continue
      if(iftype.eq.1)then
        write(*,*)' The input matrix exp(-1/2) is: '
        dscrpt = 'input matrix exp(-1/2)'
        ids = 22
      elseif(iftype.eq.2)then
        write(*,*)' The input matrix exp(+1/2) is: '
        dscrpt = 'input matrix exp(+1/2)'
        ids = 22
      elseif(iftype.eq.3)then
        write(*,*)' The input matrix exp(-1) is: '
        dscrpt = 'input matrix exp(-1)'
        ids = 20
      endif
      call tprntd(asym,msize)
      dfname = 'funct_mat.out'
      idf = 13
      fname = '                              '
      ifn = 30
      iflnum = 32
      call wfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno)
      if(yesno)then
c
c       Move the symmetric output array into a full size 
c       square array (SCR3).
c      
        ij = 0
        do 140 i = 1,msize
          do 130 j = 1,i
            ij = ij + 1
            scr3(i,j) = asym(ij)
            scr3(j,i) = asym(ij)
  130     continue
  140   continue
        call wtmatd(fname,ifn,dscrpt,ids,scr3,msize,msize,
     &              iflnum)
c
c       Close output matrix file.
c
        close(unit=32)
      else
        return
      endif
      return
      end
      subroutine funct2(x,u,scr1,scr2,msize)
      implicit real (a-h,o-z)
      logical yesno
      character*1 ans
      character*30 dfname, fname, dscrpt
      dimension x(msize,msize),
     &          u(msize,msize),
     &          scr1(msize,msize),
     &          scr2(msize,msize)
c
c     Define ACRCY ... machine dependent convergence tolerance.
c
      acrcy = 1.0e-7
   10 continue
      dscrpt = 'input matrix'
      ids = 12
      dfname = 'funct_mat.in'
      idf = 12
      fname = '                              '
      ifn = 30
      iflnum = 32
      call rfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno)
      if(yesno)then
        call rdmatd(fname,ifn,dscrpt,ids,x,msize,msize,
     &              iflnum)
c
c       Close input matrix file.
c
        close(unit=32)
      else
        call rdmatu(dscrpt,ids,x,msize,msize)
      endif
c
c     Write input array back to user.
c
      call wtmatu(dscrpt,ids,x,msize,msize)
      write(*,*)' Is this correct? <y or n> '
      read(*,'(a)')ans
      if(ans.eq.'Y'.or.ans.eq.'y')then
c
c       Put unit matrix in U and scratch array SCR1.
c
        do 40 i = 1,msize
          u(i,i) = 1.0e0
          scr1(i,i) = 1.0e0
   40   continue
c
c       Perform the expansion until convergence is reached.
c       The index K will specify expansion power.  Arbitrarily
c       choose some large number for possible expansion limit.
c
        do 120 k = 1,100
c
c         copy scr1 to scr2.
c
          do 60 i = 1,msize
            do 50 j = 1,msize
              scr2(i,j) = scr1(i,j)
   50       continue
   60     continue
c
c         Calculate the Kth order expansion term (SCR1):
c
c             scr1  = (1/k) (scr2) x (x)
c
c             where scr2 = (k-1)st expansion term
c
c         NOTE: This is forming successive powers of RKAPPA
c               in SCR1
c
          rk=float(k)
          do 90 i = 1,msize
            do 80 j = 1,msize
              scr1(i,j) = 0.0e0
              do 70 l = 1,msize
                scr1(i,j) = scr1(i,j) + scr2(i,l)*x(l,j)
   70         continue
              scr1(i,j) = scr1(i,j)/rk
   80       continue
   90     continue
          diff=0.0e0
c
c         Sum this Kth order expansion term into U.
c
          do 110 i = 1,msize
            do 100 j = 1,msize
              u(i,j) = u(i,j) + scr1(i,j)
              if(abs(scr1(i,j)).gt.diff)diff=abs(scr1(i,j))
  100       continue
  110     continue
          if(diff.le.acrcy)goto 130
  120   continue
  130   continue
        write(*,1000)acrcy,k
        dscrpt = 'e exp(input matrix)'
        ids = 19
        call  wtmatu(dscrpt,ids,u,msize,msize)
        dfname = 'funct_mat.out'
        idf = 13
        fname = '                              '
        ifn = 30
        iflnum = 32
        call wfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno)
        if(yesno)then
          call wtmatd(fname,ifn,dscrpt,ids,u,msize,msize,
     &                iflnum)
c
c         Close output matrix file.
c
          close(unit=32)
        endif
      else
        write(*,*)' Please try again. '
        goto 10
      endif
      return
 1000 format(' The matrix e exp(input matrix) stable to ',
     &       F11.8,' was ',/,' achieved after a ',I2,
     &       ' term power series expansion. ')
      end