program utmatu
c ***************** Program UTMATU *****************
c This program is designed to read in a real matrix,
c A, a real transformation matrix, B, perform the
c transformation: 
c        X = B(transpose) * A * B,
c and output the result.
c
c Copyright (C) 1992 Jeff Nichols and Jack Simons
c
c **************************************************
      implicit real (a-h,o-z)
      include "limits.h"
c
c     Set core needed to:
c         4 arrays of size (maxrow x maxcol)
c     total = 4*maxrow*maxcol
c
      dimension pcore(4*maxrow*maxcol),ipntr(10)
      itotal = 4*maxrow*maxcol + 10
      iask=1
      call modinfo('utmatu      ')
      call shwmem('utmatu      ',itotal)
      call goon
   10 continue
      write(*,*)
      write(*,*)' Input the row order of the first matrix (A). '
      read(*,*)iarow
      if(iarow.gt.maxrow)then
        write(*,1010)maxrow
        goto 10
      endif
   20 continue
      write(*,*)' Input the column order of the first matrix (A). '
      read(*,*)iacol
      if(iacol.gt.maxcol)then
        write(*,1010)maxcol
        goto 20
      endif
   30 continue
      write(*,*)' Input the row order of the second matrix (B). '
      read(*,*)ibrow
      if(ibrow.gt.maxrow)then
        write(*,1010)maxrow
        goto 30
      endif
   40 continue
      write(*,*)' Input the column order of the second matrix (B). '
      read(*,*)ibcol
      if(ibcol.gt.maxcol)then
        write(*,1010)maxcol
        goto 40
      endif
      if(iacol.ne.ibrow)then
        write(*,1020)
        call qmexit
      endif
      ipntr(1) = 1
      ipntr(2) = ipntr(1) + iarow*iacol
      ipntr(3) = ipntr(2) + ibrow*ibcol
      ipntr(4) = ipntr(3) + iarow*ibcol
      ipntr(5) = ipntr(4) + ibcol*ibcol
      call trans(pcore(ipntr(1)),pcore(ipntr(2)),pcore(ipntr(3)),
     &           pcore(ipntr(4)),iarow,iacol,ibrow,ibcol)
      call qmexit
      stop
 1010 format(' Maximum order possible is ',I3,'. ',/,
     &       ' Please try again. ')
 1020 format(' These matrices cannot be multiplied together. ',/,
     &       ' The column order of the first matrices is NOT ',/,
     &       ' compatible with the row order of the second. ',/,
     &       ' Start over. ')
      end
      subroutine trans(a,b,ab,btab,iarow,iacol,ibrow,ibcol)
      implicit real (a-h,o-z)
      logical yesno
      character*1 ans
      character*30 dfname, fname, dscrpt
      dimension a(iarow,iacol),
     &          b(ibrow,ibcol),
     &          ab(iarow,ibcol),
     &          btab(ibcol,ibcol)
   10 continue
c     
c     Read in the real matrix A (IAROW X IACOL) (disk).
c     
      dscrpt = 'input matrix to be transformed'
      ids = 30
      dfname = 'a_matrix.in'
      idf = 11
      fname = '                              '
      ifn = 30
      iflnum = 32
      call rfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno)
      if(yesno)then
        call rdmatd(fname,ifn,dscrpt,ids,a,iarow,iacol,
     &              iflnum)
c
c       Close input file.
c
        close(unit=32)
      else
        call rdmatu(dscrpt,ids,a,iarow,iacol)
      endif
c
c     Write input array A back to user.
c
      call wtmatu(dscrpt,ids,a,iarow,iacol)
      write(*,*)' Is this correct? <y or n> '
      read(*,'(a)')ans
      if(ans.eq.'Y'.or.ans.eq.'y')then
   40   continue
c     
c       Read in the real matrix B (IBROW X IBCOL) (disk).
c     
        dscrpt = 'transformation matrix'
        ids = 21
        dfname = 'b_matrix.in'
        idf = 11
        fname = '                              '
        ifn = 30
        iflnum = 33
        call rfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno)
        if(yesno)then
          call rdmatd(fname,ifn,dscrpt,ids,b,ibrow,ibcol,
     &                iflnum)
c
c         Close input file.
c
          close(unit=33)
        else
          call rdmatu(dscrpt,ids,b,ibrow,ibcol)
        endif
c
c       Write input array B back to user.
c
        call wtmatu(dscrpt,ids,b,ibrow,ibcol)
        write(*,*)' Is this correct? <y or n> '
        read(*,'(a)')ans
        if(ans.eq.'Y'.or.ans.eq.'y')then
          do 90 i = 1,iarow
            do 80 j = 1,ibcol
              ab(i,j) = 0.0e0
              do 70 k = 1,iacol
                ab(i,j) = ab(i,j) + a(i,k)*b(k,j)
   70         continue
   80       continue
   90     continue
          do 120 i = 1,ibcol
            do 110 j = 1,ibcol
              btab(i,j) = 0.0e0
              do 100 k = 1,ibrow
                btab(i,j) = btab(i,j) + b(k,i)*ab(k,j)
  100         continue
  110       continue
  120     continue
          dscrpt = 'transformed matrix'
          ids = 18
          call  wtmatu(dscrpt,ids,btab,ibcol,ibcol)
        else
          goto 40
        endif
      else
        goto 10
      endif
      dfname = 'utmatu.out'
      idf = 10
      fname = '                              '
      ifn = 30
      iflnum = 34
      call wfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno)
      if(yesno)then
        call wtmatd(fname,ifn,dscrpt,ids,btab,ibcol,ibcol,
     &              iflnum)
c
c       Close output file.
c
        close(unit=34)
      endif
      return
      end