← Back to branch summary

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

« back to all changes in this revision

Viewing changes to RBio/RBwrite_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/RBwrite_64.f
 
1
c=======================================================================
 
2
c=== RBio/RBwrite_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 RBkind: determine the type of a MATLAB matrix
 
12
c-----------------------------------------------------------------------
 
13
c
 
14
c   input: a zero-based MATLAB sparse matrix
 
15
c
 
16
c       nrow    number of rows of A
 
17
c       ncol    number of columns of A
 
18
c       Ap      size ncol+1, column pointers
 
19
c       Ai      size nnz, row indices (nnz = Ap (ncol+1))
 
20
c       Ax      size nnz, real values
 
21
c       Az      size nnz, imaginary values (not accessed if A is real)
 
22
c       cmplex  1 if A is complex, 0 otherwise
 
23
c
 
24
c   output:
 
25
c       mkind:  r: 0 (real), p: 1 (pattern), c: 2 (complex),
 
26
c               i: 3 (integer)
 
27
c       skind:  r: -1 (rectangular), u: 0 (unsymmetric), s: 1 symmetric,
 
28
c               h: 2 (Hermitian), z: 3 (skew symmetric)
 
29
c
 
30
c   workspace:
 
31
c       munch   size ncol+1, not defined on input or output
 
32
c
 
33
c   Note that the MATLAB matrix is zero-based (Ap and Ai).  1 must be
 
34
c   added whenever they are used (see "1+" in the code below).
 
35
c
 
36
c   See also SuiteSparse/CHOLMOD/MatrixOps/cholmod_symmetry.c, which
 
37
c   also determines if the diagonal is positive.
 
38
c-----------------------------------------------------------------------
 
39
 
 
40
        subroutine RBkind (nrow, ncol, Ap, Ai, Ax, Az,
 
41
     $      cmplex, mkind, skind, mtype, nnz, munch, kmin, kmax)
 
42
 
 
43
        integer*8
 
44
     $      nrow, ncol, Ap (ncol+1), Ai (*), mkind, skind,
 
45
     $      munch (ncol+1), p, i, j, pt, nnz, k, kmin, kmax
 
46
        integer*4 cmplex
 
47
        double precision Ax (*), Az (*), x_r, x_i, xtr, xti
 
48
        logical is_p, is_s, is_h, is_z, is_int
 
49
        character mtype*3
 
50
 
 
51
c       ----------------------------------------------------------------
 
52
c       determine numeric type (I*A, R*A, P*A, C*A)
 
53
c       ----------------------------------------------------------------
 
54
 
 
55
c       pattern: if real and all entries are 1.
 
56
c       integer: if real and all entries are integers.
 
57
c       complex: if cmplex is 1.
 
58
c       real: otherwise.
 
59
 
 
60
        nnz = 1+ (Ap (ncol+1) - 1)
 
61
        kmin = 0
 
62
        kmax = 0
 
63
 
 
64
        if (cmplex .eq. 1) then
 
65
c           complex matrix (C*A)
 
66
            mtype (1:1) = 'c'
 
67
            mkind = 2
 
68
        else
 
69
c           select P** format if all entries are equal to 1
 
70
c           select I** format if all entries are integer and
 
71
c           between -99,999,999 and +999,999,999
 
72
            is_p = .true.
 
73
            is_int = .true.
 
74
            k = dint (Ax (1))
 
75
            kmin = k
 
76
            kmax = k
 
77
            do 10 p = 1, nnz
 
78
                if (Ax (p) .ne. 1) then
 
79
                    is_p = .false.
 
80
                endif
 
81
                k = dint (Ax (p))
 
82
                kmin = min (kmin, k)
 
83
                kmax = max (kmax, k)
 
84
                if (k .ne. Ax (p)) then
 
85
                    is_int = .false.
 
86
                endif
 
87
                if (k .le. -99999999 .or. k .ge. 999999999) then
 
88
c                   use real format for really big integers
 
89
                    is_int = .false.
 
90
                endif
 
91
                if (.not. is_int .and. .not. is_p) then
 
92
                    goto 20
 
93
                endif
 
94
10              continue
 
95
20          continue
 
96
            if (is_p) then
 
97
c               pattern-only matrix (P*A)
 
98
                mtype (1:1) = 'p'
 
99
                mkind = 1
 
100
            elseif (is_int) then
 
101
c               integer matrix (I*A)
 
102
                mtype (1:1) = 'i'
 
103
                mkind = 3
 
104
            else
 
105
c               real matrix (R*A)
 
106
                mtype (1:1) = 'r'
 
107
                mkind = 0
 
108
            endif
 
109
        endif
 
110
 
 
111
c       only assembled matrices are handled
 
112
        mtype (3:3) = 'a'
 
113
 
 
114
c       ----------------------------------------------------------------
 
115
c       determine symmetry (*RA, *UA, *SA, *HA, *ZA)
 
116
c       ----------------------------------------------------------------
 
117
 
 
118
c       Note that A must have sorted columns for this method to work.
 
119
c       This is not checked, since all MATLAB matrices "should" have
 
120
c       sorted columns.  Use spcheck(A) to check for this, if needed.
 
121
 
 
122
        if (nrow .ne. ncol) then
 
123
c           rectangular matrix (*RA), no need to check values or pattern
 
124
            mtype (2:2) = 'r'
 
125
            skind = -1
 
126
            return
 
127
        endif
 
128
 
 
129
c       if complex, the matrix is Hermitian until proven otherwise
 
130
        is_h = (cmplex .eq. 1)
 
131
 
 
132
c       the matrix is symmetric until proven otherwise
 
133
        is_s = .true.
 
134
 
 
135
c       a non-pattern matrix is skew symmetric until proven otherwise
 
136
        is_z = (mkind .ne. 1)
 
137
 
 
138
c       if this method returns early, the matrix is unsymmetric
 
139
        mtype (2:2) = 'u'
 
140
        skind = 0
 
141
 
 
142
c       initialize the munch pointers
 
143
        do 30 j = 1, ncol
 
144
            munch (j) = 1+ (Ap (j))
 
145
30      continue
 
146
 
 
147
        do 50 j = 1, ncol
 
148
 
 
149
c           consider all entries not yet munched in column j
 
150
            do 40 p = munch (j), 1+ (Ap (j+1)-1)
 
151
 
 
152
                i = 1+ (Ai (p))
 
153
 
 
154
                if (i .lt. j) then
 
155
c                   entry A(i,j) is unmatched, matrix is unsymmetric
 
156
                    return
 
157
                endif
 
158
 
 
159
c               get the A(j,i) entry, if it exists
 
160
                pt = munch (i)
 
161
 
 
162
c               munch the A(j,i) entry
 
163
                munch (i) = pt + 1
 
164
 
 
165
                if (pt .ge. 1+ (Ap (i+1))) then
 
166
c                   entry A(j,i) doesn't exist, matrix unsymmetric
 
167
                    return
 
168
                endif
 
169
 
 
170
                if (1+ (Ai (pt)) .ne. j) then
 
171
c                   entry A(j,i) doesn't exist, matrix unsymmetric
 
172
                    return
 
173
                endif
 
174
 
 
175
c               A(j,i) exists; check its value with A(i,j)
 
176
 
 
177
                if (cmplex .eq. 1) then
 
178
 
 
179
c                   get A(i,j)
 
180
                    x_r = Ax (p)
 
181
                    x_i = Az (p)
 
182
c                   get A(j,i)
 
183
                    xtr = Ax (pt)
 
184
                    xti = Az (pt)
 
185
                    if (x_r .ne. xtr .or. x_i .ne. xti) then
 
186
c                       the matrix cannot be *SA
 
187
                        is_s = .false.
 
188
                    endif
 
189
                    if (x_r .ne. -xtr .or. x_i .ne. -xti) then
 
190
c                       the matrix cannot be *ZA
 
191
                        is_z = .false.
 
192
                    endif
 
193
                    if (x_r .ne. xtr .or. x_i .ne. -xti) then
 
194
c                       the matrix cannot be *HA
 
195
                        is_h = .false.
 
196
                    endif
 
197
 
 
198
                else
 
199
 
 
200
c                   get A(i,j)
 
201
                    x_r = Ax (p)
 
202
c                   get A(j,i)
 
203
                    xtr = Ax (pt)
 
204
                    if (x_r .ne. xtr) then
 
205
c                       the matrix cannot be *SA
 
206
                        is_s = .false.
 
207
                    endif
 
208
                    if (x_r .ne. -xtr) then
 
209
c                       the matrix cannot be *ZA
 
210
                        is_z = .false.
 
211
                    endif
 
212
 
 
213
                endif
 
214
 
 
215
                if (.not. (is_s .or. is_z .or. is_h)) then
 
216
c                   matrix is unsymmetric; terminate the test
 
217
                    return
 
218
                endif
 
219
 
 
220
40          continue
 
221
50      continue
 
222
 
 
223
c       ----------------------------------------------------------------
 
224
c       return the symmetry
 
225
c       ----------------------------------------------------------------
 
226
 
 
227
        if (is_h) then
 
228
c           Hermitian matrix (*HA)
 
229
            mtype (2:2) = 'h'
 
230
            skind = 2
 
231
        elseif (is_s) then
 
232
c           symmetric matrix (*SA)
 
233
            mtype (2:2) = 's'
 
234
            skind = 1
 
235
        elseif (is_z) then
 
236
c           skew symmetric matrix (*ZA)
 
237
            mtype (2:2) = 'z'
 
238
            skind = 3
 
239
        endif
 
240
 
 
241
        return
 
242
        end
 
243
 
 
244
 
 
245
c-----------------------------------------------------------------------
 
246
c RBformat: determine the format required for an array of values
 
247
c-----------------------------------------------------------------------
 
248
c
 
249
c This function ensures that a sufficiently wide format is used that
 
250
c can accurately represent the data.  It also ensures that when printed,
 
251
c the numerical values all have at least one blank space between them.
 
252
c This makes it trivial for a program written in C (say) to read in a
 
253
c matrix generated by RBwrite.
 
254
 
 
255
c ww, valfmt, valn, and is_int must be defined on input.  They
 
256
c are modified on output.
 
257
c-----------------------------------------------------------------------
 
258
 
 
259
        subroutine RBformat (nnz, x, ww, valfmt, valn, is_int,
 
260
     $      kmin, kmax)
 
261
        integer*8
 
262
     $      nnz, i, ww, k, nf (18), valn, nd (9), kmin, kmax
 
263
        double precision x (nnz), e, a, b
 
264
        logical is_int
 
265
        character*20 f (18), d (9), valfmt
 
266
        character*80 s
 
267
 
 
268
c       ----------------------------------------------------------------
 
269
c       define all possible formats
 
270
c       ----------------------------------------------------------------
 
271
 
 
272
        f (1)  = '(8E9.1)             '
 
273
        f (2)  = '(8E10.2)            '
 
274
        f (3)  = '(7E11.3)            '
 
275
        f (4)  = '(6E12.4)            '
 
276
        f (5)  = '(6E13.5)            '
 
277
        f (6)  = '(5E14.6)            '
 
278
        f (7)  = '(5E15.7)            '
 
279
        f (8)  = '(5E16.8)            '
 
280
        f (9)  = '(4E17.9)            '
 
281
        f (10) = '(4E18.10)           '
 
282
        f (11) = '(4E19.11)           '
 
283
        f (12) = '(4E20.12)           '
 
284
        f (13) = '(3E21.13)           '
 
285
        f (14) = '(3E22.14)           '
 
286
        f (15) = '(3E23.15)           '
 
287
        f (16) = '(3E24.16)           '
 
288
        f (17) = '(3E25.17)           '
 
289
        f (18) = '(3E26.18)           '
 
290
 
 
291
        nf (1)  = 8
 
292
        nf (2)  = 8
 
293
        nf (3)  = 7
 
294
        nf (4)  = 6
 
295
        nf (5)  = 6
 
296
        nf (6)  = 5
 
297
        nf (7)  = 5
 
298
        nf (8)  = 5
 
299
        nf (9)  = 4
 
300
        nf (10) = 4
 
301
        nf (11) = 4
 
302
        nf (12) = 4
 
303
        nf (13) = 3
 
304
        nf (14) = 3
 
305
        nf (15) = 3
 
306
        nf (16) = 3
 
307
        nf (17) = 3
 
308
        nf (18) = 3
 
309
 
 
310
        if (is_int) then
 
311
 
 
312
c           ------------------------------------------------------------
 
313
c           use an integer format
 
314
c           ------------------------------------------------------------
 
315
 
 
316
            call RBiformat (kmin, kmax, valfmt, valn, ww)
 
317
 
 
318
        else
 
319
 
 
320
c           ------------------------------------------------------------
 
321
c           determine if the matrix has huge values or NaN's
 
322
c           ------------------------------------------------------------
 
323
 
 
324
            do 10 i = 1, nnz
 
325
                a = abs (x (i))
 
326
                if (a .ne. 0) then
 
327
                    if (a .ne. a .or. a < 1e-90 .or. a > 1e90) then
 
328
                        ww = 18
 
329
                        valfmt = '(2E30.18E3)         '
 
330
                        valn = 2
 
331
                        return
 
332
                    endif
 
333
                endif
 
334
10          continue
 
335
 
 
336
c           ------------------------------------------------------------
 
337
c           find the required precision for a real or complex matrix
 
338
c           ------------------------------------------------------------
 
339
 
 
340
            do 20 i = 1, nnz
 
341
                a = x (i)
 
342
                do 30 k = ww,18
 
343
c                   write the value to a string, read back in, and check
 
344
                    write (unit = s, fmt = f(k)) a
 
345
                    read  (unit = s, fmt = f(k)) b
 
346
                    if (a .eq. b) then
 
347
                        ww = max (ww, k)
 
348
                        goto 40
 
349
                    endif
 
350
30              continue
 
351
40              continue
 
352
                ww = max (ww, k)
 
353
20          continue
 
354
 
 
355
c           valn is the number of entries per line
 
356
            valfmt = f (ww)
 
357
            valn = nf (ww)
 
358
 
 
359
        endif
 
360
 
 
361
        return
 
362
        end
 
363
 
 
364
 
 
365
c-----------------------------------------------------------------------
 
366
c RBwrite: write portions of the matrix to the file
 
367
c-----------------------------------------------------------------------
 
368
c
 
369
c   task 0: just count the total number of entries in the matrix
 
370
c   task 1: do task 0, and also construct w and cp
 
371
c   task 2: write the row indices
 
372
c   task 3: write the numerical values
 
373
c
 
374
c   Note that the MATLAB arrays A and Z are zero-based.  "1+" is added
 
375
c   to each use of Ap, Ai, Zp, and Zi.
 
376
c-----------------------------------------------------------------------
 
377
 
 
378
        subroutine RBwrite (task, nrow, ncol, skind, cmplex, doZ, Ap,
 
379
     $      Ai, Ax, Az, Zp, Zi, mkind,
 
380
     $      indfmt, indn, valfmt, valn, nnz, w, cp)
 
381
 
 
382
        integer*8
 
383
     $      task, nrow, ncol, Ap (*), Ai (*), Zp (*), Zi (*),
 
384
     $      cp (*), w (*), nnz, znz, ibuf (80), j, i, nbuf, pa, pz,
 
385
     $      cmplex, paend, pzend, ia, iz, skind, indn, valn, p, mkind
 
386
        logical doZ
 
387
        double precision xbuf (80), xr, xi, Ax (*), Az (*)
 
388
        character valfmt*20, indfmt*20
 
389
 
 
390
c       ----------------------------------------------------------------
 
391
c       determine number of entries in Z
 
392
c       ----------------------------------------------------------------
 
393
 
 
394
        if (doZ) then
 
395
            znz = 1+ (Zp (ncol+1) - 1)
 
396
        else
 
397
            znz = 0
 
398
        endif
 
399
 
 
400
c       clear the nonzero counts
 
401
        nnz = 0
 
402
        do 10 j = 1, ncol
 
403
            w (j) = 0
 
404
10      continue
 
405
 
 
406
c       start with an empty buffer
 
407
        nbuf = 0
 
408
 
 
409
        if (znz .eq. 0) then
 
410
 
 
411
c           ------------------------------------------------------------
 
412
c           no Z present
 
413
c           ------------------------------------------------------------
 
414
 
 
415
            do 30 j = 1, ncol
 
416
 
 
417
                do 20 pa = 1+ (Ap (j)), 1+ (Ap (j+1) - 1)
 
418
 
 
419
                    i = 1+ (Ai (pa))
 
420
                    xr = Ax (pa)
 
421
                    if (cmplex .eq. 1) then
 
422
                        xi = Az (pa)
 
423
                    endif
 
424
 
 
425
                    if (skind .le. 0 .or. i .ge. j) then
 
426
 
 
427
c                       consider the (i,j) entry with value (xr,xi)
 
428
                        nnz = nnz + 1
 
429
                        if (task .eq. 1) then
 
430
c                           only determining nonzero counts
 
431
                            w (j) = w (j) + 1
 
432
                        elseif (task .eq. 2) then
 
433
c                           printing the row indices
 
434
                            call RBiprint (indfmt, ibuf, nbuf, i, indn)
 
435
                        elseif (task .eq. 3) then
 
436
c                           printing the numerical values
 
437
                            call RBxprint (valfmt, xbuf, nbuf, xr,
 
438
     $                              valn, mkind)
 
439
                            if (cmplex .eq. 1) then
 
440
                                call RBxprint (valfmt, xbuf, nbuf, xi,
 
441
     $                              valn, mkind)
 
442
                            endif
 
443
                        endif
 
444
 
 
445
                    endif
 
446
 
 
447
20              continue
 
448
30          continue
 
449
 
 
450
        else
 
451
 
 
452
c           ------------------------------------------------------------
 
453
c           symmetric, unsymmetric or rectangular matrix, with Z present
 
454
c           ------------------------------------------------------------
 
455
 
 
456
            do 40 j = 1, ncol
 
457
 
 
458
c               find the set union of A (:,j) and Z (:,j)
 
459
 
 
460
                pa = 1+ (Ap (j))
 
461
                pz = 1+ (Zp (j))
 
462
                paend = 1+ (Ap (j+1) - 1)
 
463
                pzend = 1+ (Zp (j+1) - 1)
 
464
 
 
465
c               while entries appear in A or Z
 
466
70              continue
 
467
 
 
468
c                   get the next entry from A(:,j)
 
469
                    if (pa .le. paend) then
 
470
                        ia = 1+ (Ai (pa))
 
471
                    else
 
472
                        ia = 1+ nrow
 
473
                    endif
 
474
 
 
475
c                   get the next entry from Z(:,j)
 
476
                    if (pz .le. pzend) then
 
477
                        iz = 1+ (Zi (pz))
 
478
                    else
 
479
                        iz = 1+ nrow
 
480
                    endif
 
481
 
 
482
c                   exit loop if neither entry is present
 
483
                    if (ia .gt. nrow .and. iz .gt. nrow) goto 80
 
484
 
 
485
                    if (ia .lt. iz) then
 
486
c                       get A (i,j)
 
487
                        i = ia
 
488
                        xr = Ax (pa)
 
489
                        if (cmplex .eq. 1) then
 
490
                            xi = Az (pa)
 
491
                        endif
 
492
                        pa = pa + 1
 
493
                    else if (iz .lt. ia) then
 
494
c                       get Z (i,j)
 
495
                        i = iz
 
496
                        xr = 0
 
497
                        xi = 0
 
498
                        pz = pz + 1
 
499
                    else
 
500
c                       get A (i,j), and delete its matched Z(i,j)
 
501
                        i = ia
 
502
                        xr = Ax (pa)
 
503
                        if (cmplex .eq. 1) then
 
504
                            xi = Az (pa)
 
505
                        endif
 
506
                        pa = pa + 1
 
507
                        pz = pz + 1
 
508
                    endif
 
509
 
 
510
                    if (skind .le. 0 .or. i .ge. j) then
 
511
 
 
512
c                       consider the (i,j) entry with value (xr,xi)
 
513
                        nnz = nnz + 1
 
514
                        if (task .eq. 1) then
 
515
c                           only determining nonzero counts
 
516
                            w (j) = w (j) + 1
 
517
                        elseif (task .eq. 2) then
 
518
c                           printing the row indices
 
519
                            call RBiprint (indfmt, ibuf, nbuf, i, indn)
 
520
                        elseif (task .eq. 3) then
 
521
c                           printing the numerical values
 
522
                            call RBxprint (valfmt, xbuf, nbuf, xr,
 
523
     $                              valn, mkind)
 
524
                            if (cmplex .eq. 1) then
 
525
                                call RBxprint (valfmt, xbuf, nbuf, xi,
 
526
     $                              valn, mkind)
 
527
                            endif
 
528
                        endif
 
529
 
 
530
                    endif
 
531
 
 
532
                    goto 70
 
533
 
 
534
c               end of while loop
 
535
80              continue
 
536
 
 
537
40          continue
 
538
 
 
539
        endif
 
540
 
 
541
c       ----------------------------------------------------------------
 
542
c       determine the new column pointers, or finish printing
 
543
c       ----------------------------------------------------------------
 
544
 
 
545
        if (task .eq. 1) then
 
546
 
 
547
            cp (1) = 1
 
548
            do 100 j = 2, ncol+1
 
549
                cp (j) = cp (j-1) + w (j-1)
 
550
100         continue
 
551
 
 
552
        else if (task .eq. 2) then
 
553
 
 
554
            call RBiflush (indfmt, ibuf, nbuf)
 
555
 
 
556
        elseif (task .eq. 3) then
 
557
 
 
558
            call RBxflush (valfmt, xbuf, nbuf, mkind)
 
559
 
 
560
        endif
 
561
 
 
562
        return
 
563
        end
 
564
 
 
565
 
 
566
c-----------------------------------------------------------------------
 
567
c RBiprint: print a single integer to the file, flush buffer if needed
 
568
c-----------------------------------------------------------------------
 
569
 
 
570
        subroutine RBiprint (indfmt, ibuf, nbuf, i, indn)
 
571
        character indfmt*20
 
572
        integer*8
 
573
     $      ibuf (80), nbuf, i, indn
 
574
        if (nbuf .ge. indn) then
 
575
            call RBiflush (indfmt, ibuf, nbuf)
 
576
            nbuf = 0
 
577
        endif
 
578
        nbuf = nbuf + 1
 
579
        ibuf (nbuf) = i
 
580
        return
 
581
        end
 
582
 
 
583
 
 
584
c-----------------------------------------------------------------------
 
585
c RBiflush: flush the integer buffer to the file
 
586
c-----------------------------------------------------------------------
 
587
 
 
588
        subroutine RBiflush (indfmt, ibuf, nbuf)
 
589
        character indfmt*20
 
590
        integer*8
 
591
     $      ibuf (*), nbuf, k
 
592
        write (unit = 7, fmt = indfmt, err = 999) (ibuf (k), k = 1,nbuf)
 
593
        return
 
594
999     call mexErrMsgTxt ('error writing ints')
 
595
        return
 
596
        end
 
597
 
 
598
 
 
599
c-----------------------------------------------------------------------
 
600
c RBxprint: print a single real to the file, flush the buffer if needed
 
601
c-----------------------------------------------------------------------
 
602
 
 
603
        subroutine RBxprint (valfmt, xbuf, nbuf, x, valn, mkind)
 
604
        character valfmt*20
 
605
        integer*8
 
606
     $      nbuf, valn, mkind
 
607
        double precision xbuf (80), x
 
608
        if (nbuf .ge. valn) then
 
609
            call RBxflush (valfmt, xbuf, nbuf, mkind)
 
610
            nbuf = 0
 
611
        endif
 
612
        nbuf = nbuf + 1
 
613
        xbuf (nbuf) = x
 
614
        return
 
615
        end
 
616
 
 
617
 
 
618
c-----------------------------------------------------------------------
 
619
c RBxflush: flush the real buffer to the file
 
620
c-----------------------------------------------------------------------
 
621
 
 
622
        subroutine RBxflush (valfmt, xbuf, nbuf, mkind)
 
623
        character valfmt*20
 
624
        integer*8
 
625
     $      nbuf, k, ibuf (80), mkind
 
626
        double precision xbuf (80)
 
627
        if (mkind .eq. 3) then
 
628
c           convert to integer first; valfmt is (10I8), for example
 
629
            do 10 k = 1,nbuf
 
630
                ibuf (k) = dint (xbuf (k))
 
631
10          continue
 
632
            write (unit = 7, fmt = valfmt, err = 999)
 
633
     $          (ibuf (k), k = 1,nbuf)
 
634
        else
 
635
            write (unit = 7, fmt = valfmt, err = 999)
 
636
     $          (xbuf (k), k = 1,nbuf)
 
637
        endif
 
638
        return
 
639
999     call mexErrMsgTxt ('error writing numerical values')
 
640
        return
 
641
        end
 
642
 
 
643
 
 
644
c-----------------------------------------------------------------------
 
645
c RBiformat: determine format for printing an integer
 
646
c-----------------------------------------------------------------------
 
647
 
 
648
        subroutine RBiformat (kmin, kmax, indfmt, indn, ww)
 
649
        integer*8
 
650
     $      n, indn, kmin, kmax, ww
 
651
        character*20 indfmt
 
652
 
 
653
        if (kmin .ge. 0. and. kmax .le. 9) then
 
654
            indfmt = '(40I2)              '
 
655
            ww = 2
 
656
            indn = 40
 
657
        elseif (kmin .ge. -9 .and. kmax .le. 99) then
 
658
            indfmt = '(26I3)              '
 
659
            ww = 3
 
660
            indn = 26
 
661
        elseif (kmin .ge. -99 .and. kmax .le. 999) then
 
662
            indfmt = '(20I4)              '
 
663
            ww = 4
 
664
            indn = 20
 
665
        elseif (kmin .ge. -999 .and. kmax .le. 9999) then
 
666
            indfmt = '(16I5)              '
 
667
            ww = 5
 
668
            indn = 16
 
669
        elseif (kmin .ge. -9999 .and. kmax .le. 99999) then
 
670
            indfmt = '(13I6)              '
 
671
            ww = 6
 
672
            indn = 13
 
673
        elseif (kmin .ge. -99999 .and. kmax .le. 999999) then
 
674
            indfmt = '(11I7)              '
 
675
            ww = 7
 
676
            indn = 11
 
677
        elseif (kmin .ge. -999999 .and. kmax .le. 9999999) then
 
678
            indfmt = '(10I8)              '
 
679
            ww = 8
 
680
            indn = 10
 
681
        elseif (kmin .ge. -9999999 .and. kmax .le. 99999999) then
 
682
            indfmt = '(8I9)               '
 
683
            ww = 9
 
684
            indn = 8
 
685
        elseif (kmin .ge. -99999999 .and. kmax .le. 999999999) then
 
686
            indfmt = '(8I10)               '
 
687
            ww = 10
 
688
            indn = 8
 
689
        else
 
690
            indfmt = '(5I15)               '
 
691
            ww = 15
 
692
            indn = 5
 
693
        endif
 
694
        return
 
695
        end
 
696
 
 
697
 
 
698
c-----------------------------------------------------------------------
 
699
c RBcards: determine number of cards required
 
700
c-----------------------------------------------------------------------
 
701
 
 
702
        subroutine RBcards (nitems, nperline, ncards)
 
703
        integer*8
 
704
     $      nitems, nperline, ncards
 
705
        if (nitems .eq. 0) then
 
706
            ncards = 0
 
707
        else
 
708
            ncards = ((nitems-1) / nperline) + 1
 
709
        endif
 
710
        return
 
711
        end
 
712