program matxmat
c ***************** Program MATXMAT ****************
c This program is designed to read in two real
c matrices; A and B, and to multiply them together: 
c             AB = 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         3 arrays of size (maxrow x maxcol)
c     total = 3*maxrow*maxcol)
c
      dimension pcore(3*maxrow*maxcol),ipntr(10)
      itotal = 3*maxrow*maxcol + 10
      call modinfo('matxmat     ')
      call shwmem('matxmat     ',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
      call multiply(pcore(ipntr(1)),pcore(ipntr(2)),pcore(ipntr(3)),
     &              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 multiply(a,b,prod,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),
     &          prod(iarow,ibcol)
   10 continue
c     
c     Read in the real matrix A (IAROW X IACOL) (disk).
c     
      dscrpt = 'input matrix, A,'
      ids = 16
      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 = 'input matrix, B,'
        ids = 16
        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
              prod(i,j) = 0.0d0
              do 70 k = 1,iacol
                prod(i,j) = prod(i,j) + a(i,k)*b(k,j)
   70         continue
   80       continue
   90     continue
          dscrpt = 'product matrix (A X B)'
          ids = 22
          call  wtmatu(dscrpt,ids,prod,iarow,ibcol)
        else
          goto 40
        endif
      else
        goto 10
      endif
      dfname = 'matxmat.out'
      idf = 11
      fname = '                              '
      ifn = 30
      iflnum = 34
      call wfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno)
      if(yesno)then
        call wtmatd(fname,ifn,dscrpt,ids,prod,iarow,ibcol,
     &              iflnum)
c
c       Close output file.
c
        close(unit=34)
      endif
      return
      end