← Back to branch summary

~ubuntu-branches/ubuntu/natty/suitesparse/natty

« back to all changes in this revision

Viewing changes to RBio/RBread.f

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme
  • Date: 2006-12-22 10:16:15 UTC
  • Revision ID: james.westby@ubuntu.com-20061222101615-2ohaj8902oix2rnk
Tags: upstream-2.3.1
Import upstream version 2.3.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
RBio/RBread.f
 
1
c=======================================================================
 
2
c=== RBio/RBread =======================================================
 
3
c=======================================================================
 
4
 
 
5
c RBio: a MATLAB toolbox for reading and writing sparse matrices in
 
6
c Rutherford/Boeing format.
 
7
c Copyright (c) 2006, Timothy A. Davis, Univ. Florida.  Version 1.0.
 
8
 
 
9
 
 
10
c-----------------------------------------------------------------------
 
11
c RBheader:  read Rutherford/Boeing header lines
 
12
c-----------------------------------------------------------------------
 
13
c
 
14
c Rutherford/Boeing file type is a 3-character string:
 
15
c
 
16
c   (1) R: real, C: complex, P: pattern only, I: integer
 
17
c       mkind: R: 0, P: 1, C: 2, I: 3
 
18
c
 
19
c   (2) S: symmetric, U: unsymmetric, H: Hermitian, Z: skew symmetric,
 
20
c       R: rectangular
 
21
c       skind: R: -1, U: 0, S: 1, H: 2, Z: 3
 
22
c       
 
23
c   (3) A: assembled, E: element form
 
24
c       nelnz = 0 for A, number of elements for E
 
25
c
 
26
c       pattern matrices are given numerical values of 1 (except PZA).
 
27
c       PZA matrices have +1 in the lower triangular part and -1 in
 
28
c       the upper triangular part.
 
29
c
 
30
c   The matrix is nrow-by-ncol with nnz entries.
 
31
c   For symmetric matrices, Ai and Ax are of size 2*nnz (the upper
 
32
c       triangular part is constructed).  To skip this construction,
 
33
c       pass skind = 0 to RBpattern.
 
34
c-----------------------------------------------------------------------
 
35
 
 
36
        subroutine RBheader (title, key, mtype, nrow, ncol, nnz,
 
37
     $      ptrfmt, indfmt, valfmt, mkind, cmplex, skind, nelnz, info)
 
38
 
 
39
        integer
 
40
     $      nrow, ncol, nnz, totcrd, ptrcrd, nelnz,
 
41
     $      indcrd, valcrd, mkind, skind, cmplex, info
 
42
        character title*72, key*8, mtype*3, ptrfmt*16, indfmt*16,
 
43
     $      valfmt*20, rhstyp*3
 
44
        logical fem
 
45
 
 
46
c       ----------------------------------------------------------------
 
47
c       read header lines 1-4
 
48
c       ----------------------------------------------------------------
 
49
 
 
50
        read (7, 10, err = 91, end = 91)
 
51
     $          title, key,
 
52
     $          totcrd, ptrcrd, indcrd, valcrd,
 
53
     $          mtype, nrow, ncol, nnz, nelnz,
 
54
     $          ptrfmt, indfmt, valfmt
 
55
10      format (a72, a8 / 4i14 / a3, 11x, 4i14 / 2a16, a20)
 
56
 
 
57
        if (nrow .lt. 0 .or. ncol .lt. 0 .or. nnz .lt. 0) then
 
58
c           error: invalid matrix dimensions
 
59
            info = -6
 
60
            return
 
61
        endif
 
62
 
 
63
c       ----------------------------------------------------------------
 
64
c       skip the Harwell/Boeing header line 5, if present
 
65
c       ----------------------------------------------------------------
 
66
 
 
67
        read (7, 20, err = 91, end = 91) rhstyp
 
68
20      format (a3)
 
69
        if ((rhstyp (1:1) .eq. 'F' .or. rhstyp (1:1) .eq. 'f' .or.
 
70
     $      rhstyp (1:1) .eq. 'M' .or. rhstyp (1:1) .eq. 'm')) then
 
71
c           This is the 5th line Harwell/Boeing format.  Ignore it.
 
72
            call mexErrMsgTxt ('Harwell/Boeing RHS ignored')
 
73
        else
 
74
c           Backspace one record, since we just read in one row of
 
75
c           the column pointers.
 
76
            backspace (unit = 7)
 
77
        endif
 
78
 
 
79
c       ----------------------------------------------------------------
 
80
c       determine if real, pattern, integer, or complex
 
81
c       ----------------------------------------------------------------
 
82
 
 
83
        if (mtype (1:1) .eq. 'R' .or. mtype (1:1) .eq. 'r') then
 
84
 
 
85
c           real
 
86
            mkind = 0
 
87
            cmplex = 0
 
88
 
 
89
        elseif (mtype (1:1) .eq. 'P' .or. mtype (1:1) .eq. 'p') then
 
90
 
 
91
c           pattern
 
92
            mkind = 1
 
93
            cmplex = 0
 
94
 
 
95
        elseif (mtype (1:1) .eq. 'C' .or. mtype (1:1) .eq. 'c') then
 
96
 
 
97
c           complex
 
98
            mkind = 2
 
99
            cmplex = 1
 
100
 
 
101
        elseif (mtype (1:1) .eq. 'I' .or. mtype (1:1) .eq. 'i') then
 
102
 
 
103
c           integer
 
104
            mkind = 3
 
105
            cmplex = 0
 
106
 
 
107
        else
 
108
 
 
109
c           error: invalid matrix type
 
110
            info = -5
 
111
            return
 
112
 
 
113
        endif
 
114
 
 
115
c       ----------------------------------------------------------------
 
116
c       determine if the upper part must be constructed
 
117
c       ----------------------------------------------------------------
 
118
 
 
119
        if (mtype (2:2) .eq. 'R' .or. mtype (2:2) .eq. 'r') then
 
120
 
 
121
c           rectangular: RRA, PRA, IRA, and CRA matrices
 
122
            skind = -1
 
123
 
 
124
        elseif (mtype (2:2) .eq. 'U' .or. mtype (2:2) .eq. 'u') then
 
125
 
 
126
c           unsymmetric: RUA, PUA, IUA, and CUA matrices
 
127
            skind = 0
 
128
 
 
129
        elseif (mtype (2:2) .eq. 'S' .or. mtype (2:2) .eq. 's') then
 
130
 
 
131
c           symmetric: RSA, PSA, ISA, and CSA matrices
 
132
            skind = 1
 
133
 
 
134
        elseif (mtype (2:2) .eq. 'H' .or. mtype (2:2) .eq. 'h') then
 
135
 
 
136
c           Hermitian: CHA (PHA, IHA, and RHA are valid, but atypical)
 
137
            skind = 2
 
138
 
 
139
        elseif (mtype (2:2) .eq. 'Z' .or. mtype (2:2) .eq. 'z') then
 
140
 
 
141
c           skew symmetric: RZA, PZA, IZA, and CZA
 
142
            skind = 3
 
143
 
 
144
        else
 
145
 
 
146
c           error: invalid matrix type
 
147
            info = -5
 
148
            return
 
149
 
 
150
        endif
 
151
 
 
152
c       ----------------------------------------------------------------
 
153
c       assembled matrices or elemental matrices (**A, **E)
 
154
c       ----------------------------------------------------------------
 
155
 
 
156
        if (mtype (3:3) .eq. 'A' .or. mtype (3:3) .eq. 'a') then
 
157
 
 
158
c           assembled - ignore nelnz
 
159
            fem = .false.
 
160
            nelnz = 0
 
161
 
 
162
        elseif (mtype (3:3) .eq. 'E' .or. mtype (3:3) .eq. 'e') then
 
163
 
 
164
c           finite-element
 
165
            fem = .true.
 
166
            continue
 
167
 
 
168
        else
 
169
 
 
170
c           error: invalid matrix type
 
171
            info = -5
 
172
            return
 
173
 
 
174
        endif
 
175
 
 
176
c       ----------------------------------------------------------------
 
177
c       assembled matrices must be square if skind is not R
 
178
c       ----------------------------------------------------------------
 
179
 
 
180
        if (.not. fem .and. skind .ne. -1 .and. nrow .ne. ncol) then
 
181
 
 
182
c           error: invalid matrix dimensions
 
183
            info = -6
 
184
            return
 
185
 
 
186
        endif
 
187
 
 
188
c       ----------------------------------------------------------------
 
189
c       matrix is valid
 
190
c       ----------------------------------------------------------------
 
191
 
 
192
        info = 0
 
193
        return
 
194
 
 
195
c       ----------------------------------------------------------------
 
196
c       error reading file
 
197
c       ----------------------------------------------------------------
 
198
 
 
199
91      info = -91
 
200
        return
 
201
        end
 
202
 
 
203
 
 
204
c-----------------------------------------------------------------------
 
205
c RBpattern: read the column pointers and row indices
 
206
c-----------------------------------------------------------------------
 
207
 
 
208
c   w and cp are both of size ncol+1 (undefined on input).
 
209
c
 
210
c   The matrix is contained in Ap, Ai, and Ax on output (undefined on
 
211
c   input).  It has nzeros explicit zero entries.  cp (1..ncol+1) are
 
212
c   the column pointers for the matrix Z that will contain all the
 
213
c   explicit zero entries.  Ax is not read in (see RBrread and
 
214
c   RBcread).
 
215
c
 
216
c   info is returned as:
 
217
c       0       ok
 
218
c       -1      invalid column pointers
 
219
c       -2      row index out of range
 
220
c       -3      duplicate entry
 
221
c       -4      entries in upper triangular part of symmetric matrix
 
222
c       -5      invalid matrix type
 
223
c       -6      invalid dimensions
 
224
c       -7      matrix contains unsorted columns
 
225
c       -91     error reading file (header) 
 
226
c       -92     error reading file (column pointers) 
 
227
c       -93     error reading file (row indices) 
 
228
c       -94     error reading file (numerical values: A, or sparse b)
 
229
c-----------------------------------------------------------------------
 
230
 
 
231
        subroutine RBpattern (ptrfmt, indfmt, nrow, ncol, nnz, skind,
 
232
     $      Ap, Ai, w, cp, info, nw)
 
233
        integer
 
234
     $      nrow, ncol, nnz, skind, info, Ap (ncol+1), Ai (nnz),
 
235
     $      nw, w (nw), cp (ncol+1)
 
236
        character ptrfmt*16, indfmt*16
 
237
        integer
 
238
     $      j, i, p, ilast
 
239
 
 
240
c       ----------------------------------------------------------------
 
241
c       read the pointers and check them
 
242
c       ----------------------------------------------------------------
 
243
 
 
244
        call RBiread (ptrfmt, ncol+1, Ap, info)
 
245
        if (info .lt. 0) then
 
246
            return
 
247
        endif
 
248
 
 
249
        if (Ap (1) .ne. 1 .or. Ap (ncol+1) - 1 .ne. nnz) then
 
250
c           error: invalid matrix (col pointers)
 
251
            info = -1
 
252
            return
 
253
        endif
 
254
        do 10 j = 2, ncol+1
 
255
            if (Ap (j) .lt. Ap (j-1)) then
 
256
c               error: invalid matrix (col pointers)
 
257
                info = -1
 
258
                return
 
259
            endif
 
260
10      continue
 
261
 
 
262
c       ----------------------------------------------------------------
 
263
c       read the row indices and check them
 
264
c       ----------------------------------------------------------------
 
265
 
 
266
        call RBiread (indfmt, nnz, Ai, info)
 
267
        if (info .lt. 0) then
 
268
            info = -93
 
269
            return
 
270
        endif
 
271
 
 
272
        do 20 i = 1, nrow
 
273
            w (i) = -1
 
274
20      continue
 
275
 
 
276
        do 40 j = 1, ncol
 
277
            ilast = 0
 
278
            do 30 p = Ap (j), Ap (j+1) - 1
 
279
                i = Ai (p)
 
280
                if (i .lt. 1 .or. i .gt. nrow) then
 
281
c                   error: row index out of range
 
282
c                   print *, 'column j, rows!', j, i, nrow
 
283
                    info = -2
 
284
                    return
 
285
                endif
 
286
                if (w (i) .eq. j) then
 
287
c                   error: duplicate entry in matrix
 
288
                    print *, 'column j, duplicate!', j, i
 
289
                    info = -3
 
290
                    return
 
291
                endif
 
292
                w (i) = j
 
293
                if (i .lt. ilast) then
 
294
c                   error: matrix contains unsorted columns
 
295
                    print *, 'column j, unsorted!', j, i, ilast
 
296
                    info = -7
 
297
                    return
 
298
                endif
 
299
                ilast = i
 
300
30          continue
 
301
40      continue
 
302
 
 
303
c       ----------------------------------------------------------------
 
304
c       construct new column pointers for symmetric matrices
 
305
c       ----------------------------------------------------------------
 
306
 
 
307
        if (skind .gt. 0) then
 
308
 
 
309
c           ------------------------------------------------------------
 
310
c           compute the column counts for the whole matrix
 
311
c           ------------------------------------------------------------
 
312
 
 
313
            do 50 j = 1, ncol+1
 
314
                w (j) = 0
 
315
50          continue
 
316
 
 
317
            do 70 j = 1, ncol
 
318
                do 60 p = Ap (j), Ap (j+1)-1
 
319
                    i = Ai (p)
 
320
                    if (i .eq. j) then
 
321
c                       diagonal entry, only appears as A(j,j)
 
322
                        w (j) = w (j) + 1
 
323
                    elseif (i .gt. j) then
 
324
c                       entry in lower triangular part, A(i,j) will be
 
325
c                       duplicated as A(j,i), so count it in both cols
 
326
                        w (i) = w (i) + 1
 
327
                        w (j) = w (j) + 1
 
328
                    else
 
329
c                       error: entry in upper triangular part
 
330
                        info = -4
 
331
c                       print *, 'column j, entries in upper!',j
 
332
                        return
 
333
                    endif
 
334
60              continue
 
335
70          continue
 
336
 
 
337
c           ------------------------------------------------------------
 
338
c           compute the new column pointers
 
339
c           ------------------------------------------------------------
 
340
 
 
341
            cp (1) = 1
 
342
            do 80 j = 2, ncol+1
 
343
                cp (j) = cp (j-1) + w (j-1)
 
344
80          continue
 
345
 
 
346
        endif
 
347
 
 
348
c       ----------------------------------------------------------------
 
349
c       matrix is valid
 
350
c       ----------------------------------------------------------------
 
351
 
 
352
        info = 0
 
353
        return
 
354
        end
 
355
 
 
356
 
 
357
c-----------------------------------------------------------------------
 
358
c RBmangle: convert 1-based matrix into 0-based
 
359
c-----------------------------------------------------------------------
 
360
 
 
361
        subroutine RBmangle (ncol, Ap, Ai, nnz)
 
362
        integer
 
363
     $      ncol, nnz, p, j, Ap (ncol+1), Ai (*)
 
364
        nnz = Ap (ncol + 1) - 1
 
365
 
 
366
        do 10 j = 1, ncol+1
 
367
            Ap (j) = Ap (j) - 1
 
368
10      continue
 
369
        do 20 p = 1, nnz
 
370
            Ai (p) = Ai (p) - 1
 
371
20      continue
 
372
 
 
373
        return
 
374
        end
 
375
 
 
376
c-----------------------------------------------------------------------
 
377
 
 
378
        subroutine RBiread (ifmt, n, I, info)
 
379
        integer
 
380
     $      n, I (n), info, p
 
381
        character ifmt*16
 
382
        info = 0
 
383
        read (7, ifmt, err = 92, end = 92) (I (p), p = 1, n)
 
384
        return
 
385
92      info = -92
 
386
        return
 
387
        end
 
388
 
 
389
c-----------------------------------------------------------------------
 
390
 
 
391
        subroutine RBxread (xfmt, n, X, info)
 
392
        integer
 
393
     $      mkind, n, info, p
 
394
        double precision X (n)
 
395
        character xfmt*20
 
396
        info = 0
 
397
        read (7, xfmt, err = 94, end = 94) (X (p), p = 1, n)
 
398
        return
 
399
94      info = -94
 
400
        return
 
401
        end
 
402
 
 
403
 
 
404
c-----------------------------------------------------------------------
 
405
c RBerr: report an error to MATLAB
 
406
c-----------------------------------------------------------------------
 
407
c
 
408
c   info = 0 is OK, info < 0 is a fatal error, info > 0 is a warning
 
409
 
 
410
        subroutine RBerr (info)
 
411
        integer
 
412
     $      info
 
413
 
 
414
        if (info .eq. -7) then
 
415
            call mexErrMsgTxt ('matrix contains unsorted columns')
 
416
 
 
417
        elseif (info .eq. -1) then
 
418
            call mexErrMsgTxt ('invalid matrix (col pointers)')
 
419
 
 
420
        elseif (info .eq. -2) then
 
421
            call mexErrMsgTxt ('row index out of range)')
 
422
 
 
423
        elseif (info .eq. -3) then
 
424
            call mexErrMsgTxt ('duplicate entry in matrix')
 
425
 
 
426
        elseif (info .eq. -4) then
 
427
            call mexErrMsgTxt ('invalid symmetric matrix')
 
428
 
 
429
        elseif (info .eq. -5) then
 
430
            call mexErrMsgTxt ('invalid matrix type')
 
431
 
 
432
        elseif (info .eq. -6) then
 
433
            call mexErrMsgTxt ('invalid matrix dimensions')
 
434
 
 
435
        elseif (info .eq. -911) then
 
436
            call mexErrMsgTxt ('finite-element form not supported')
 
437
 
 
438
        elseif (info .eq. -91) then
 
439
            call mexErrMsgTxt ('error reading file (header)')
 
440
 
 
441
        elseif (info .eq. -92) then
 
442
            call mexErrMsgTxt ('error reading file (column pointers)')
 
443
 
 
444
        elseif (info .eq. -93) then
 
445
            call mexErrMsgTxt ('error reading file (row indices)')
 
446
 
 
447
        elseif (info .eq. -94) then
 
448
            call mexErrMsgTxt ('error reading file (numerical values)')
 
449
 
 
450
        elseif (info .eq. -95) then
 
451
            call mexErrMsgTxt ('error reading file (right-hand-side)')
 
452
 
 
453
        elseif (info .lt. 0) then
 
454
            print *, info
 
455
            call mexErrMsgTxt ('error (unspecified)')
 
456
 
 
457
        elseif (info .gt. 0) then
 
458
            print *, info
 
459
            call mexErrMsgTxt ('warning (unspecified)')
 
460
 
 
461
        endif
 
462
        return
 
463
        end
 
464