1c=======================================================================
2c== readhb =============================================================
3c=======================================================================
4
5c-----------------------------------------------------------------------
6c UMFPACK Copyright (c) 2005-2012 by Timothy A. Davis,
7c http://www.suitesparse.com.  All Rights Reserved.
8c See ../Doc/License.txt for License.
9c-----------------------------------------------------------------------
10
11c readhb:
12c       read a sparse matrix in the Harwell/Boeing format and
13c       output a matrix in triplet format.
14c
15c usage (for example):
16c
17c       in a Unix shell:
18c       readhb < HB/arc130.rua > tmp/A
19c
20c       Then, in MATLAB, you can do the following:
21c       >> load tmp/A
22c       >> A = spconvert (A) ;
23c       >> spy (A)
24
25        integer nzmax, nmax
26        parameter (nzmax = 20000000, nmax = 250000)
27        integer Ptr (nmax), Index (nzmax), n, nz, totcrd, ptrcrd,
28     $          indcrd, valcrd, rhscrd, ncol, nrow, nrhs, row, col, p
29        character title*72, key*30, type*3, ptrfmt*16,
30     $          indfmt*16, valfmt*20, rhsfmt*20
31        logical sym
32        double precision Value (nzmax), skew
33        character rhstyp*3
34        integer nzrhs, nel
35
36        integer ne, nnz
37
38c-----------------------------------------------------------------------
39
40c       read header information from Harwell/Boeing matrix
41
42        read (5, 10, err = 998)
43     $          title, key,
44     $          totcrd, ptrcrd, indcrd, valcrd, rhscrd,
45     $          type, nrow, ncol, nz, nel,
46     $          ptrfmt, indfmt, valfmt, rhsfmt
47        if (rhscrd .gt. 0) then
48c          new Harwell/Boeing format:
49           read (5, 20, err = 998) rhstyp,nrhs,nzrhs
50           endif
5110      format (a72, a8 / 5i14 / a3, 11x, 4i14 / 2a16, 2a20)
5220      format (a3, 11x, 2i14)
53
54        skew = 0.0
55        if (type (2:2) .eq. 'Z' .or. type (2:2) .eq. 'z') skew = -1.0
56        if (type (2:2) .eq. 'S' .or. type (2:2) .eq. 's') skew =  1.0
57        sym = skew .ne. 0.0
58
59        write (0, 31) key
6031      format ('Matrix key: ', a8)
61
62        n = max (nrow, ncol)
63
64        if (n .ge. nmax .or. nz .gt. nzmax) then
65           write (0, *) 'Matrix too big!'
66           write (0, *) '(recompile readhb.f with larger nzmax, nmax)'
67           stop
68        endif
69
70        read (5, ptrfmt, err = 998) (Ptr (p), p = 1, ncol+1)
71        read (5, indfmt, err = 998) (Index (p), p = 1, nz)
72
73        do 55 col = ncol+2, n+1
74           Ptr (col) = Ptr (ncol+1)
7555      continue
76
77c       read the values
78        if (valcrd .gt. 0) then
79           read (5, valfmt, err = 998) (Value (p), p = 1, nz)
80        else
81           do 50 p = 1, nz
82              Value (p) = 1
8350         continue
84        endif
85
86c  create the triplet form of the input matrix
87
88        ne = 0
89        nnz = 0
90        do 100 col = 1, n
91           do 90 p = Ptr (col), Ptr (col+1) - 1
92              row = Index (p)
93
94              ne = ne + 1
95                  nnz = nnz + 1
96                  write (6, 200) row, col, Value (p)
97
98              if (sym .and. row .ne. col) then
99                 ne = ne + 1
100                 if (Value (p) .ne. 0) then
101                    nnz = nnz + 1
102                    write (6, 200) col, row, skew * Value (p)
103                 endif
104              endif
105
10690            continue
107100        continue
108200     format (2i7, e30.18e3)
109
110c       write (0,*) 'Number of entries: ',ne,' True nonzeros: ', nnz
111        stop
112
113998     write (0,*) 'Read error: Harwell/Boeing matrix'
114        stop
115        end
116