1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
|
c-----------------------------------------------------------------------
c Read a sparse matrix in the Harwell/Boeing format and output a matrix
c in zero-based triplet format. Only the lower triangular part of a
c symmetric matrix is provided. Does not handle skew-symmetric
c matrices.
c-----------------------------------------------------------------------
integer nzmax, nmax
parameter (nzmax = 100000000, nmax = 1000000)
integer Ptr (nmax), Index (nzmax), n, nz, totcrd, ptrcrd,
$ indcrd, valcrd, rhscrd, ncol, nrow, nrhs, row, col, p
character title*72, key*30, type*3, ptrfmt*16,
$ indfmt*16, valfmt*20, rhsfmt*20
logical sym
double precision Value (nzmax), skew
character rhstyp*3
integer nzrhs, nel, stype
c-----------------------------------------------------------------------
c read header information from Harwell/Boeing matrix
read (5, 10, err = 998)
$ title, key,
$ totcrd, ptrcrd, indcrd, valcrd, rhscrd,
$ type, nrow, ncol, nz, nel,
$ ptrfmt, indfmt, valfmt, rhsfmt
if (rhscrd .gt. 0) then
c new Harwell/Boeing format:
read (5, 20, err = 998) rhstyp,nrhs,nzrhs
endif
10 format (a72, a8 / 5i14 / a3, 11x, 4i14 / 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) title, key, type, nrow, ncol, nz
if (rhscrd .gt. 0) then
c new Harwell/Boeing format:
c write (0, 40) rhstyp,nrhs,nzrhs
endif
30 format ( /
$ ' title: ', a72 /
$ ' key: ', a8 /
$ ' type: ', a3, ' nrow: ', i14, ' ncol: ', i14 /
$ ' nz: ', i14)
40 format (' rhstyp: ', a3, ' nrhs: ', i14, ' nzrhs: ', i14)
* write (0, *) ' sym: ', sym, ' skew: ', skew
if (skew .eq. -1) then
write (0, *) 'Cannot handle skew-symmetric matrices'
stop
endif
n = max (nrow, ncol)
if (ncol .ge. nmax .or. nz .gt. nzmax) then
write (0, *) ' Matrix too big!'
stop
endif
read (5, ptrfmt, err = 998) (Ptr (p), p = 1, ncol+1)
read (5, indfmt, err = 998) (Index (p), p = 1, nz)
c read the values
if (valcrd .gt. 0) then
read (5, valfmt, err = 998) (Value (p), p = 1, nz)
endif
c create the triplet form of the input matrix
c stype = 0: unsymmetric
c stype = -1: symmetric, lower triangular part present
stype = -skew
do 100 col = 1, ncol
do 90 p = Ptr (col), Ptr (col+1) - 1
row = Index (p)
if (valcrd .gt. 0) then
write (6, 200) row-1, col-1, Value (p)
c if (sym .and. row .ne. col) then
c write (6, 200) col-1, row-1, skew * Value (p)
c endif
else
write (6, 201) row-1, col-1
endif
90 continue
100 continue
200 format (2i8, e30.18e3)
201 format (2i8, ' 1')
stop
998 write (0,*) 'Read error: Harwell/Boeing matrix'
stop
end
|
