← Back to branch summary

~logan/ubuntu/trusty/suitesparse/4.2.1-3ubuntu1

« back to all changes in this revision

Viewing changes to RBio/RBread_64.f

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme
  • Date: 2007-05-29 09:36:29 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20070529093629-zowquo0b7slkk6nc
Tags: 3.0.0-2
http://bugs.debian.org/426349
http://bugs.debian.org/426354
* suitesparse builds properly twice in a row
* Bug fix: "suitesparse - FTBFS: Broken build depens: libgfortran1-dev",
  thanks to Bastian Blank (Closes: #426349).
* Bug fix: "suitesparse_3.0.0-1: FTBFS: build-depends on
  libgfortran1-dev", thanks to Steve Langasek (Closes: #426354).

Show diffs side-by-side

added added

removed removed

Lines of Context:
RBio/RBread_64.f
 
1
c=======================================================================
 
2
c=== RBio/RBread_64 ====================================================
 
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) 2007, Timothy A. Davis, Univ. of Florida
 
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*8
 
40
     $      nrow, ncol, nnz, totcrd, ptrcrd, nelnz,
 
41
     $      indcrd, valcrd, mkind, skind, info
 
42
        integer*4 cmplex
 
43
        character title*72, key*8, mtype*3, ptrfmt*16, indfmt*16,
 
44
     $      valfmt*20, rhstyp*3
 
45
        logical fem
 
46
 
 
47
c       ----------------------------------------------------------------
 
48
c       read header lines 1-4
 
49
c       ----------------------------------------------------------------
 
50
 
 
51
        read (7, 10, err = 91, end = 91)
 
52
     $          title, key,
 
53
     $          totcrd, ptrcrd, indcrd, valcrd,
 
54
     $          mtype, nrow, ncol, nnz, nelnz,
 
55
     $          ptrfmt, indfmt, valfmt
 
56
10      format (a72, a8 / 4i14 / a3, 11x, 4i14 / 2a16, a20)
 
57
 
 
58
        if (nrow .lt. 0 .or. ncol .lt. 0 .or. nnz .lt. 0) then
 
59
c           error: invalid matrix dimensions
 
60
            info = -6
 
61
            return
 
62
        endif
 
63
 
 
64
c       ----------------------------------------------------------------
 
65
c       skip the Harwell/Boeing header line 5, if present
 
66
c       ----------------------------------------------------------------
 
67
 
 
68
        read (7, 20, err = 91, end = 91) rhstyp
 
69
20      format (a3)
 
70
        if ((rhstyp (1:1) .eq. 'F' .or. rhstyp (1:1) .eq. 'f' .or.
 
71
     $      rhstyp (1:1) .eq. 'M' .or. rhstyp (1:1) .eq. 'm')) then
 
72
c           This is the 5th line Harwell/Boeing format.  Ignore it.
 
73
            call mexErrMsgTxt ('Harwell/Boeing RHS ignored')
 
74
        else
 
75
c           Backspace one record, since we just read in one row of
 
76
c           the column pointers.
 
77
            backspace (unit = 7)
 
78
        endif
 
79
 
 
80
c       ----------------------------------------------------------------
 
81
c       determine if real, pattern, integer, or complex
 
82
c       ----------------------------------------------------------------
 
83
 
 
84
        if (mtype (1:1) .eq. 'R' .or. mtype (1:1) .eq. 'r') then
 
85
 
 
86
c           real
 
87
            mkind = 0
 
88
            cmplex = 0
 
89
 
 
90
        elseif (mtype (1:1) .eq. 'P' .or. mtype (1:1) .eq. 'p') then
 
91
 
 
92
c           pattern
 
93
            mkind = 1
 
94
            cmplex = 0
 
95
 
 
96
        elseif (mtype (1:1) .eq. 'C' .or. mtype (1:1) .eq. 'c') then
 
97
 
 
98
c           complex
 
99
            mkind = 2
 
100
            cmplex = 1
 
101
 
 
102
        elseif (mtype (1:1) .eq. 'I' .or. mtype (1:1) .eq. 'i') then
 
103
 
 
104
c           integer
 
105
            mkind = 3
 
106
            cmplex = 0
 
107
 
 
108
        else
 
109
 
 
110
c           error: invalid matrix type
 
111
            info = -5
 
112
            return
 
113
 
 
114
        endif
 
115
 
 
116
c       ----------------------------------------------------------------
 
117
c       determine if the upper part must be constructed
 
118
c       ----------------------------------------------------------------
 
119
 
 
120
        if (mtype (2:2) .eq. 'R' .or. mtype (2:2) .eq. 'r') then
 
121
 
 
122
c           rectangular: RRA, PRA, IRA, and CRA matrices
 
123
            skind = -1
 
124
 
 
125
        elseif (mtype (2:2) .eq. 'U' .or. mtype (2:2) .eq. 'u') then
 
126
 
 
127
c           unsymmetric: RUA, PUA, IUA, and CUA matrices
 
128
            skind = 0
 
129
 
 
130
        elseif (mtype (2:2) .eq. 'S' .or. mtype (2:2) .eq. 's') then
 
131
 
 
132
c           symmetric: RSA, PSA, ISA, and CSA matrices
 
133
            skind = 1
 
134
 
 
135
        elseif (mtype (2:2) .eq. 'H' .or. mtype (2:2) .eq. 'h') then
 
136
 
 
137
c           Hermitian: CHA (PHA, IHA, and RHA are valid, but atypical)
 
138
            skind = 2
 
139
 
 
140
        elseif (mtype (2:2) .eq. 'Z' .or. mtype (2:2) .eq. 'z') then
 
141
 
 
142
c           skew symmetric: RZA, PZA, IZA, and CZA
 
143
            skind = 3
 
144
 
 
145
        else
 
146
 
 
147
c           error: invalid matrix type
 
148
            info = -5
 
149
            return
 
150
 
 
151
        endif
 
152
 
 
153
c       ----------------------------------------------------------------
 
154
c       assembled matrices or elemental matrices (**A, **E)
 
155
c       ----------------------------------------------------------------
 
156
 
 
157
        if (mtype (3:3) .eq. 'A' .or. mtype (3:3) .eq. 'a') then
 
158
 
 
159
c           assembled - ignore nelnz
 
160
            fem = .false.
 
161
            nelnz = 0
 
162
 
 
163
        elseif (mtype (3:3) .eq. 'E' .or. mtype (3:3) .eq. 'e') then
 
164
 
 
165
c           finite-element
 
166
            fem = .true.
 
167
            continue
 
168
 
 
169
        else
 
170
 
 
171
c           error: invalid matrix type
 
172
            info = -5
 
173
            return
 
174
 
 
175
        endif
 
176
 
 
177
c       ----------------------------------------------------------------
 
178
c       assembled matrices must be square if skind is not R
 
179
c       ----------------------------------------------------------------
 
180
 
 
181
        if (.not. fem .and. skind .ne. -1 .and. nrow .ne. ncol) then
 
182
 
 
183
c           error: invalid matrix dimensions
 
184
            info = -6
 
185
            return
 
186
 
 
187
        endif
 
188
 
 
189
c       ----------------------------------------------------------------
 
190
c       matrix is valid
 
191
c       ----------------------------------------------------------------
 
192
 
 
193
        info = 0
 
194
        return
 
195
 
 
196
c       ----------------------------------------------------------------
 
197
c       error reading file
 
198
c       ----------------------------------------------------------------
 
199
 
 
200
91      info = -91
 
201
        return
 
202
        end
 
203
 
 
204
 
 
205
c-----------------------------------------------------------------------
 
206
c RBpattern: read the column pointers and row indices
 
207
c-----------------------------------------------------------------------
 
208
 
 
209
c   w and cp are both of size ncol+1 (undefined on input).
 
210
c
 
211
c   The matrix is contained in Ap, Ai, and Ax on output (undefined on
 
212
c   input).  It has nzeros explicit zero entries.  cp (1..ncol+1) are
 
213
c   the column pointers for the matrix Z that will contain all the
 
214
c   explicit zero entries.  Ax is not read in (see RBrread and
 
215
c   RBcread).
 
216
c
 
217
c   info is returned as:
 
218
c       0       ok
 
219
c       -1      invalid column pointers
 
220
c       -2      row index out of range
 
221
c       -3      duplicate entry
 
222
c       -4      entries in upper triangular part of symmetric matrix
 
223
c       -5      invalid matrix type
 
224
c       -6      invalid dimensions
 
225
c       -7      matrix contains unsorted columns
 
226
c       -91     error reading file (header) 
 
227
c       -92     error reading file (column pointers) 
 
228
c       -93     error reading file (row indices) 
 
229
c       -94     error reading file (numerical values: A, or sparse b)
 
230
c-----------------------------------------------------------------------
 
231
 
 
232
        subroutine RBpattern (ptrfmt, indfmt, nrow, ncol, nnz, skind,
 
233
     $      Ap, Ai, w, cp, info, nw)
 
234
        integer*8
 
235
     $      nrow, ncol, nnz, skind, info, Ap (ncol+1), Ai (nnz),
 
236
     $      nw, w (nw), cp (ncol+1)
 
237
        character ptrfmt*16, indfmt*16
 
238
        integer*8
 
239
     $      j, i, p, ilast
 
240
 
 
241
c       ----------------------------------------------------------------
 
242
c       read the pointers and check them
 
243
c       ----------------------------------------------------------------
 
244
 
 
245
        call RBiread (ptrfmt, ncol+1, Ap, info)
 
246
        if (info .lt. 0) then
 
247
            return
 
248
        endif
 
249
 
 
250
        if (Ap (1) .ne. 1 .or. Ap (ncol+1) - 1 .ne. nnz) then
 
251
c           error: invalid matrix (col pointers)
 
252
            info = -1
 
253
            return
 
254
        endif
 
255
        do 10 j = 2, ncol+1
 
256
            if (Ap (j) .lt. Ap (j-1)) then
 
257
c               error: invalid matrix (col pointers)
 
258
                info = -1
 
259
                return
 
260
            endif
 
261
10      continue
 
262
 
 
263
c       ----------------------------------------------------------------
 
264
c       read the row indices and check them
 
265
c       ----------------------------------------------------------------
 
266
 
 
267
        call RBiread (indfmt, nnz, Ai, info)
 
268
        if (info .lt. 0) then
 
269
            info = -93
 
270
            return
 
271
        endif
 
272
 
 
273
        do 20 i = 1, nrow
 
274
            w (i) = -1
 
275
20      continue
 
276
 
 
277
        do 40 j = 1, ncol
 
278
            ilast = 0
 
279
            do 30 p = Ap (j), Ap (j+1) - 1
 
280
                i = Ai (p)
 
281
                if (i .lt. 1 .or. i .gt. nrow) then
 
282
c                   error: row index out of range
 
283
c                   print *, 'column j, rows!', j, i, nrow
 
284
                    info = -2
 
285
                    return
 
286
                endif
 
287
                if (w (i) .eq. j) then
 
288
c                   error: duplicate entry in matrix
 
289
c                   print *, 'column j, duplicate!', j, i
 
290
                    info = -3
 
291
                    return
 
292
                endif
 
293
                w (i) = j
 
294
                if (i .lt. ilast) then
 
295
c                   error: matrix contains unsorted columns
 
296
c                   print *, 'column j, unsorted!', j, i, ilast
 
297
                    info = -7
 
298
                    return
 
299
                endif
 
300
                ilast = i
 
301
30          continue
 
302
40      continue
 
303
 
 
304
c       ----------------------------------------------------------------
 
305
c       construct new column pointers for symmetric matrices
 
306
c       ----------------------------------------------------------------
 
307
 
 
308
        if (skind .gt. 0) then
 
309
 
 
310
c           ------------------------------------------------------------
 
311
c           compute the column counts for the whole matrix
 
312
c           ------------------------------------------------------------
 
313
 
 
314
            do 50 j = 1, ncol+1
 
315
                w (j) = 0
 
316
50          continue
 
317
 
 
318
            do 70 j = 1, ncol
 
319
                do 60 p = Ap (j), Ap (j+1)-1
 
320
                    i = Ai (p)
 
321
                    if (i .eq. j) then
 
322
c                       diagonal entry, only appears as A(j,j)
 
323
                        w (j) = w (j) + 1
 
324
                    elseif (i .gt. j) then
 
325
c                       entry in lower triangular part, A(i,j) will be
 
326
c                       duplicated as A(j,i), so count it in both cols
 
327
                        w (i) = w (i) + 1
 
328
                        w (j) = w (j) + 1
 
329
                    else
 
330
c                       error: entry in upper triangular part
 
331
                        info = -4
 
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*8
 
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*8
 
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*8
 
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*8
 
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