c     ******************************************************************
c     MATRIXIO.F
c     Contains the following subroutines for simple formatted I/O of
c     vectors and matrices, where the prefix d_ stands for double
c     precision double precision and i_ stands for integer:
c     
c     d_read_matrix, d_write_matrix, dg_read_matrix, dg_write_matrix,
c     i_read_matrix, i_write_matrix.
c     
c     These routines can be used for I/O of vectors (one-dimensional
c     matrices) as follows:
c     -- {d,i}_write_matrix called with dimensions (1,n) prints a row
c        vector; when called with dimension (n,1), it writes a 
c        column vector.
c     -- {d,i}_read_matrix called with dimensions (1,n) or (n,1) will
c        read a vector.
c
c     itoa is a utility subroutine included to insert the value of
c     n into the passed format string.  For example, when n is 3, the
c     passed format '(f7.4tr2)' is expanded locally into '(3(f7.4 tr2))'.
c
c     No arguments in this file are modified upon return, except
c     string_n in itoa.
c     
c     There is currently no support for zero-dimensional matrices.
c
c     ******************************************************************
c     John Kerl
c     8/18/93
c     ******************************************************************

c     ******************************************************************
c     D_WRITE_MATRIX(matrix, m, n, form, len_form, type, offset)
c
c     ARGS:
c     matrix   -- The matrix to be written, of type double precision.
c     m,n      -- Dimensions of the matrix, with m the number of rows
c                 and n the number of columns.  If the matrix is upper-
c                 or lower-triangular, only m is used (since -- in a
c                 great violence to geometrical usage -- triangular
c                 matrices are assumed to be square).
c     form     -- Edit descriptor for matrix element, e.g. '(e11.4 tr1)'.
c     len_form -- May be passed as, for example, len ('(e11.4 tr1)').
c
c     (NOTE:  If either form is '' or len_form is 0, then a default is
c     used.  For d_write_matrix, this is '(f7.4 tr2)'; for
c     i_write_matrix, '(i4 tr2)'.)
c
c     type     -- 'ut' for upper-triangular; 'lt' for lower-triangular;
c                 'g' (or anything else) for a general matrix.  Perhaps
c                 I should also add support for band / tridiagonal
c                 matrices, and/or strides.
c     offset   -- Used only for UT matrices.  Number of spaces to indent
c                 to account for non-printed entries.  Aesthetically,
c                 should be equal to the float field width plus spacing
c                 between entries in the format string.  For example,
c                 with format '(f7.4 tr2)', offset should be 9.
c

c     ******************************************************************
c     This is a stub for general matrices.
      subroutine dg_write_matrix(matrix, m, n, form, len_form)
      implicit none
      integer m, n, len_form
      double precision matrix(m,n)
      character form
      call d_write_matrix(matrix, m, n, form, len_form, 'g', 0)
      end ! subroutine dg_write_matrix

c     ******************************************************************
      subroutine d_write_matrix(matrix, m, n, orig_form, orig_len_form,
     &     type, orig_offset)
      implicit none
      integer m, n, i, j, orig_len_form, len_form, len_n, orig_offset,
     &     offset
      double precision matrix(m,n)
      character orig_form, form*40, type*2
      character*20 string_n ! adequate for integers up to 64 bits
      character*60 local_form

c     orig_form and orig_len_form are used so as not to modify arguments.
c     form and len_form are local and we can play with them freely.

      len_form=orig_len_form
      form(1:len_form)=orig_form(1:len_form)
      offset=orig_offset

      if ((form .eq. '') .or. (len_form .eq. 0)) then
         len_form=10
         form(1:len_form)='(f7.4tr2)'
      end if
      if (offset .eq. -1) offset=9

c     All code up to the do-loop serves merely to insert the value of
c     n and parentheses into the format string: e.g. '(e11.4 tr2)' 
c     becomes '(3(e11.4 tr2))'.  Allowing the format string to be passed
c     without the dimension prefix allows flexibile, repeated use of
c     a single format string within the calling routine.

      call itoa(n, string_n, len_n)

      local_form(1               :1               ) = '('
      local_form(2               :len_n+1         ) = string_n
      local_form(len_n+2         :len_n+len_form+1) = form(1:len_form)
      local_form(len_n+len_form+2:len_n+len_form+2) = ')'

c     write (*, '(a a)')
c    &  'Local format string is ', local_form(1:len_form+len_n+2)

      if (type .eq. 'lt') then
         do i=1,m
            write (*, local_form) (matrix(i,j), j=1,i)
         end do
      else if (type .eq. 'ut') then
         do i=1,m
            write (*, '(<(i-1)*offset>a1$)') (' ', j=1, (i-1)*offset)
            write (*, local_form) (matrix(i,j), j=i,m)
         end do
      else
         do i=1,m
            write (*, local_form) (matrix(i,j), j=1,n)
         end do
      end if

      end ! subroutine d_write_matrix

c     ******************************************************************
c     This is a stub for general matrices.
      subroutine dg_read_matrix(matrix, m, n)
      implicit none
      integer m, n
      double precision matrix(m,n)
      call d_read_matrix(matrix, m, n, 'g')
      end ! subroutine dg_read_matrix
c     ******************************************************************
      subroutine d_read_matrix(matrix, m, n, type)
      implicit none
c     Just a plain-vanilla list-directed read; coded as a function
c     call in order to hide low-level detail (e.g. loop integers).

      integer m,n,i,j
      double precision matrix(m,n)
      character*2 type

      if (type .eq. 'ut') then
         read (*,*) ((matrix(i,j), j=i,n), i=1,n)
      else if (type .eq. 'lt') then
         read (*,*) ((matrix(i,j), j=1,i), i=1,m)
      else
         read (*,*) ((matrix(i,j), j=1,n), i=1,m)
      end if

      end ! subroutine d_read_matrix

c     ******************************************************************
      subroutine i_write_matrix(matrix, m, n, orig_form, orig_len_form,
     &     type, offset)
      implicit none
      integer m, n, i, j, orig_len_form, len_form, len_n, offset
      integer matrix(m,n)
      character orig_form, form*40, type*2
      character*20 string_n ! adequate for integers up to 64 bits
      character*60 local_form

c     orig_form and orig_len_form are used so as not to modify arguments.
c     form and len_form are local and we can play with them freely.

      len_form=orig_len_form
      form(1:len_form)=orig_form(1:len_form)

      if ((form .eq. '') .or. (len_form .eq. 0)) then
         len_form=10
         form(1:len_form)='(f7.4tr2)'
      end if

c     All code up to the do-loop serves merely to insert the value of
c     n and parentheses into the format string: e.g. '(e11.4 tr2)' 
c     becomes '(3(e11.4 tr2))'.  Allowing the format string to be passed
c     without the dimension prefix allows flexibile, repeated use of
c     a single format string within the calling routine.

      call itoa(n, string_n, len_n)

      local_form(1               :1               ) = '('
      local_form(2               :len_n+1         ) = string_n
      local_form(len_n+2         :len_n+len_form+1) = form(1:len_form)
      local_form(len_n+len_form+2:len_n+len_form+2) = ')'

c     write (*, '(a a)')
c    &  'Local format string is ', local_form(1:len_form+len_n+2)

      if (type .eq. 'lt') then
         do i=1,m
            write (*, local_form) (matrix(i,j), j=1,i)
         end do
      else if (type .eq. 'ut') then
         do i=1,m
            write (*, '(<(i-1)*offset>a1$)') (' ', j=1, (i-1)*offset)
            write (*, local_form) (matrix(i,j), j=1,m)
         end do
      else
         do i=1,m
            write (*, local_form) (matrix(i,j), j=1,n)
         end do
      end if

      end ! subroutine i_write_matrix

c     ******************************************************************
      subroutine i_read_matrix(matrix, m, n, type)
      implicit none
      integer m,n,i,j
      integer matrix(m,n)
      character*2 type

      if (type .eq. 'ut') then
         read (*,*) ((matrix(i,j), j=i,n), i=1,n)
      else if (type .eq. 'lt') then
         read (*,*) ((matrix(i,j), j=1,i), i=1,m)
      else
         read (*,*) ((matrix(i,j), j=1,n), i=1,m)
      end if

      end ! subroutine i_read_matrix

c     ******************************************************************
c     ITOA
c     Converts an integer to a character string.
c     
c     ARGS:
c     int_n    -- The integer to be converted.
c     string_n -- Holds the result.
c     length   -- Holds the number of digits in int_n.  Thus the result
c                 may be referenced in the calling routine as 
c                 string_n(1:length)
c     ******************************************************************
      subroutine itoa(int_n, string_n, length)
      implicit none
      integer int_n, length
      character*20 string_n ! adequate for up to 64-bit integers
      integer local_n, i

      length = int(log10(dble(int_n))) + 1
      local_n = int_n
      do i=length-1, 0, -1
         string_n(length-i:length-i) = 
     &        achar((local_n / 10**i) + ichar('0'))
         local_n = mod(local_n, 10**i)
      end do

      end ! subroutine itoa

c     ** end of file ***************************************************
