c $Id: inout.f,v 1.1 2008-04-11 06:01:06 geuzaine Exp $
c----------------------------------------------------------------------c
      subroutine psplot (ncol,ja,ia,iunt,mode)
      integer ja(*),ia(*),iunt,ncol,id,mode
      call pspltm (ncol, ncol, mode, ja, ia, ' ',
     &     0, 5.0, "in", 0, id, iunt)
      return
      end
c----------------------------------------------------------------------c
c                          S P A R S K I T                             c
c----------------------------------------------------------------------c
C                        INPUT-OUTPUT MODULE                           c
c----------------------------------------------------------------------c
c contents:                                                            c
c----------                                                            c
c  readmt : reads matrices in the Boeing/Harwell format.               c
c  prtmt  : prints matrices in the Boeing/Harwell format.              c
c  dump   : outputs matrix rows in a simple format (debugging purposes)c 
c  pspltm : generates a post-script plot of the non-zero pattern of A  c
c  pltmt  : produces a 'pic' file for plotting a sparse matrix         c
c  smms   : write the matrx in a format used in SMMS package           c
c  readsm : reads matrics in coordinate format (as in SMMS package)    c
c  readsk : reads matrices in CSR format (simplified H/B formate).     c
c  skit   : writes matrics to a file, format same as above.            c
c  prtunf : writes matrics (in CSR format) unformatted                 c
c  readunf: reads unformatted data of matrics (in CSR format)          c
c----------------------------------------------------------------------c
	subroutine readmt (nmax,nzmax,job,iounit,a,ja,ia,rhs,nrhs,
     *                     guesol,nrow,ncol,nnz,title,key,type,ierr)
c-----------------------------------------------------------------------
c this subroutine reads  a boeing/harwell matrix. handles right hand 
c sides in full format only (no sparse right hand sides).
c Also the matrix must be in assembled forms.
c Author: Youcef Saad - Date: Sept., 1989
c         updated Oct 31, 1989.
c-----------------------------------------------------------------------
c on entry:
c---------
c nmax 	 =  max column dimension  allowed for matrix. The array ia should 
c	    be of length at least ncol+1 (see below) if job.gt.0
c nzmax	 = max number of nonzeros elements allowed. the arrays a, 
c          and ja should be of length equal to nnz (see below) if these
c          arrays are to be read (see job).
c          
c job	 = integer to indicate what is to be read. (note: job is an
c          input and output parameter, it can be modified on return)
c          job = 0    read the values of ncol, nrow, nnz, title, key,
c                     type and return. matrix is not read and arrays
c                     a, ja, ia, rhs are not touched.
c          job = 1    read srtucture only, i.e., the arrays ja and ia.
c          job = 2    read matrix including values, i.e., a, ja, ia
c          job = 3    read matrix and right hand sides: a,ja,ia,rhs.
c		      rhs may contain initial guesses and exact 
c                     solutions appended to the actual right hand sides.
c		      this will be indicated by the output parameter
c                     guesol [see below]. 
c                     
c nrhs   = integer. nrhs is an input as well as ouput parameter.
c          at input nrhs contains the total length of the array rhs.
c          See also ierr and nrhs in output parameters.
c
c iounit = logical unit number where to read the matrix from.
c
c on return:
c---------- 
c job    = on return job may be modified to the highest job it could
c          do: if job=2 on entry but no matrix values are available it
c          is reset to job=1 on return. Similarly of job=3 but no rhs 
c          is provided then it is rest to job=2 or job=1 depending on 
c          whether or not matrix values are provided.
c          Note that no error message is triggered (i.e. ierr = 0 
c          on return in these cases. It is therefore important to
c          compare the values of job on entry and return ).
c
c a	 = the a matrix in the a, ia, ja (column) storage format
c ja 	 = column number of element a(i,j) in array a.
c ia     = pointer  array. ia(i) points to the beginning of column i.
c
c rhs    = real array of size nrow + 1 if available (see job)
c
c nrhs   = integer containing the number of right-hand sides found
c          each right hand side may be accompanied with an intial guess
c          and also the exact solution.
c
c guesol = a 2-character string indicating whether an initial guess 
c          (1-st character) and / or the exact solution (2-nd
c          character) is provided with the right hand side.
c	   if the first character of guesol is 'G' it means that an
c          an intial guess is provided for each right-hand side.
c          These are appended to the right hand-sides in the array rhs.
c	   if the second character of guesol is 'X' it means that an
c          exact solution is provided for each right-hand side.
c          These are  appended to the right hand-sides 
c          and the initial guesses (if any) in the array rhs.
c
c nrow   = number of rows in matrix
c ncol	 = number of columns in matrix 
c nnz	 = number of nonzero elements in A. This info is returned
c          even if there is not enough space in a, ja, ia, in order
c          to determine the minimum storage needed. 
c
c title  = character*72 = title of matrix test ( character a*72). 
c key    = character*8  = key of matrix 
c type   = charatcer*3  = type of matrix.
c          for meaning of title, key and type refer to documentation 
c          Harwell/Boeing matrices.
c
c ierr   = integer used for error messages 
c         * ierr  =  0 means that  the matrix has been read normally. 
c         * ierr  =  1 means that  the array matrix could not be read 
c         because ncol+1 .gt. nmax
c         * ierr  =  2 means that  the array matrix could not be read 
c         because nnz .gt. nzmax 
c         * ierr  =  3 means that  the array matrix could not be read 
c         because both (ncol+1 .gt. nmax) and  (nnz .gt. nzmax )
c         * ierr  =  4 means that  the right hand side (s) initial 
c         guesse (s) and exact solution (s)   could  not be
c         read because they are stored in sparse format (not handled
c         by this routine ...) 
c         * ierr  =  5 means that the right-hand-sides, initial guesses
c         and exact solutions could not be read because the length of 
c         rhs as specified by the input value of nrhs is not 
c         sufficient to store them. The rest of the matrix may have
c         been read normally.
c 
c Notes:
c-------
c 1) The file inout must be open (and possibly rewound if necessary)
c    prior to calling readmt.
c 2) Refer to the documentation on the Harwell-Boeing formats
c    for details on the format assumed by readmt.
c    We summarize the format here for convenience.
c  
c    a) all lines in inout are assumed to be 80 character long.
c    b) the file consists of a header followed by the block of the 
c       column start pointers followed by the block of the
c       row indices, followed by the block of the real values and
c       finally the numerical values of the right-hand-side if a 
c       right hand side is supplied. 
c    c) the file starts by a header which contains four lines if no
c       right hand side is supplied and five lines otherwise.
c       * first line contains the title (72 characters long) followed by
c         the 8-character identifier (name of the matrix, called key)
c        [ A72,A8 ]
c       * second line contains the number of lines for each
c         of the following data blocks (4 of them) and the total number 
c         of lines excluding the header.
c        [5i4]
c       * the third line contains a three character string identifying
c         the type of matrices as they are referenced in the Harwell
c         Boeing documentation [e.g., rua, rsa,..] and the number of
c         rows, columns, nonzero entries.
c         [A3,11X,4I14]
c       * The fourth line contains the variable fortran format
c         for the following data blocks.
c         [2A16,2A20] 
c       * The fifth line is only present if right-hand-sides are 
c         supplied. It consists of three one character-strings containing
c         the storage format for the right-hand-sides 
c         ('F'= full,'M'=sparse=same as matrix), an initial guess 
c         indicator ('G' for yes), an exact solution indicator 
c         ('X' for yes), followed by the number of right-hand-sides
c         and then the number of row indices. 
c         [A3,11X,2I14] 
c     d) The three following blocks follow the header as described 
c        above.
c     e) In case the right hand-side are in sparse formats then 
c        the fourth block uses the same storage format as for the matrix
c        to describe the NRHS right hand sides provided, with a column
c        being replaced by a right hand side.
c-----------------------------------------------------------------------
      character title*72, key*8, type*3, ptrfmt*16, indfmt*16,
     1       valfmt*20, rhsfmt*20, rhstyp*3, guesol*2
      integer totcrd, ptrcrd, indcrd, valcrd, rhscrd, nrow, ncol,
     1     nnz, neltvl, nrhs, nmax, nzmax, nrwindx
      integer ia (nmax+1), ja (nzmax) 
      real*8 a(nzmax), rhs(*) 
c-----------------------------------------------------------------------
      ierr = 0
      lenrhs = nrhs
c
      read (iounit,10) title, key, totcrd, ptrcrd, indcrd, valcrd, 
     1     rhscrd, type, nrow, ncol, nnz, neltvl, ptrfmt, indfmt, 
     2     valfmt, rhsfmt
 10   format (a72, a8 / 5i14 / a3, 11x, 4i14 / 2a16, 2a20)
c
      if (rhscrd .gt. 0) read (iounit,11) rhstyp, nrhs, nrwindx
 11   format (a3,11x,i14,i14)
c
c anything else to read ?
c
      if (job .le. 0) return
c     ---- check whether matrix is readable ------ 
      n = ncol
      if (ncol .gt. nmax) ierr = 1
      if (nnz .gt. nzmax) ierr = ierr + 2
      if (ierr .ne. 0) return
c     ---- read pointer and row numbers ---------- 
      read (iounit,ptrfmt) (ia (i), i = 1, n+1)
      read (iounit,indfmt) (ja (i), i = 1, nnz)
c     --- reading values of matrix if required....
      if (job .le. 1)  return
c     --- and if available ----------------------- 
      if (valcrd .le. 0) then
	 job = 1
	 return
      endif
      read (iounit,valfmt) (a(i), i = 1, nnz)
c     --- reading rhs if required ---------------- 
      if (job .le. 2)  return
c     --- and if available ----------------------- 
      if ( rhscrd .le. 0) then
	 job = 2
	 return
      endif
c     
c     --- read right-hand-side.-------------------- 
c     
      if (rhstyp(1:1) .eq. 'M') then 
         ierr = 4
         return
      endif
c
      guesol = rhstyp(2:3) 
c     
      nvec = 1 
      if (guesol(1:1) .eq. 'G' .or. guesol(1:1) .eq. 'g') nvec=nvec+1
      if (guesol(2:2) .eq. 'X' .or. guesol(2:2) .eq. 'x') nvec=nvec+1
c     
      len = nrhs*nrow 
c     
      if (len*nvec .gt. lenrhs) then
         ierr = 5
         return
      endif
c
c read right-hand-sides
c
      next = 1
      iend = len
      read(iounit,rhsfmt) (rhs(i), i = next, iend)
c
c read initial guesses if available
c
      if (guesol(1:1) .eq. 'G' .or. guesol(1:1) .eq. 'g') then
         next = next+len
         iend = iend+ len
         read(iounit,valfmt) (rhs(i), i = next, iend)
      endif
c     
c read exact solutions if available
c
      if (guesol(2:2) .eq. 'X' .or. guesol(2:2) .eq. 'x') then
         next = next+len
         iend = iend+ len
         read(iounit,valfmt) (rhs(i), i = next, iend)
      endif
c     
      return
c--------- end of readmt -----------------------------------------------
c----------------------------------------------------------------------- 
      end
c-----------------------------------------------------------------------
      subroutine prtmt (nrow,ncol,a,ja,ia,rhs,guesol,title,key,type,
     1     ifmt,job,iounit)
c-----------------------------------------------------------------------
c writes a matrix in Harwell-Boeing format into a file.
c assumes that the matrix is stored in COMPRESSED SPARSE COLUMN FORMAT.
c some limited functionality for right hand sides. 
c Author: Youcef Saad - Date: Sept., 1989 - updated Oct. 31, 1989 to
c cope with new format. 
c-----------------------------------------------------------------------
c on entry:
c---------
c nrow   = number of rows in matrix
c ncol	 = number of columns in matrix 
c a	 = real*8 array containing the values of the matrix stored 
c          columnwise
c ja 	 = integer array of the same length as a containing the column
c          indices of the corresponding matrix elements of array a.
c ia     = integer array of containing the pointers to the beginning of 
c	   the row in arrays a and ja.
c rhs    = real array  containing the right-hand-side (s) and optionally
c          the associated initial guesses and/or exact solutions
c          in this order. See also guesol for details. the vector rhs will
c          be used only if job .gt. 2 (see below). Only full storage for
c          the right hand sides is supported. 
c
c guesol = a 2-character string indicating whether an initial guess 
c          (1-st character) and / or the exact solution (2-nd)
c          character) is provided with the right hand side.
c	   if the first character of guesol is 'G' it means that an
c          an intial guess is provided for each right-hand sides. 
c          These are assumed to be appended to the right hand-sides in 
c          the array rhs.
c	   if the second character of guesol is 'X' it means that an
c          exact solution is provided for each right-hand side.
c          These are assumed to be appended to the right hand-sides 
c          and the initial guesses (if any) in the array rhs.
c
c title  = character*72 = title of matrix test ( character a*72 ).
c key    = character*8  = key of matrix 
c type   = charatcer*3  = type of matrix.
c
c ifmt	 = integer specifying the format chosen for the real values
c	   to be output (i.e., for a, and for rhs-guess-sol if 
c          applicable). The meaning of ifmt is as follows.
c	  * if (ifmt .lt. 100) then the D descriptor is used,
c           format Dd.m, in which the length (m) of the mantissa is 
c           precisely the integer ifmt (and d = ifmt+6)
c	  * if (ifmt .gt. 100) then prtmt will use the 
c           F- descriptor (format Fd.m) in which the length of the 
c           mantissa (m) is the integer mod(ifmt,100) and the length 
c           of the integer part is k=ifmt/100 (and d = k+m+2)
c	    Thus  ifmt= 4   means  D10.4  +.xxxxD+ee    while
c	          ifmt=104  means  F7.4   +x.xxxx
c	          ifmt=205  means  F9.5   +xx.xxxxx
c	    Note: formats for ja, and ia are internally computed.
c
c job	 = integer to indicate whether matrix values and
c	   a right-hand-side is available to be written
c          job = 1   write srtucture only, i.e., the arrays ja and ia.
c          job = 2   write matrix including values, i.e., a, ja, ia
c          job = 3   write matrix and one right hand side: a,ja,ia,rhs.
c	   job = nrhs+2 write matrix and nrhs successive right hand sides
c	   Note that there cannot be any right-hand-side if the matrix
c	   has no values. Also the initial guess and exact solutions when 
c          provided are for each right hand side. For example if nrhs=2 
c          and guesol='GX' there are 6 vectors to write.
c          
c
c iounit = logical unit number where to write the matrix into.
c
c on return:
c---------- 
c the matrix a, ja, ia will be written in output unit iounit
c in the Harwell-Boeing format. None of the inputs is modofied.
c  
c Notes: 1) This code attempts to pack as many elements as possible per
c        80-character line. 
c        2) this code attempts to avoid as much as possible to put
c        blanks in the formats that are written in the 4-line header
c	 (This is done for purely esthetical reasons since blanks
c        are ignored in format descriptors.)
c        3) sparse formats for right hand sides and guesses are not
c        supported.
c-----------------------------------------------------------------------
      character title*72,key*8,type*3,ptrfmt*16,indfmt*16,valfmt*20,
     *	        guesol*2, rhstyp*3
      integer totcrd, ptrcrd, indcrd, valcrd, rhscrd, nrow, ncol,
     1     nnz, nrhs, len, nperli, nrwindx
      integer ja(*), ia(*) 	
      real*8 a(*),rhs(*)
c--------------
c     compute pointer format
c--------------
      nnz    = ia(ncol+1) -1
      len    = int ( alog10(0.1+real(nnz+1))) + 1
      nperli = 80/len
      ptrcrd = ncol/nperli + 1
      if (len .gt. 9) then
         assign 101 to ix
      else
         assign 100 to ix
      endif
      write (ptrfmt,ix) nperli,len
 100  format(1h(,i2,1HI,i1,1h) )
 101  format(1h(,i2,1HI,i2,1h) )
c----------------------------
c compute ROW index format
c----------------------------
      len    = int ( alog10(0.1+real(nrow) )) + 1
      nperli = min0(80/len,nnz)
      indcrd = (nnz-1)/nperli+1
      write (indfmt,100) nperli,len
c---------------
c compute values and rhs format (using the same for both)
c--------------- 
      valcrd	= 0
      rhscrd  = 0
c quit this part if no values provided.
      if (job .le. 1) goto 20
c     
      if (ifmt .ge. 100) then
         ihead = ifmt/100
         ifmt = ifmt-100*ihead
         len = ihead+ifmt+2
         nperli = 80/len
c     
         if (len .le. 9 ) then
            assign 102 to ix
         elseif (ifmt .le. 9) then
            assign 103 to ix
         else 
            assign 104 to ix
         endif
c     
         write(valfmt,ix) nperli,len,ifmt
 102     format(1h(,i2,1hF,i1,1h.,i1,1h) )
 103     format(1h(,i2,1hF,i2,1h.,i1,1h) )
 104     format(1h(,i2,1hF,i2,1h.,i2,1h) )
C
      else
         len = ifmt + 6
         nperli = 80/len
c     try to minimize the blanks in the format strings.
         if (nperli .le. 9) then
	    if (len .le. 9 ) then
	       assign 105 to ix
	    elseif (ifmt .le. 9) then
	       assign 106 to ix
	    else 
	       assign 107 to ix
	    endif
	 else 
	    if (len .le. 9 ) then
	       assign 108 to ix
	    elseif (ifmt .le. 9) then
	       assign 109 to ix
	    else 
               assign 110 to ix
            endif
         endif
c-----------
         write(valfmt,ix) nperli,len,ifmt
 105     format(1h(,i1,1hD,i1,1h.,i1,1h) )
 106     format(1h(,i1,1hD,i2,1h.,i1,1h) )
 107     format(1h(,i1,1hD,i2,1h.,i2,1h) )
 108     format(1h(,i2,1hD,i1,1h.,i1,1h) )
 109     format(1h(,i2,1hD,i2,1h.,i1,1h) )
 110     format(1h(,i2,1hD,i2,1h.,i2,1h) )
c     
      endif 	    
      valcrd = (nnz-1)/nperli+1
      nrhs   = job -2
      if (nrhs .ge. 1) then
         i = (nrhs*nrow-1)/nperli+1
         rhscrd = i
         if (guesol(1:1) .eq. 'G' .or. guesol(1:1) .eq. 'g')
     +      rhscrd = rhscrd+i
         if (guesol(2:2) .eq. 'X' .or. guesol(2:2) .eq. 'x')
     +      rhscrd = rhscrd+i
         rhstyp = 'F'//guesol
      endif 
 20   continue
c     
      totcrd = ptrcrd+indcrd+valcrd+rhscrd
c     write 4-line or five line header
      write(iounit,10) title,key,totcrd,ptrcrd,indcrd,valcrd,
     1     rhscrd,type,nrow,ncol,nnz,nrhs,ptrfmt,indfmt,valfmt,valfmt
c-----------------------------------------------------------------------
      nrwindx = 0
      if (nrhs .ge. 1) write (iounit,11) rhstyp, nrhs, nrwindx
 10   format (a72, a8 / 5i14 / a3, 11x, 4i14 / 2a16, 2a20)
 11   format(A3,11x,i14,i14)
c     
      write(iounit,ptrfmt) (ia (i), i = 1, ncol+1)
      write(iounit,indfmt) (ja (i), i = 1, nnz)
      if (job .le. 1) return
      write(iounit,valfmt) (a(i), i = 1, nnz)
      if (job .le. 2) return 
      len = nrow*nrhs 
      next = 1
      iend = len
      write(iounit,valfmt) (rhs(i), i = next, iend)
c     
c     write initial guesses if available
c     
      if (guesol(1:1) .eq. 'G' .or. guesol(1:1) .eq. 'g') then
         next = next+len
         iend = iend+ len
         write(iounit,valfmt) (rhs(i), i = next, iend)
      endif
c     
c     write exact solutions if available
c     
      if (guesol(2:2) .eq. 'X' .or. guesol(2:2) .eq. 'x') then
         next = next+len
         iend = iend+ len
         write(iounit,valfmt) (rhs(i), i = next, iend)
      endif
c     
      return
c----------end of prtmt ------------------------------------------------ 
c-----------------------------------------------------------------------
      end
c----------------------------------------------------------------------- 
      subroutine dump (i1,i2,values,a,ja,ia,iout)
      integer i1, i2, ia(*), ja(*), iout
      real*8 a(*) 
      logical values 
c-----------------------------------------------------------------------
c outputs rows i1 through i2 of a sparse matrix stored in CSR format 
c (or columns i1 through i2 of a matrix stored in CSC format) in a file, 
c one (column) row at a time in a nice readable format. 
c This is a simple routine which is useful for debugging. 
c-----------------------------------------------------------------------
c on entry:
c---------
c i1    = first row (column) to print out
c i2    = last row (column) to print out 
c values= logical. indicates whether or not to print real values.
c         if value = .false. only the pattern will be output.
c a,
c ja, 
c ia    =  matrix in CSR format (or CSC format) 
c iout  = logical unit number for output.
c---------- 
c the output file iout will have written in it the rows or columns 
c of the matrix in one of two possible formats (depending on the max 
c number of elements per row. The values are output with only 
c two digits of accuracy (D9.2). )
c-----------------------------------------------------------------------
c     local variables
      integer maxr, i, k1, k2 
c
c select mode horizontal or vertical 
c
        maxr = 0
        do 1 i=i1, i2
           maxr = max0(maxr,ia(i+1)-ia(i))
 1      continue
        
        if (maxr .le. 8) then
c
c able to do one row acros line
c
        do 2 i=i1, i2
           write(iout,100) i
	   k1=ia(i)
	   k2 = ia(i+1)-1
	   write (iout,101) (ja(k),k=k1,k2)
	   if (values) write (iout,102) (a(k),k=k1,k2)
 2      continue
      else 
c
c unable to one row acros line. do three items at a time
c across a line 
         do 3 i=i1, i2
            if (values) then
               write(iout,200) i
            else
               write(iout,203) i               
            endif
            k1=ia(i)
            k2 = ia(i+1)-1
            if (values) then
               write (iout,201) (ja(k),a(k),k=k1,k2)
            else
                write (iout,202) (ja(k),k=k1,k2)
             endif
 3       continue
      endif 
c
c formats :
c
 100  format (1h ,34(1h-),' row',i6,1x,34(1h-) )
 101  format(' col:',8(i5,6h     : ))
 102  format(' val:',8(D9.2,2h :) )
 200  format (1h ,30(1h-),' row',i3,1x,30(1h-),/
     *     3('  columns :    values  * ') )
c-------------xiiiiiihhhhhhddddddddd-*-
 201  format(3(1h ,i6,6h   :  ,D9.2,3h * ) )
 202  format(6(1h ,i5,6h  *    ) ) 
 203  format (1h ,30(1h-),' row',i3,1x,30(1h-),/
     *     3('  column  :  column   *') )
      return
c----end-of-dump--------------------------------------------------------
c-----------------------------------------------------------------------
      end
c-----------------------------------------------------------------------
      subroutine pspltm(nrow,ncol,mode,ja,ia,title,ptitle,size,munt,
     *     nlines,lines,iunt)
c-----------------------------------------------------------------------
      integer nrow,ncol,ptitle,mode,iunt, ja(*), ia(*), lines(nlines) 
      real size
      character title*(*), munt*2 
c----------------------------------------------------------------------- 
c PSPLTM - PostScript PLoTer of a (sparse) Matrix
c This version by loris renggli (renggli@masg1.epfl.ch), Dec 1991
c and Youcef Saad 
c------
c Loris RENGGLI, Swiss Federal Institute of Technology, Math. Dept
c CH-1015 Lausanne (Switzerland)  -- e-mail:  renggli@masg1.epfl.ch
c Modified by Youcef Saad -- June 24, 1992 to add a few features:
c separation lines + acceptance of MSR format.
c-----------------------------------------------------------------------
c input arguments description :
c
c nrow   = number of rows in matrix
c
c ncol   = number of columns in matrix 
c
c mode   = integer indicating whether the matrix is stored in 
c           CSR mode (mode=0) or CSC mode (mode=1) or MSR mode (mode=2) 
c
c ja     = column indices of nonzero elements when matrix is
c          stored rowise. Row indices if stores column-wise.
c ia     = integer array of containing the pointers to the 
c          beginning of the columns in arrays a, ja.
c
c title  = character*(*). a title of arbitrary length to be printed 
c          as a caption to the figure. Can be a blank character if no
c          caption is desired.
c
c ptitle = position of title; 0 under the drawing, else above
c
c size   = size of the drawing  
c
c munt   = units used for size : 'cm' or 'in'
c
c nlines = number of separation lines to draw for showing a partionning
c          of the matrix. enter zero if no partition lines are wanted.
c
c lines  = integer array of length nlines containing the coordinates of 
c          the desired partition lines . The partitioning is symmetric: 
c          a horizontal line across the matrix will be drawn in 
c          between rows lines(i) and lines(i)+1 for i=1, 2, ..., nlines
c          an a vertical line will be similarly drawn between columns
c          lines(i) and lines(i)+1 for i=1,2,...,nlines 
c
c iunt   = logical unit number where to write the matrix into.
c----------------------------------------------------------------------- 
c additional note: use of 'cm' assumes european format for paper size
c (21cm wide) and use of 'in' assumes american format (8.5in wide).
c The correct centering of the figure depends on the proper choice. Y.S.
c-----------------------------------------------------------------------
c external 
      integer LENSTR
      external LENSTR
c local variables ---------------------------------------------------
      integer n,nr,nc,maxdim,istart,ilast,ii,k,ltit
      real lrmrgn,botmrgn,xtit,ytit,ytitof,fnstit,siz
      real xl,xr, yb,yt, scfct,u2dot,frlw,delt,paperx,conv,xx,yy
      logical square 
c change square to .true. if you prefer a square frame around
c a rectangular matrix
      data haf /0.5/, zero/0.0/, conv/2.54/,square/.false./
c-----------------------------------------------------------------------
      siz = size
      nr = nrow
      nc = ncol
      n = nc
      if (mode .eq. 0) n = nr
c      nnz = ia(n+1) - ia(1) 
      maxdim = max(nrow, ncol)
      m = 1 + maxdim
      nc = nc+1
      nr = nr+1
c
c units (cm or in) to dot conversion factor and paper size
c 
      if (munt.eq.'cm' .or. munt.eq.'CM') then
         u2dot = 72.0/conv
        paperx = 21.0
      else
        u2dot = 72.0
        paperx = 8.5*conv
        siz = siz*conv
      end if
c
c left and right margins (drawing is centered)
c 
      lrmrgn = (paperx-siz)/2.0
c
c bottom margin : 2 cm
c
      botmrgn = 2.0
c scaling factor
      scfct = siz*u2dot/m
c matrix frame line witdh
      frlw = 0.25
c font size for title (cm)
      fnstit = 0.5
      ltit = LENSTR(title)
c position of title : centered horizontally
c                     at 1.0 cm vertically over the drawing
      ytitof = 1.0
      xtit = paperx/2.0
      ytit = botmrgn+siz*nr/m + ytitof
c almost exact bounding box
      xl = lrmrgn*u2dot - scfct*frlw/2
      xr = (lrmrgn+siz)*u2dot + scfct*frlw/2
      yb = botmrgn*u2dot - scfct*frlw/2
      yt = (botmrgn+siz*nr/m)*u2dot + scfct*frlw/2
      if (ltit.gt.0) then
        yt = yt + (ytitof+fnstit*0.70)*u2dot
      end if
c add some room to bounding box
      delt = 10.0
      xl = xl-delt
      xr = xr+delt
      yb = yb-delt
      yt = yt+delt
c
c correction for title under the drawing
      if (ptitle.eq.0 .and. ltit.gt.0) then
      ytit = botmrgn + fnstit*0.3
      botmrgn = botmrgn + ytitof + fnstit*0.7
      end if
c begin of output
c
      write(iunt,10) '%!'
      write(iunt,10) '%%Creator: PSPLTM routine'
      write(iunt,12) '%%BoundingBox:',xl,yb,xr,yt
      write(iunt,10) '%%EndComments'
      write(iunt,10) '/cm {72 mul 2.54 div} def'
      write(iunt,10) '/mc {72 div 2.54 mul} def'
      write(iunt,10) '/pnum { 72 div 2.54 mul 20 string'
      write(iunt,10) 'cvs print ( ) print} def'
      write(iunt,10)
     1  '/Cshow {dup stringwidth pop -2 div 0 rmoveto show} def'
c
c we leave margins etc. in cm so it is easy to modify them if
c needed by editing the output file
      write(iunt,10) 'gsave'
      if (ltit.gt.0) then
      write(iunt,*) '/Helvetica findfont ',fnstit,
     &             ' cm scalefont setfont '
      write(iunt,*) xtit,' cm ',ytit,' cm moveto '
      write(iunt,'(3A)') '(',title(1:ltit),') Cshow'
      end if
      write(iunt,*) lrmrgn,' cm ',botmrgn,' cm translate'
      write(iunt,*) siz,' cm ',m,' div dup scale '
c------- 
c draw a frame around the matrix
      write(iunt,*) frlw,' setlinewidth'
      write(iunt,10) 'newpath'
      write(iunt,11) 0, 0, ' moveto'
      if (square) then
      write(iunt,11) m,0,' lineto'
      write(iunt,11) m, m, ' lineto'
      write(iunt,11) 0,m,' lineto'
      else
      write(iunt,11) nc,0,' lineto'
      write(iunt,11) nc,nr,' lineto'
      write(iunt,11) 0,nr,' lineto'
      end if
      write(iunt,10) 'closepath stroke'
c
c     drawing the separation lines 
c 
      write(iunt,*)  ' 0.2 setlinewidth'
      do 22 kol=1, nlines 
         isep = lines(kol) 
c
c     horizontal lines 
c
         yy =  real(nrow-isep) + haf 
         xx = real(ncol+1) 
         write(iunt,13) zero, yy, ' moveto '
         write(iunt,13)  xx, yy, ' lineto stroke '
c
c vertical lines 
c
         xx = real(isep) + haf 
         yy = real(nrow+1)  
         write(iunt,13) xx, zero,' moveto '
         write(iunt,13) xx, yy, ' lineto stroke '             
 22     continue
c 
c----------- plotting loop ---------------------------------------------
c
      write(iunt,10) '1 1 translate'
      write(iunt,10) '0.8 setlinewidth'
      write(iunt,10) '/p {moveto 0 -.40 rmoveto '
      write(iunt,10) '           0  .80 rlineto stroke} def'
c     
      do 1 ii=1, n
        istart = ia(ii)
        ilast  = ia(ii+1)-1 
        if (mode .eq. 1) then
          do 2 k=istart, ilast
            write(iunt,11) ii-1, nrow-ja(k), ' p'
 2        continue 
        else
          do 3 k=istart, ilast
            write(iunt,11) ja(k)-1, nrow-ii, ' p'
 3        continue          
c add diagonal element if MSR mode.
          if (mode .eq. 2) 
     *         write(iunt,11) ii-1, nrow-ii, ' p' 
c
        endif
 1    continue
c-----------------------------------------------------------------------
      write(iunt,10) 'showpage'
      return
c
 10   format (A)
 11   format (2(I6,1x),A)
 12   format (A,4(1x,F9.2))
 13   format (2(F9.2,1x),A)
c-----------------------------------------------------------------------
      end
c
      integer function lenstr(s)
c-----------------------------------------------------------------------
c return length of the string S
c-----------------------------------------------------------------------
      character*(*) s
      integer len
      intrinsic len
      integer n
c----------------------------------------------------------------------- 
      n = len(s)
10    continue
        if (s(n:n).eq.' ') then
          n = n-1
          if (n.gt.0) go to 10
        end if
      lenstr = n
c
      return
c--------end-of-pspltm--------------------------------------------------
c-----------------------------------------------------------------------
      end
c----------------------------------------------------------------------- 
      subroutine pltmt (nrow,ncol,mode,ja,ia,title,key,type,
     1     job, iounit)
c-----------------------------------------------------------------------
c this subroutine creates a 'pic' file for plotting the pattern of 
c a sparse matrix stored in general sparse format. it is not intended
c to be a means of plotting large matrices (it is very inefficient).
c It is however useful for small matrices and can be used for example
c for inserting matrix plots in a text. The size of the plot can be
c 7in x 7in or 5 in x 5in .. There is also an option for writing a 
c 3-line header in troff (see description of parameter job).
c Author: Youcef Saad - Date: Sept., 1989
c See SPARSKIT/UNSUPP/ for a version of this to produce a post-script
c file. 
c-----------------------------------------------------------------------
c nrow   = number of rows in matrix
c
c ncol	 = number of columns in matrix 
c
c mode   = integer indicating whether the matrix is stored
c          row-wise (mode = 0) or column-wise (mode=1)
c
c ja     = column indices of nonzero elements when matrix is
c	   stored rowise. Row indices if stores column-wise.
c ia     = integer array of containing the pointers to the 
c	   beginning of the columns in arrays a, ja.
c
c title  = character*71 = title of matrix test ( character a*71 ).
c key    = character*8  = key of matrix 
c type   = character*3  = type of matrix. 
c
c job    = this integer parameter allows to set a few minor 
c          options. First it tells pltmt whether or not to
c          reduce the plot. The standard size of 7in is then
c          replaced by a 5in plot. It also tells pltmt whether or
c          not to append to the pic file a few 'troff' lines that 
c          produce a centered caption includingg the title, key and 
c          types as well as the size and number of nonzero elements.
c          job =  0 : do not reduce and do not make caption.
c          job =  1 : reduce and do not make caption.
c          job = 10 : do not reduce and make caption
c          job = 11 : reduce and make caption.
c          (i.e. trailing digit for reduction, leading digit for caption)
c
c iounit = logical unit number where to write the matrix into.
c
c-----------------------------------------------------------------------
c example of usage . 
c-----------------
c In the fortran code:
c  a) read a Harwell/Boeing matrix
c          call readmt (.....)
c	   iout = 13
c  b) generate pic file:
c          call  pltmt (nrow,ncol,mode,ja,ia,title,key,type,iout)
c	   stop
c ---------
c Then in a unix environment plot the matrix by the command
c
c	pic FOR013.DAT | troff -me | lpr -Ppsx
c
c-----------------------------------------------------------------------
c notes: 1) Plots square as well as rectangular matrices.
c            (however not as much tested with rectangular matrices.)
c	  2) the dot-size is adapted according to the size of the
c            matrix.
c	  3) This is not meant at all as a way of plotting large
c            matrices. The pic file generaled will have one line for
c            each nonzero element. It is  only meant for use in
c	     such things as document poreparations etc..
c         4) The caption written will print the 71 character long
c            title. This may not be centered correctly if the
c            title has trailing blanks (a problem with Troff).
c            if you want the title centered then you can center
c            the string in title before calling pltmt. 
c       
c-----------------------------------------------------------------------
      integer ja(*), ia(*)
      character key*8,title*72,type*3
      real x, y
c-------
      n = ncol
      if (mode .eq. 0) n = nrow
      nnz = ia(n+1) - ia(1) 
      maxdim = max0 (nrow, ncol)
      xnrow = real(nrow)
      ptsize = 0.08
      hscale = (7.0 -2.0*ptsize)/real(maxdim-1) 
      vscale = hscale 
      xwid  = ptsize + real(ncol-1)*hscale + ptsize
      xht   = ptsize + real(nrow-1)*vscale + ptsize
      xshift = (7.0-xwid)/2.0
      yshift = (7.0-xht)/2.0 
c------
      if (mod(job,10) .eq. 1) then
         write (iounit,88)
      else
         write (iounit,89)
      endif
 88   format('.PS 5in',/,'.po 1.8i')
 89   format('.PS',/,'.po 0.7i')
      write(iounit,90) 
 90   format('box invisible wid 7.0 ht 7.0 with .sw at (0.0,0.0) ') 
      write(iounit,91) xwid, xht, xshift, yshift
 91   format('box wid ',f5.2,' ht ',f5.2,
     *     ' with .sw at (',f5.2,',',f5.2,')' )
c     
c     shift points slightly to account for size of dot , etc..
c     
      tiny = 0.03
      if (mod(job,10) .eq. 1) tiny = 0.05
      xshift = xshift + ptsize - tiny
      yshift = yshift + ptsize + tiny
c     
c-----------------------------------------------------------------------
c     
      ips = 8
      if (maxdim .le. 500) ips = 10
      if (maxdim .le. 300) ips = 12
      if (maxdim .le. 100) ips = 16
      if (maxdim .lt. 50) ips = 24
      write(iounit,92) ips
 92   format ('.ps ',i2)
c     
c-----------plottingloop --------------------------------------------- 
c     
      do 1 ii=1, n
         istart = ia(ii)
         ilast  = ia(ii+1)-1 
         if (mode .ne. 0) then
            x = real(ii-1)
            do 2 k=istart, ilast
               y = xnrow-real(ja(k))
               write(iounit,128) xshift+x*hscale, yshift+y*vscale
 2          continue 
         else
            y = xnrow - real(ii)
            do 3 k=istart, ilast
               x = real(ja(k)-1)
               write(iounit,128) xshift+x*hscale, yshift+y*vscale
 3          continue	    
         endif
 1    continue
c-----------------------------------------------------------------------
 128  format(7h"." at ,f6.3,1h,,f6.3,8h ljust  )
      write (iounit, 129)
 129  format('.PE')
c     quit if caption not desired. 
      if ( (job/10) .ne. 1)  return
c     
      write(iounit,127) key, type, title
      write(iounit,130) nrow,ncol,nnz
 127  format('.sp 4'/'.ll 7i'/'.ps 12'/'.po 0.7i'/'.ce 3'/,
     *     'Matrix:  ',a8,',  Type:  ',a3,/,a72)
 130  format('Dimension: ',i4,' x ',i4,',  Nonzero elements: ',i5)
      return
c----------------end-of-pltmt ------------------------------------------
c----------------------------------------------------------------------- 
      end
c-----------------------------------------------------------------------
      subroutine smms (n,first,last,mode,a,ja,ia,iout)
      integer ia(*), ja(*), n, first, last, mode, iout
      real*8 a(*)
c-----------------------------------------------------------------------
c writes a matrix in Coordinate (SMMS) format -- 
c-----------------------------------------------------------------------
c on entry:
c---------
c n     = integer = size of matrix -- number of rows (columns if matrix
c         is stored columnwise) 
c first  = first row (column) to be output. This routine will output 
c          rows (colums) first to last. 
c last   = last row (column) to be output. 
c mode   = integer giving some information about the storage of the 
c          matrix. A 3-digit decimal number. 'htu' 
c         * u = 0 means that matrix is stored row-wise 
c         * u = 1 means that matrix is stored column-wise 
c         * t = 0 indicates that the matrix is stored in CSR format 
c         * t = 1 indicates that the matrix is stored in MSR format. 
c         * h = ... to be added. 
c a,
c ja,
c ia    =  matrix in CSR or MSR format (see mode) 
c iout  = output unit number.
c
c on return:
c----------
c the output file iout will have written in it the matrix in smms
c (coordinate format)
c
c-----------------------------------------------------------------------
        logical msr, csc 
c
c determine mode ( msr or csr )
c
        msr = .false.
        csc = .false. 
	if (mod(mode,10) .eq. 1) csc = .true.
        if ( (mode/10) .eq. 1) msr = .true. 
        
        write (iout,*) n
      do 2 i=first, last 
         k1=ia(i)
         k2 = ia(i+1)-1
c     write (iout,*) ' row ', i 
         if (msr) write(iout,'(2i6,e22.14)')  i, i, a(i) 
         do 10 k=k1, k2
            if (csc) then
               write(iout,'(2i6,e22.14)')  ja(k), i, a(k)
            else
               write(iout,'(2i6,e22.14)')  i, ja(k), a(k)
            endif 
 10      continue
 2    continue
c----end-of-smms--------------------------------------------------------
c-----------------------------------------------------------------------
      end
c-----------------------------------------------------------------------
      subroutine readsm (nmax,nzmax,n,nnz,ia,ja,a,iout,ierr)
      integer nmax, nzmax, row, n, iout, i, j, k, ierr
      integer ia(nmax+1), ja(nzmax)
      real*8  a(nzmax), x
c-----------------------------------------------------------------------
c     read a matrix in coordinate format as is used in the SMMS
c     package (F. Alvarado), i.e. the row is in ascending order.
c     Outputs the matrix in CSR format.
c-----------------------------------------------------------------------
c coded by Kesheng Wu on Oct 21, 1991 with the supervision of Y. Saad
c-----------------------------------------------------------------------
c on entry:
c---------
c nmax  = the maximum size of array
c nzmax = the maximum number of nonzeros
c iout  = the I/O unit that has the data file
c
c on return:
c---------- 
c n     = integer = size of matrix
c nnz   = number of non-zero entries in the matrix
c a,
c ja, 
c ia    = matrix in CSR format
c ierr  = error code,
c         0 -- subroutine end with intended job done
c         1 -- error in I/O unit iout
c         2 -- end-of-file reached while reading n, i.e. a empty data file
c         3 -- n non-positive or too large
c         4 -- nnz is zero or larger than nzmax
c         5 -- data file is not orgnized in the order of ascending
c              row indices
c
c in case of errors:
c   n will be set to zero (0). In case the data file has more than nzmax
c   number of entries, the first nzmax entries will be read, and are not 
c   cleared on return. The total number of entry is determined.
c   Ierr is set.
c-----------------------------------------------------------------------
c
      rewind(iout)
      nnz = 0
      ia(1) = 1
      row = 1
c
      read (iout,*, err=1000, end=1010) n
      if ((n.le.0) .or. (n.gt.nmax)) goto 1020
c
 10   nnz = nnz + 1
      read (iout, *, err=1000, end=100) i, j, x

c     set the pointers when needed
      if (i.gt.row) then
         do 20 k = row+1, i
            ia(k) = nnz
 20      continue
         row = i
      else if (i.lt.row) then
         goto 1040
      endif

      ja(nnz) = j
      a (nnz) = x

      if (nnz.lt.nzmax) then
         goto 10
      else
         goto 1030
      endif

c     normal return -- end of file reached
 100  ia(row+1) = nnz
      nnz = nnz - 1
      if (nnz.eq.0) goto 1030
c
c     everything seems to be OK.
c
      ierr = 0
      return
c
c     error handling code
c
c     error in reading data entries
c
 1000 ierr = 1
      goto 2000
c
c     empty file
c
 1010 ierr  = 2
      goto 2000
c
c     problem with n
c
 1020 ierr = 3
      goto 2000
c
c     problem with nnz
c
 1030 ierr = 4
c
c     try to determine the real number of entries, in case needed
c
      if (nnz.ge.nzmax) then
 200     read(iout, *, err=210, end=210) i, j, x
         nnz = nnz + 1
         goto 200
 210     continue
      endif
      goto 2000
c
c     data entries not ordered
c
 1040 ierr = 5
 2000 n = 0
      return
c----end-of-readsm------------------------------------------------------
c-----------------------------------------------------------------------
      end
c-----------------------------------------------------------------------
      subroutine readsk (nmax,nzmax,n,nnz,a,ja,ia,iounit,ierr)
      integer nmax, nzmax, iounit, n, nnz, i, ierr
      integer ia(nmax+1), ja(nzmax) 
      real*8 a(nzmax)
c-----------------------------------------------------------------------
c Reads matrix in Compressed Saprse Row format. The data is supposed to
c appear in the following order -- n, ia, ja, a
c Only square matrices accepted. Format has following features
c (1) each number is separated by at least one space (or end-of-line), 
c (2) each array starts with a new line.
c-----------------------------------------------------------------------
c coded by Kesheng Wu on Oct 21, 1991 with supervision of Y. Saad
c-----------------------------------------------------------------------
c on entry:
c---------
c nmax 	 = max column dimension  allowed for matrix.
c nzmax	 = max number of nonzeros elements allowed. the arrays a, 
c          and ja should be of length equal to nnz (see below).
c iounit = logical unit number where to read the matrix from.
c
c on return:
c---------- 
c ia,
c ja,
c a      = matrx in CSR format
c n      = number of rows(columns) in matrix
c nnz	 = number of nonzero elements in A. This info is returned
c          even if there is not enough space in a, ja, ia, in order
c          to determine the minimum storage needed.
c ierr   = error code,
c          0 : OK;
c          1 : error when try to read the specified I/O unit.
c          2 : end-of-file reached during reading of data file.
c          3 : array size in data file is negtive or larger than nmax;
c          4 : nunmer of nonzeros in data file is negtive or larger than nzmax
c in case of errors:
c---------
c     n is set to 0 (zero), at the same time ierr is set.
c-----------------------------------------------------------------------
c     
c     read the size of the matrix
c
      rewind(iounit)
      read (iounit, *, err=1000, end=1010) n
      if ((n.le.0).or.(n.gt.nmax)) goto 1020
c     
c     read the pointer array ia(*)
c     
      read (iounit, *, err=1000, end=1010) (ia(i), i=1, n+1)
c     
c     Number of None-Zeros
c     
      nnz = ia(n+1) - 1
      if ((nnz.le.0).or.(nnz.gt.nzmax)) goto 1030
c     
c     read the column indices array
c     
      read (iounit, *, err=1000, end=1010) (ja(i), i=1, nnz)
c     
c     read the matrix elements
c     
      read (iounit, *, err=1000, end=1010) (a(i), i=1, nnz)
c     
c     normal return
c
      ierr = 0
      return
c     
c     error handling code
c     
c     error in reading I/O unit
 1000 ierr = 1
      goto 2000
c
c     EOF reached in reading
 1010 ierr =2
      goto 2000
c
c     n non-positive or too large
 1020 ierr = 3
      n = 0
      goto 2000
c
c     NNZ non-positive or too large
 1030 ierr = 4
c     
c     the real return statement
c     
 2000 n = 0
      return
c---------end of readsk ------------------------------------------------
c-----------------------------------------------------------------------
      end
c-----------------------------------------------------------------------
      subroutine skit (n, a, ja, ia, ifmt, iounit, ierr)
c-----------------------------------------------------------------------
c     Writes a matrix in Compressed Sparse Row format to an I/O unit.
c     It tryes to pack as many number as possible into lines of less than
c     80 characters. Space is inserted in between numbers for separation
c     to avoid carrying a header in the data file. This can be viewed
c     as a simplified Harwell-Boeing format. 
c-----------------------------------------------------------------------
c Modified from subroutine prtmt written by Y. Saad
c-----------------------------------------------------------------------
c on entry:
c---------
c n      = number of rows(columns) in matrix
c a      = real*8 array containing the values of the matrix stored 
c          columnwise
c ja     = integer array of the same length as a containing the column
c          indices of the corresponding matrix elements of array a.
c ia     = integer array of containing the pointers to the beginning of 
c          the row in arrays a and ja.
c ifmt   = integer specifying the format chosen for the real values
c          to be output (i.e., for a, and for rhs-guess-sol if 
c          applicable). The meaning of ifmt is as follows.
c          * if (ifmt .lt. 100) then the D descriptor is used,
c          format Dd.m, in which the length (m) of the mantissa is 
c          precisely the integer ifmt (and d = ifmt+6)
c          * if (ifmt .gt. 100) then prtmt will use the 
c          F- descriptor (format Fd.m) in which the length of the 
c          mantissa (m) is the integer mod(ifmt,100) and the length 
c          of the integer part is k=ifmt/100 (and d = k+m+2)
c          Thus  ifmt= 4   means  D10.4  +.xxxxD+ee    while
c          ifmt=104  means  F7.4   +x.xxxx
c          ifmt=205  means  F9.5   +xx.xxxxx
c          Note: formats for ja, and ia are internally computed.
c     
c iounit = logical unit number where to write the matrix into.
c     
c on return:
c----------
c ierr   = error code, 0 for normal 1 for error in writing to iounit.
c
c on error:
c--------
c     If error is encontacted when writing the matrix, the whole matrix
c     is written to the standard output.
c     ierr is set to 1.
c-----------------------------------------------------------------------
      character ptrfmt*16,indfmt*16,valfmt*20
      integer iounit, n, ifmt, len, nperli, nnz, i, ihead
      integer ja(*), ia(*), ierr
      real*8 a(*)
c--------------
c     compute pointer format
c--------------
      nnz    = ia(n+1)
      len    = int ( alog10(0.1+real(nnz))) + 2
      nnz    = nnz - 1
      nperli = 80/len
      
c      print *, ' skit entries:', n, nnz, len, nperli

      if (len .gt. 9) then
         assign 101 to ix
      else
         assign 100 to ix
      endif
      write (ptrfmt,ix) nperli,len
 100  format(1h(,i2,1HI,i1,1h) )
 101  format(1h(,i2,1HI,i2,1h) )
c----------------------------
c     compute ROW index format
c----------------------------
      len    = int ( alog10(0.1+real(n) )) + 2
      nperli = min0(80/len,nnz)
      write (indfmt,100) nperli,len
c---------------------------
c     compute value format
c---------------------------
      if (ifmt .ge. 100) then
         ihead = ifmt/100
         ifmt = ifmt-100*ihead
         len = ihead+ifmt+3
         nperli = 80/len
c     
         if (len .le. 9 ) then
            assign 102 to ix
         elseif (ifmt .le. 9) then
            assign 103 to ix
         else 
            assign 104 to ix
         endif
c     
         write(valfmt,ix) nperli,len,ifmt
 102     format(1h(,i2,1hF,i1,1h.,i1,1h) )
 103     format(1h(,i2,1hF,i2,1h.,i1,1h) )
 104     format(1h(,i2,1hF,i2,1h.,i2,1h) )
C     
      else
         len = ifmt + 7
         nperli = 80/len
c     try to minimize the blanks in the format strings.
         if (nperli .le. 9) then
	    if (len .le. 9 ) then
	       assign 105 to ix
	    elseif (ifmt .le. 9) then
	       assign 106 to ix
	    else 
	       assign 107 to ix
	    endif
	 else 
	    if (len .le. 9 ) then
	       assign 108 to ix
	    elseif (ifmt .le. 9) then
	       assign 109 to ix
	    else 
               assign 110 to ix
            endif
         endif
c-----------
         write(valfmt,ix) nperli,len,ifmt
 105     format(1h(,i1,1hD,i1,1h.,i1,1h) )
 106     format(1h(,i1,1hD,i2,1h.,i1,1h) )
 107     format(1h(,i1,1hD,i2,1h.,i2,1h) )
 108     format(1h(,i2,1hD,i1,1h.,i1,1h) )
 109     format(1h(,i2,1hD,i2,1h.,i1,1h) )
 110     format(1h(,i2,1hD,i2,1h.,i2,1h) )
c     
      endif 	    
c     
c     output the data
c     
      write(iounit, *) n
      write(iounit,ptrfmt,err=1000) (ia(i), i = 1, n+1)
      write(iounit,indfmt,err=1000) (ja(i), i = 1, nnz)
      write(iounit,valfmt,err=1000) ( a(i), i = 1, nnz)
c
c     done, if no trouble is encounted in writing data
c
      ierr = 0
      return
c     
c     if can't write the data to the I/O unit specified, should be able to
c     write everything to standard output (unit 6)
c     
 1000 write(0, *) 'Error, Can''t write data to sepcified unit',iounit
      write(0, *) 'Write the matrix into standard output instead!'
      ierr = 1
      write(6,*) n
      write(6,ptrfmt) (ia(i), i=1, n+1)
      write(6,indfmt) (ja(i), i=1, nnz)
      write(6,valfmt) ( a(i), i=1, nnz)
      return
c----------end of skit ------------------------------------------------- 
c-----------------------------------------------------------------------
      end
c-----------------------------------------------------------------------
      subroutine prtunf(n, a, ja, ia, iout, ierr)
c-----------------------------------------------------------------------
c This subroutine dumps the arrays used for storing sparse compressed row
c format in machine code, i.e. unformatted using standard FORTRAN term.
c-----------------------------------------------------------------------
c First coded by Kesheng Wu on Oct 21, 1991 under the instruction of
c Prof. Y. Saad
c-----------------------------------------------------------------------
c On entry:
c     n: the size of the matrix (matrix is n X n)
c    ia: integer array stores the stariting position of each row.
c    ja: integer array stores the column indices of each entry.
c     a: the non-zero entries of the matrix.
c  iout: the unit number opened for storing the matrix.
c On return:
c  ierr: a error, 0 if everything's OK, else 1 if error in writing data.
c On error:
c  set ierr to 1.
c  No redirection is made, since direct the machine code to the standard
c output may cause unpridictable consequences.
c-----------------------------------------------------------------------
      integer iout, n, nnz, ierr, ia(*), ja(*)
      real*8  a(*)
      nnz = ia(n+1)-ia(1) 
c
      write(unit=iout, err=1000)  n
      write(unit=iout, err=1000) (ia(k),k=1,n+1) 
      if (nnz .gt. 0) then
         write(unit=iout, err=1000) (ja(k),k=1,nnz)
         write(unit=iout, err=1000) ( a(k),k=1,nnz)
      endif 
c
      ierr = 0
      return
c
 1000 ierr = 1
      return
      end
c---------end of prtunf ------------------------------------------------
c
c-----------------------------------------------------------------------
      subroutine readunf(nmax,nzmax,n,nnz,a,ja,ia,iounit,ierr)
c-----------------------------------------------------------------------
c This subroutine reads a matix store in machine code (FORTRAN
c unformatted form). The matrix is in CSR format.
c-----------------------------------------------------------------------
c First coded by Kesheng Wu on Oct 21, 1991 under the instruction of
c Prof. Y. Saad
c-----------------------------------------------------------------------
c On entry:
c    nmax: the maximum value of matrix size.
c   nzmax: the maximum number of non-zero entries.
c  iounit: the I/O unit that opened for reading.
c On return:
c       n: the actual size of array.
c     nnz: the actual number of non-zero entries.
c ia,ja,a: the matrix in CSR format.
c    ierr: a error code, it's same as that used in reaadsk
c          0 -- OK
c          1 -- error in reading iounit
c          2 -- end-of-file reached while reading data file
c          3 -- n is non-positive or too large
c          4 -- nnz is non-positive or too large
c On error:
c     return with n set to 0 (zero). nnz is kept if it's set already,
c     in case one want to use it to determine the size of array needed
c     to hold the data.
c-----------------------------------------------------------------------
c
      integer nmax, nzmax, n, iounit, nnz, k
      integer ia(nmax+1), ja(nzmax)
      real*8  a(nzmax)
c
      rewind iounit
c      
      read (unit=iounit, err=1000, end=1010) n
      if ((n.le.0) .or. (n.gt.nmax)) goto 1020
c
      read(unit=iounit, err=1000, end=1010) (ia(k),k=1,n+1)
c
      nnz = ia(n+1) - 1
      if ((nnz.le.0) .or. (nnz.gt.nzmax)) goto 1030
c
      read(unit=iounit, err=1000, end=1010) (ja(k),k=1,nnz)
      read(unit=iounit, err=1000, end=1010) (a(k),k=1,nnz)
c
c     everything seems to be OK.
c
      ierr = 0
      return
c
c     error handling
c
 1000 ierr = 1
      goto 2000
 1010 ierr = 2
      goto 2000
 1020 ierr = 3
      goto 2000
 1030 ierr = 4
 2000 n = 0
      return
      end
c---------end of readunf ----------------------------------------------