← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to Lib/sparse/sparsetools/spconv.f.src

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/sparse/sparsetools/spconv.f.src
1
 
c    -*- fortran -*-
2
 
c    convert csc matrix a to csr matrix b
3
 
c    convert csr matrix a to csc matrix b
4
 
c    transpose csc matrix
5
 
c    transpose csr matrix
6
 
c    all depends on how you interpret the results
7
 
         
8
 
      subroutine <_c>transp(m,n,a,rowa,ptra,nnzamax,b,colb,ptrb)
9
 
      <_t> a(0:nnzamax-1), b(0:nnzamax-1)
10
 
      integer rowa(0:nnzamax-1), colb(0:nnzamax-1)
11
 
      integer ptra(0:n), ptrb(0:m)
12
 
      integer j, bnum, cola, ia, ra
13
 
 
14
 
      ptrb(0) = 0
15
 
      bnum = 0
16
 
      do j = 0, m-1
17
 
         do cola = 0, n-1
18
 
            do ia = ptra(cola), ptra(cola+1)-1
19
 
               ra = rowa(ia)
20
 
               if (ra.eq.j) then
21
 
                  b(bnum) = a(ia)
22
 
                  colb(bnum) = cola
23
 
                  bnum = bnum+1
24
 
               end if
25
 
            end do
26
 
         end do
27
 
         ptrb(j+1) = bnum
28
 
      end do
29
 
 
30
 
      return
31
 
      end subroutine <_c>transp
32
 
 
33
 
c      returns the index into the data array
34
 
c      and the value
35
 
      subroutine <_c>cscgetel(a,rowa,ptra,nnzamax,n,row,col,ind,val)
36
 
      
37
 
      <_t> a(0:nnzamax-1), val
38
 
      integer rowa(0:nnzamax-1), ptra(0:n)
39
 
      integer nnzamax, row, col, ind
40
 
 
41
 
      integer ia
42
 
 
43
 
      ind = -1
44
 
      val = 0.0
45
 
      do ia = ptra(col), ptra(col+1)-1
46
 
         if (rowa(ia).eq.row) then
47
 
            ind = ia
48
 
            val = a(ia)
49
 
            goto 10
50
 
         end if
51
 
      end do
52
 
 
53
 
 10   continue
54
 
      return 
55
 
 
56
 
      end subroutine <_c>cscgetel
57
 
 
58
 
      
59
 
c      set an element into the data array (assumes there is room)
60
 
c
61
 
      subroutine <_c>cscsetel(a,rowa,ptra,nnzamax,n,row,col,val)
62
 
 
63
 
      <_t> a(0:nnzamax-1), val
64
 
      integer rowa(0:nnzamax-1), ptra(0:n)
65
 
      integer nnza, n, row, col, nnzamax
66
 
 
67
 
      integer ra, newia
68
 
 
69
 
      nnza = ptra(n)
70
 
      newia = ptra(col)
71
 
      do ia = newia, ptra(col+1)-1
72
 
         ra = rowa(ia)
73
 
         if (ra.eq.row) then
74
 
c     just replace and be done
75
 
            a(ia) = val
76
 
            goto 10
77
 
         else if (ra.gt.row) then 
78
 
            newia = ia
79
 
            goto 5
80
 
         end if
81
 
      end do
82
 
 
83
 
c     we are past where it should be stored
84
 
c         (assumes sorted) so make room for this new element
85
 
c     here so that it will run if this is the first entry in the
86
 
c         column
87
 
 5    continue
88
 
      ia = newia
89
 
      do iia = nnza, ia+1, -1
90
 
         a(iia) = a(iia-1)
91
 
         rowa(iia) = rowa(iia-1)
92
 
      end do
93
 
      a(ia) = val
94
 
      rowa(ia) = row
95
 
      do iia = col+1, n
96
 
         ptra(iia) = ptra(iia)+1
97
 
      end do
98
 
 
99
 
      return
100
 
 
101
 
 10   continue
102
 
      
103
 
      return
104
 
      end subroutine <_c>cscsetel
105
 
 
106
 
c     assumes sorted by column and then row 
107
 
 
108
 
      subroutine <_c>cootocsc(n,vals,row,col,nnz,a,rowa,ptra,nnzamax,
109
 
     $     ierr)
110
 
      
111
 
      <_t> vals(0:nnz-1), a(0:nnzamax-1)
112
 
      integer rowa(0:nnzamax-1), ptra(0:n)
113
 
      integer row(0:nnz-1), col(0:nnz-1)
114
 
      integer n, nnzamax, nnz, ierr
115
 
 
116
 
      <_t> val
117
 
      integer j, k, cumsum
118
 
 
119
 
      if (nnz.gt.nnzamax) goto 999
120
 
      ierr=0
121
 
      nnza=0
122
 
      do k = 0, n-1
123
 
         ptra(k) = 0
124
 
      end do
125
 
      do ia = 0, nnzamax-1
126
 
         a(ia) = 0.0
127
 
         rowa(ia) = 0
128
 
      end do      
129
 
      do k = 0, nnz-1
130
 
         j = col(k)
131
 
         val = vals(k)
132
 
         if (val.ne.0.0) then
133
 
            a(nnza) = val
134
 
            rowa(nnza) = row(k)
135
 
            ptra(j+1) = ptra(j+1) + 1
136
 
            nnza = nnza + 1
137
 
         end if
138
 
      end do
139
 
 
140
 
c   successful completion (fix the ptr array)
141
 
      cumsum = 0
142
 
      do k = 1, n
143
 
         cumsum = cumsum + ptra(k)
144
 
         ptra(k) = cumsum         
145
 
      end do
146
 
      return      
147
 
 
148
 
      return
149
 
 
150
 
 999  continue
151
 
      ierr=1
152
 
      return
153
 
      end subroutine <_c>cootocsc
154
 
 
155
 
      subroutine <_c>csctocoo(n, vals, row, col, a, rowa, ptra, nnzamax)
156
 
      <_t> vals(0:nnzamax-1), a(0:nnzamax-1)
157
 
      integer rowa(0:nnzamax-1), ptra(0:n)
158
 
      integer row(0:nnzamax-1), col(0:nnzamax-1)
159
 
      integer n, nnzamax
160
 
 
161
 
      integer i, j, k
162
 
 
163
 
      j = 0
164
 
      do k = 0, n-1
165
 
         do i = ptra(k), ptra(k+1)-1
166
 
            row(j) = rowa(i)
167
 
            col(j) = k
168
 
            vals(j) = a(i)
169
 
            j = j + 1
170
 
         end do
171
 
      end do
172
 
 
173
 
      return
174
 
      end subroutine <_c>csctocoo
175
 
 
176
 
c     intended for re-entry in case a is too small
177
 
 
178
 
      subroutine <_c>fulltocsc(m,n,fulla,a,rowa,ptra,nnzamax,
179
 
     $     irow,jcol,ierr)
180
 
      
181
 
      <_t> fulla(0:m-1,0:n-1), a(0:nnzamax-1)
182
 
      integer rowa(0:nnzamax-1), ptra(0:n)
183
 
      integer m, n, nnzamax, irow, jcol, ierr
184
 
 
185
 
      <_t> val
186
 
      integer nnza, i, j, k, cumsum
187
 
 
188
 
      ierr=0
189
 
      nnza = ierr
190
 
      do j = jcol, n-1
191
 
         do i = irow, m-1
192
 
            val = fulla(i,j)
193
 
            if (val.ne.0.0) then
194
 
               if (nnza.ge.nnzamax) goto 999
195
 
               a(nnza) = val
196
 
               rowa(nnza) = i
197
 
               ptra(j+1) = ptra(j+1) + 1
198
 
               nnza = nnza + 1
199
 
            end if
200
 
         end do
201
 
      end do      
202
 
 
203
 
c   successful completion (fix the ptr array)
204
 
      cumsum = 0
205
 
      do k = 1, n
206
 
         cumsum = cumsum + ptra(k)
207
 
         ptra(k) = cumsum         
208
 
      end do
209
 
      return
210
 
 
211
 
 999  continue
212
 
      ierr=nnza
213
 
      irow = i
214
 
      jcol = j
215
 
      return
216
 
      end subroutine <_c>fulltocsc
217
 
 
218
 
 
219
 
      subroutine <_c>csctofull(m, n, fulla, a, rowa, ptra, nnzamax)
220
 
 
221
 
      <_t> fulla(0:m-1, 0:n-1), a(0:nnzamax-1)
222
 
      integer rowa(0:nnzamax-1), ptra(0:n)
223
 
      integer m, n, nnzamax
224
 
 
225
 
      integer ia, k
226
 
c      integer i, j
227
 
 
228
 
c      do j = 0, n-1
229
 
c         do i = 0, m-1
230
 
c            fulla(i,j) = 0.0
231
 
c         end do
232
 
c      end do
233
 
 
234
 
      do k = 0, n-1
235
 
         do ia = ptra(k), ptra(k+1)-1
236
 
            fulla(rowa(ia),k) = a(ia)
237
 
         end do
238
 
      end do
239
 
      
240
 
      return
241
 
      end subroutine <_c>csctofull
242
 
 
243
 
 
244
 
c  extract a sub-matrix from a sparse matrix
245
 
c     c = a(ibeg:iend, jbeg:jend)  (inclusive of end-points)
246
 
c
247
 
c     intended for re-entry if nnzcmax is not big enough 
248
 
c      irow and jcol should initially be 0 and 0
249
 
 
250
 
      subroutine <_c>cscextract(n,a,rowa,ptra,nnzamax,ibeg,iend,jbeg,
251
 
     $     jend, c, rowc, ptrc, nnzcmax, irow, jcol, ierr)
252
 
 
253
 
      <_t> a(0:nnzamax-1), c(0:nnzcmax-1)
254
 
      integer rowa(0:nnzamax-1), rowc(0:nnzamax-1)
255
 
      integer ptra(0:n), ptrc(0:jend-jbeg+1)
256
 
      integer ibeg, iend, jbeg, jend
257
 
      integer nnzamax, nnzcmax, irow, jcol, ierr
258
 
 
259
 
      integer nnza, nc, j, iabeg, ia, k, cumsum
260
 
 
261
 
      nc = jend-jbeg+1
262
 
 
263
 
      nnza = ptra(n)
264
 
      nnzc = ierr
265
 
 
266
 
      if (jcol .lt. jbeg) jcol = jbeg
267
 
 
268
 
      do j = jcol, jend
269
 
         jc = j - jbeg
270
 
c   look through row indices in this column, copy all those that are
271
 
c      in required range.
272
 
         iabeg = max(irow,ptra(j))
273
 
         do ia = iabeg, ptra(j+1)-1
274
 
            ra = rowa(ia)
275
 
            if ((ra.le.iend).and.(ra.ge.ibeg)) then
276
 
               if (nnzc.ge.nnzcmax) goto 999
277
 
               c(nnzc) = a(ia)
278
 
               rowc(nnzc) = ra-ibeg
279
 
               ptrc(jc+1) = ptrc(jc+1) + 1
280
 
               nnzc = nnzc + 1
281
 
            end if
282
 
         end do
283
 
      end do
284
 
 
285
 
c   successful completion (fix the ptr array)
286
 
      cumsum = 0
287
 
      do k = 1, n
288
 
         cumsum = cumsum + ptrc(k)
289
 
         ptrc(k) = cumsum         
290
 
      end do
291
 
      return
292
 
 
293
 
 
294
 
      return
295
 
 
296
 
 999  continue
297
 
      ierr = nnzc
298
 
      irow = ia
299
 
      jcol = j
300
 
      return
301
 
      end subroutine <_c>cscextract
302
 
 
303
 
 
304
 
      
305
 
      subroutine <_c>diatocsc(m,n,diags,numdia,diasize,offsets,
306
 
     $     a,rowa,ptra,nzmax, ierr)
307
 
      
308
 
      integer n, numdia, diasize, nzmax, ierr
309
 
      <_t> diags(0:numdia-1, 0:diasize-1), a(0:nzmax-1)
310
 
      <_t> val
311
 
      integer offsets(0:numdia-1), rowa(0:nzmax-1)
312
 
      integer ptra(0:n)
313
 
 
314
 
      integer col, nnza, jj, row, idiag, ia
315
 
 
316
 
      nnza = 0
317
 
      ierr = 0
318
 
c      write (*,*) "SDSD", m, n, numdia, diasize, nzmax
319
 
      do col = 0, n-1
320
 
c   Loop through the offsets 
321
 
         do jj = 0, numdia-1
322
 
            row = col - offsets(jj)
323
 
            if ((row.ge.0).and.(row.lt.m)) then
324
 
               idiag = min(row, col)
325
 
               val = diags(jj, idiag)
326
 
               if (val.ne.0.0) then
327
 
                  if (nnza.ge.nzmax) goto 999
328
 
c   find place to insert this row (inserted in sorted order)
329
 
                  ia = ptra(col)
330
 
 10               if ((ia.lt.ptra(col+1)).and.(rowa(ia).lt.row)) then
331
 
                     ia = ia + 1
332
 
                     goto 10
333
 
                  end if
334
 
                  do iia = nnza, ia+1, -1
335
 
                     a(iia) = a(iia-1)
336
 
                     rowa(iia) = rowa(iia-1)
337
 
                  end do
338
 
                  a(ia) = val
339
 
                  rowa(ia) = row
340
 
                  nnza = nnza + 1
341
 
               end if
342
 
            end if
343
 
         end do
344
 
         ptra(col+1) = nnza
345
 
      end do
346
 
 
347
 
      return 
348
 
 
349
 
 999  continue
350
 
      ierr = 1
351
 
      return
352
 
      
353
 
 
354
 
      end subroutine <_c>diatocsc