← Back to branch summary

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

« back to all changes in this revision

Viewing changes to RBio/RBcread.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/RBcread.f
1
 
c=======================================================================
2
 
c=== RBio/RBcread ======================================================
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 RB*read:  read a Rutherford/Boeing matrix
12
 
c-----------------------------------------------------------------------
13
 
 
14
 
        subroutine RBcread
15
 
     $      (nrow, ncol, nnz, ptrfmt, indfmt, valfmt,
16
 
     $      mkind, skind, Ap, Ai, Ax, nzeros, w, cp, info, nw, nnz1)
17
 
        integer
18
 
     $      nrow, ncol, nnz, mkind, skind, p, j, i, alen, llen, k,
19
 
     $      ilast, nzeros, info, nw, nnz1
20
 
        integer
21
 
     $      Ap (ncol+1), Ai (nnz1), w (nw), cp (ncol+1)
22
 
        complex*16 Ax (nnz1), x
23
 
        character ptrfmt*16, indfmt*16, valfmt*20
24
 
 
25
 
c       ----------------------------------------------------------------
26
 
c       read the column pointers and row indices
27
 
c       ----------------------------------------------------------------
28
 
 
29
 
        call RBpattern (ptrfmt, indfmt, nrow, ncol, nnz, skind,
30
 
     $      Ap, Ai, w, cp, info, nw)
31
 
        if (info .ne. 0) then
32
 
c           error: pattern is invalid
33
 
            return
34
 
        endif
35
 
 
36
 
c       ----------------------------------------------------------------
37
 
c       read the values
38
 
c       ----------------------------------------------------------------
39
 
 
40
 
c           read nnz values
41
 
            read (7, valfmt, err = 94, end = 94) (Ax (p), p = 1, nnz)
42
 
 
43
 
c       ----------------------------------------------------------------
44
 
c       construct upper triangular part for symmetric matrices
45
 
c       ----------------------------------------------------------------
46
 
 
47
 
c       If skind is zero, then the upper part is not constructed.  This
48
 
c       allows the caller to skip this part, and create the upper part
49
 
c       of a symmetric (S,H,Z) matrix.  Just pass skind = 0.
50
 
 
51
 
        if (skind .gt. 0) then
52
 
 
53
 
c           ------------------------------------------------------------
54
 
c           shift the matrix by adding gaps to the top of each column
55
 
c           ------------------------------------------------------------
56
 
 
57
 
            do 30 j = ncol, 1, -1
58
 
 
59
 
c               number of entries in lower tri. part (incl. diagonal)
60
 
                llen = Ap (j+1) - Ap (j)
61
 
 
62
 
c               number of entries in entire column
63
 
                alen = cp (j+1) - cp (j)
64
 
 
65
 
c               move the column from Ai (Ap(j) ... Ap(j+1)-1)
66
 
c               down to Ai (cp(j+1)-llen ... cp(j+1)-1), leaving a gap
67
 
c               at Ai (Ap(j) ... cp(j+1)-llen)
68
 
 
69
 
                do 20 k = 1, llen
70
 
                    Ai (cp (j+1) - k) = Ai (Ap (j+1) - k)
71
 
                    Ax (cp (j+1) - k) = Ax (Ap (j+1) - k)
72
 
20              continue
73
 
30          continue
74
 
 
75
 
c           ------------------------------------------------------------
76
 
c           populate the upper triangular part
77
 
c           ------------------------------------------------------------
78
 
 
79
 
c           create temporary column pointers to point to the gaps
80
 
            do 40 j = 1, ncol
81
 
                w (j) = cp (j)
82
 
40          continue
83
 
 
84
 
            do 60 j = 1, ncol
85
 
 
86
 
c               scan the entries in the lower tri. part, in
87
 
c               Ai (cp(j+1)-llen ... cp(j+1)-1)
88
 
                llen = Ap (j+1) - Ap (j)
89
 
                do 50 k = 1, llen
90
 
 
91
 
c                   get the A(i,j) entry in the lower triangular part
92
 
                    i = Ai (cp (j+1) - k)
93
 
                    x = Ax (cp (j+1) - k)
94
 
 
95
 
c                   add A(j,i) as the next entry in column i (excl diag)
96
 
                    if (i .ne. j) then
97
 
                        p = w (i)
98
 
                        w (i) = w (i) + 1
99
 
                        Ai (p) = j
100
 
 
101
 
                        if (skind .eq. 1) then
102
 
c                           *SA matrix
103
 
                            Ax (p) = x
104
 
                        elseif (skind .eq. 2) then
105
 
c                           *HA matrix
106
 
                            Ax (p) = dconjg (x)
107
 
                        else
108
 
c                           *ZA matrix
109
 
                            Ax (p) = -x
110
 
                        endif
111
 
 
112
 
                    endif
113
 
 
114
 
50              continue
115
 
60          continue
116
 
 
117
 
c           finalize the column pointers
118
 
            do 70 j = 1, ncol+1
119
 
                Ap (j) = cp (j)
120
 
70          continue
121
 
 
122
 
        endif
123
 
 
124
 
c       ----------------------------------------------------------------
125
 
c       count the number of explicit zeros
126
 
c       ----------------------------------------------------------------
127
 
 
128
 
        nzeros = 0
129
 
        do 90 j = 1, ncol
130
 
            cp (j) = nzeros + 1
131
 
            do 80 p = Ap (j), Ap (j+1)-1
132
 
                if (Ax (p) .eq. 0) then
133
 
                    nzeros = nzeros + 1
134
 
                endif
135
 
80          continue
136
 
90      continue
137
 
        cp (ncol+1) = nzeros + 1
138
 
 
139
 
c       ----------------------------------------------------------------
140
 
c       matrix is valid
141
 
c       ----------------------------------------------------------------
142
 
 
143
 
        info = 0
144
 
        return
145
 
 
146
 
c       ----------------------------------------------------------------
147
 
c       error return
148
 
c       ----------------------------------------------------------------
149
 
 
150
 
94      info = -94
151
 
        return
152
 
        end
153
 
 
154
 
 
155
 
c-----------------------------------------------------------------------
156
 
c RB*zeros: extract explicit zero entries
157
 
c-----------------------------------------------------------------------
158
 
c
159
 
c   nrow-by-ncol: size of A and Z
160
 
c   cp: column pointers of Z on input
161
 
c   Ap, Ai, Ax: matrix with zeros on input, pruned on output
162
 
c   Zp, Zi, Zx: empty matrix on input, pattern of zeros on output
163
 
c
164
 
c-----------------------------------------------------------------------
165
 
 
166
 
        subroutine RBczeros
167
 
     $      (nrow, ncol, cp, Ap, Ai, Ax, Zp, Zi, Zx)
168
 
        integer
169
 
     $      nrow, ncol, Ap (ncol+1), Ai (*), Zp (ncol+1), Zi (*),
170
 
     $      cp (*), i, j, p, pa, pz, p1
171
 
        complex*16 Ax (*), x
172
 
        double precision Zx (*)
173
 
 
174
 
c       ----------------------------------------------------------------
175
 
c       copy the column pointers if Z is being constructed
176
 
c       ----------------------------------------------------------------
177
 
 
178
 
        do 10 j = 1, ncol+1
179
 
            Zp (j) = cp (j)
180
 
10      continue
181
 
 
182
 
c       ----------------------------------------------------------------
183
 
c       split the matrix
184
 
c       ----------------------------------------------------------------
185
 
 
186
 
        pa = 1
187
 
        pz = 1
188
 
        do 30 j = 1, ncol
189
 
c           save the new start of column j of A
190
 
            p1 = Ap (j)
191
 
            Ap (j) = pa
192
 
            pz = Zp (j)
193
 
c           split column j of A
194
 
            do 20 p = p1, Ap (j+1)-1
195
 
                i = Ai (p)
196
 
                x = Ax (p)
197
 
                if (x .eq. 0) then
198
 
c                   copy into Z
199
 
                    Zi (pz) = i
200
 
                    Zx (pz) = 1
201
 
                    pz = pz + 1
202
 
                else
203
 
c                   copy into A
204
 
                    Ai (pa) = i
205
 
                    Ax (pa) = x
206
 
                    pa = pa + 1
207
 
                endif
208
 
20          continue
209
 
30      continue
210
 
        Ap (ncol+1) = pa
211
 
 
212
 
        return
213
 
        end
214
 
 
215
 
 
216
 
c-----------------------------------------------------------------------
217
 
c RB*prune: discard explicit zero entries
218
 
c-----------------------------------------------------------------------
219
 
c
220
 
c   nrow-by-ncol: size of A
221
 
c   Ap, Ai, Ax: matrix with zeros on input, pruned on output
222
 
c
223
 
c-----------------------------------------------------------------------
224
 
 
225
 
        subroutine RBcprune
226
 
     $      (nrow, ncol, Ap, Ai, Ax)
227
 
        integer
228
 
     $      nrow, ncol, Ap (ncol+1), Ai (*), i, j, p, pa, pz, p1
229
 
        complex*16 Ax (*), x
230
 
 
231
 
c       ----------------------------------------------------------------
232
 
c       prune the matrix
233
 
c       ----------------------------------------------------------------
234
 
 
235
 
        pa = 1
236
 
        do 20 j = 1, ncol
237
 
c           save the new start of column j of A
238
 
            p1 = Ap (j)
239
 
            Ap (j) = pa
240
 
c           prune column j of A
241
 
            do 10 p = p1, Ap (j+1)-1
242
 
                i = Ai (p)
243
 
                x = Ax (p)
244
 
                if (x .ne. 0) then
245
 
c                   copy into A
246
 
                    Ai (pa) = i
247
 
                    Ax (pa) = x
248
 
                    pa = pa + 1
249
 
                endif
250
 
10          continue
251
 
20      continue
252
 
        Ap (ncol+1) = pa
253
 
 
254
 
        return
255
 
        end