subroutine ireadhb(fname,type,
     $                     nrows,ncols,nnz)

      
        implicit none
        integer nrows, ncols, nnz
        integer totcrd, ptrcrd,
     $		indcrd, valcrd, rhscrd, nrhs, row, col, p
        character title*72, key*30, type*3, ptrfmt*16,
     $          indfmt*16, valfmt*20, rhsfmt*20
        character fname*256
        logical sym
        double precision skew, myrand
        character rhstyp*3
        integer nzrhs, nel

c-----------------------------------------------------------------------

c       read header information from Harwell/Boeing matrix
        open (99, file=fname, err=999, status="OLD")

        read (99, 10, err = 999)
     $          title, key,
     $          totcrd, ptrcrd, indcrd, valcrd, rhscrd,
     $          type, nrows, ncols, nnz, nel
10      format (a72, a8 / 5i14 / a3, 11x, 4i14)

        write (0, 30)
     $          title, key,
     $          totcrd, ptrcrd, indcrd, valcrd, rhscrd,
     $          type, nrows, ncols, nnz, nel
30      format (
     $          ' title: ', a72 /
     $          ' key: ', a8 /
     $          ' Lines: tot: ', i14,' ptr: ',i14,' ind: ',i14 /
     $          '        val: ', i14,' rhs: ',i14 /
     $          ' type: ', a3, ' nrow: ', i14, ' ncol: ', i14 /
     $          ' nz: ', i14, ' elements: ', i14)

        close (99)
        return

999     write (0,*) 'Read error: Harwell/Boeing matrix'
        stop
        end

        subroutine dreadhb(fname,
     $                     nrows,ncols,nnz,
     $                     Ptr,Index,Value)

        implicit none
        integer nrows, ncols, nnz, fnamel
        integer Ptr(*), Index(*), totcrd, ptrcrd,
     $		indcrd, valcrd, rhscrd, nrhs, row, col, p
        character title*72, key*30, type*3, ptrfmt*16,
     $          indfmt*16, valfmt*20, rhsfmt*20
        character fname*256
        logical sym
        double precision Value (*), skew, myrand
        character rhstyp*3
        integer nzrhs, nel

c-----------------------------------------------------------------------

c       read header information from Harwell/Boeing matrix

        open (99, file=fname, err=998, status="OLD")

        read (99, 5, err = 998)
     $          title, key,
     $          totcrd, ptrcrd, indcrd, valcrd, rhscrd,
     $          type, nrows, ncols, nnz, nel
5       format (a72, a8 / 5i14 / a3, 11x, 4i14)

        read (99, 10, err = 998)
     $          ptrfmt, indfmt, valfmt, rhsfmt
        if (rhscrd .gt. 0) then
c          new Harwell/Boeing format:
           read (99, 20, err = 998) rhstyp,nrhs,nzrhs
           endif
10      format (2a16, 2a20)
20      format (a3, 11x, 2i14)

        skew = 0.0
        if (type (2:2) .eq. 'Z' .or. type (2:2) .eq. 'z') skew = -1.0
        if (type (2:2) .eq. 'S' .or. type (2:2) .eq. 's') skew =  1.0
        sym = skew .ne. 0.0

        write (0, 30)
     $          ptrfmt, indfmt, valfmt, rhsfmt
        if (rhscrd .gt. 0) then
c          new Harwell/Boeing format:
           write (0, 40) rhstyp,nrhs,nzrhs
           endif
30      format (
     $          ' ptrfmt: ', a20, ' rowfmt: ', a20, /
     $          ' valfmt: ', a20, ' rhsfmt: ', a20)
40      format (' rhstyp: ', a3, ' nrhs: ', i14, ' nzrhs: ', i14)
        write (0, *) ' sym: ', sym, ' skew: ', skew

        print *,'reading colptr'

        read (99, ptrfmt, err = 998) (Ptr (p), p = 1, ncols+1)
        print *,'reading rowind'
        read (99, indfmt, err = 998) (Index (p), p = 1, nnz)

c      what's this? maybe for rectangualr matrices

        do 55 col = ncols+2, ncols+1
           Ptr (col) = Ptr (ncols+1)
55         continue

        print *,'reading values'
c       read the values, or create random-valued matrix
        if (valcrd .gt. 0) then
           read (99, valfmt, err = 998) (Value (p), p = 1, nnz)
        else
          if (sym) then
            do 57 col = 1, ncols
              do 56 p = Ptr(col), Ptr(col+1)-1
                row = Index(p)
                if (row .eq. col) then
                  Value(p) = ncols
                else
                  Value(p) = -1.0
                endif
56            continue
57          continue
          else
            Value (1) = myrand (0)
            do 50 p = 1, nnz
               Value (p) = myrand (-1)
50          continue
          endif
        endif

c  create the triplet form of the input matrix

c        do 100 col = 1, n
c           do 90 p = Ptr (col), Ptr (col+1) - 1
c              row = Index (p)
c              write (6, 200) row, col, Value (p)
c              if (sym .and. row .ne. col) then
c		 write (6, 200) col, row, skew * Value (p)
c		 endif
c90            continue
c100        continue
c200	format (2i7, e26.16e3)

        close (99)
        return

998     write (0,*) 'Read error: Harwell/Boeing matrix'
        stop
        end

c=== Myrand ============================================================
c
c  Derived from the FA01 routines in the MUPS package (CERFACS and/or
c  Harwell).  CERFACS and/or Harwell copyrights may apply.  Permission
c  granted to use this routine in the DEMO PROGRAM only.
c
c  DEMO PROGRAM.
c
c  random number generator
c  i = 0:  reinitialize the sequence
c  i >=0:  return 0 < x < 1
c  i < 0:  return -1 < x < 1

        double precision function myrand (i)
        integer i
        double precision seed, start, mfac, d2to32
        common /mrand/ seed
        parameter (start = 1431655765.d0,
     $             d2to32 = 4294967296.d0, mfac = 9228907.d0)

        if (i .eq. 0) then
c          reinitialize to known sequence
           seed = start
           endif
        seed = dmod (seed * mfac, d2to32)

        if (i .ge. 0) then
           myrand = (seed/d2to32)
        else
           myrand = 2 * (seed/d2to32) - 1
           endif
        return
        end