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