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
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
|
// Copyright Bruno Pinçon, ESIAL-IECN, Inria CORIDA project
// <bruno.pincon@iecn.u-nancy.fr>
//
// This set of scilab 's macros provide a few sparse utilities.
//
// This software is governed by the CeCILL license under French law and
// abiding by the rules of distribution of free software. You can use,
// modify and/ or redistribute the software under the terms of the CeCILL
// license as circulated by CEA, CNRS and INRIA at the following URL
// "http://www.cecill.info".
//
// As a counterpart to the access to the source code and rights to copy,
// modify and redistribute granted by the license, users are provided only
// with a limited warranty and the software's author, the holder of the
// economic rights, and the successive licensors have only limited
// liability.
//
// In this respect, the user's attention is drawn to the risks associated
// with loading, using, modifying and/or developing or reproducing the
// software by the user in light of its specific status of free software,
// that may mean that it is complicated to manipulate, and that also
// therefore means that it is reserved for developers and experienced
// professionals having in-depth computer knowledge. Users are therefore
// encouraged to load and test the software's suitability as regards their
// requirements in conditions enabling the security of their systems and/or
// data to be ensured and, more generally, to use and operate it in the
// same conditions as regards security.
//
// The fact that you are presently reading this means that you have had
// knowledge of the CeCILL license and that you accept its terms.
function [A,description,ref,mtype] = ReadHBSparse(filename)
//
// PURPOSE
// An utility to read the Harwell-Boeing sparse matrix
// format. Currently don't work for unassembled matrix.
// Also possible rhs presents in the file are not red.
//
// ARGUMENTS
// filename : (optional) a string given the filename
// (eventually preceeding by the path), if
// filename is not given then the function
// use uigetfile to get filename
// A : the sparse matrix
// description: a string given some information about the
// matrix
// ref : a string given the reference of the matrix
// mtype : a string given the "type" of the matrix
//
// COMMENTS
// Generally the file name is of the form ref.matrixtype
// where mtype is a 3 letters word given some
// informations (already inside the file) on the matrix :
// 1st letter : R|C|P for real|complex|pattern (no values given)
// 2d letter : S|H|Z|U for symmetric|hermitian|skew symmetric|unsymmetric
// 3d letter : A|E for assembled|unassembled matrix
// (case E is not treated by this func)
// REFERENCES
// Users' Guide for the Harwell-Boeing Sparse Matrix Collection
// Iain S. Duff, Roger G. Grimes, John G. Lewis
//
// You may found this guide and numerous sparse
// matrices (in the Harwell-Boeing format) at the
// University of Florida Sparse Matrix Collection
// web site :
//
// http://www.cise.ufl.edu/research/sparse/matrices/
//
// maintained by Tim Davis <http://www.cise.ufl.edu/~davis/>
//
// AUTHOR
// Bruno Pincon <Bruno.Pincon@iecn.u-nancy.fr>
//
[lhs, rhs] = argn()
if rhs == 0 then
filename = uigetfile("*.[rc][shzu]a", title=[gettext("Choose a sparse matrix"); ...
gettext(" then click on OK ")])
elseif rhs == 1 then
if typeof(filename) ~= "string" then
error(msprintf(gettext("%s: Wrong type for input argument #%d: String expected.\n"),"ReadHBSparse",1));
end
else
error(msprintf(gettext("%s: Wrong number of input argument(s): At least %d expected.\n"),"ReadHBSparse",0));
end
unit = file("open", filename, "old")
// 1) read the 4 or 5 header lines
line1 = read(unit,1,1,"(A)")
ref = stripblanks(part(line1,73:80))
ref = convstr(ref)
description = stripblanks(part(line1,1:72))
description = convstr(description)
// normaly the line2 contains 5 numbers and if the last is zero
// the file contains no rhs but in this case in some HB file
// we have only 4 numbers. So the following is a trick to take
// into account this pb
line2 = read(unit,1,1,"(A)") // read the line 2 as a string
line2 = evstr(line2) // this string is then transform is a row vector
if length(line2) < 5 then
Rhs_in_file = %f
else
if line2(5) == 0 then
Rhs_in_file = %f
else
Rhs_in_file = %t
warning(msprintf(gettext("%s: The file contains a rhs but it will not be read !"),"ReadHBSparse"));
end
end
line3 = read(unit,1,1,"(A)")
mtype = convstr(part(line3,1:3))
Dimensions = evstr(part(line3,4:80))
if part(mtype,3)=="e" then
error(msprintf(gettext("%s: currently don''t read unassembled (elemental) sparse matrix."),"ReadHBSparse"));
end
TypeValues = part(mtype,1) // r for real, c for complex, p for pattern
m = Dimensions(1) // number of rows
n = Dimensions(2) // number of columns
nb_nz = Dimensions(3) // number of non zeros
line4 = read(unit,1,1,"(A)") // these are the formats
form1 = stripblanks(part(line4,1:16))
form2 = stripblanks(part(line4,17:32))
form3 = stripblanks(part(line4,33:52))
// for the 2 first replace Ix by Fx.0 : the read func uses only float formats
form1 = replace_Ix_by_Fx(form1)
form2 = replace_Ix_by_Fx(form2)
if Rhs_in_file then // a 5 header line to read (but ignored)
line5 = read(unit,1,1,"(A)")
end
// 2) read the datas
col_ptr = read(unit,1,n+1,form1)
ind_row = read(unit,1,nb_nz,form2)
select TypeValues
case "p" // values given
warning(msprintf(gettext("%s: No values for this matrix (only non zero pattern) : put some 1."),"ReadHBSparse"));
val = ones(1,nb_nz)
case "r" // values are real
val = read(unit,1,nb_nz,form3)
case "c" // values are complex
valc = matrix( read(unit,1,2*nb_nz,form3) , 2, nb_nz )
val = valc(1,:) + %i*(valc(2,:))
clear valc
end
file("close", unit)
// 3) form the sparse scilab matrix
// 3-1/ form the basic matrix
ind_col = ones(1,nb_nz)
for i = 2:n
ind_col(col_ptr(i):col_ptr(i+1)-1) = i
end
A = sparse([ind_row' ind_col'], val, [m n])
clear ind_row ind_col col_ptr val // to regain some memory
// 3-2/ complete the matrix depending the symetry property
MatrixSymetry = part(mtype,2)
select MatrixSymetry
case "s" // (real or complex) symmetric matrix
A = A - diag(diag(A)) + A.'
case "h" // complex hermitian matrix
A = A - diag(diag(A)) + A'
case "z" // skew symmetric matrix
A = A - A'
end
endfunction
|
