← Back to branch summary

~ubuntu-branches/ubuntu/karmic/python-scipy/karmic

« back to all changes in this revision

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

  • Committer: Bazaar Package Importer
  • Author(s): Daniel T. Chen (new)
  • Date: 2005-03-16 02:15:29 UTC
  • Revision ID: james.westby@ubuntu.com-20050316021529-xrjlowsejs0cijig
Tags: upstream-0.3.2
Import upstream version 0.3.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/sparse/sparsetools/spblas.f.src
 
1
c    -*- fortran -*-
 
2
c    author:  travis e. oliphant, 2004
 
3
c
 
4
c    computational tools for sparse matrices in <_t>
 
5
c
 
6
c    add two matrices in compressed sparse column format
 
7
c        this assumes that the rowa indices are sorted for each column
 
8
c        and indexes are 0-based
 
9
 
 
10
      subroutine <_c>cscadd(n,a,rowa,ptra,nnzamax,b,rowb,ptrb,nnzbmax,
 
11
     $     c,rowc,ptrc,nnzcmax,ierr)
 
12
      <_t> a(0:nnzamax-1), b(0:nnzbmax-1), c(0:nnzcmax-1), val
 
13
      integer n, rowa(0:nnzamax-1), rowb(0:nnzbmax-1), rowc(0:nnzcmax-1)
 
14
      integer ptra(0:n), ptrb(0:n), ptrc(0:n), nnzamax, nnzbmax
 
15
      integer ia, ib, ra, rb, j, nnzcmax
 
16
      integer ierr, nnzc
 
17
      
 
18
      ierr=0
 
19
      nnzc=0
 
20
      ia=ptra(0)
 
21
      ib=ptrb(0)
 
22
 
 
23
c  is there a need to initialize the output arrays?  assume no for now.
 
24
c      do j = 0, n
 
25
c         ptrc(j) = 0
 
26
c      end do
 
27
c      do j = 0, nnzcmax-1
 
28
c         rowc(j) = 0
 
29
c         c(j) = 0.0
 
30
c      end do
 
31
 
 
32
      do j = 0, n-1
 
33
         iamax = ptra(j+1)
 
34
         ibmax = ptrb(j+1)
 
35
         do while ((ia.lt.iamax).and.(ib.lt.ibmax))
 
36
            ra = rowa(ia)
 
37
            rb = rowb(ib)
 
38
            if (ra.eq.rb) then
 
39
c     this is a common element
 
40
               val = a(ia) + b(ib)
 
41
               ia = ia + 1
 
42
               ib = ib + 1
 
43
               if (val.eq.0.0) goto 10
 
44
               if (nnzc.ge.nnzcmax) goto 999
 
45
               c(nnzc) = val
 
46
               rowc(nnzc) = ra
 
47
            else if (ra .lt. rb) then
 
48
c     a has this but not b
 
49
               val = a(ia)
 
50
               ia = ia + 1 
 
51
               if (val.eq.0.0) goto 10 
 
52
               if (nnzc.ge.nnzcmax) goto 999
 
53
               c(nnzc) = val
 
54
               rowc(nnzc) = ra
 
55
            else
 
56
c     b has this but not a
 
57
               val = b(ib)
 
58
               ib = ib + 1
 
59
               if (val.eq.0.0) goto 10
 
60
               if (nnzc.ge.nnzcmax) goto 999
 
61
               c(nnzc) = val
 
62
               rowc(nnzc) = rb
 
63
            end if
 
64
            ptrc(j+1) = ptrc(j+1) + 1
 
65
            nnzc = nnzc + 1
 
66
 10         continue
 
67
         end do
 
68
         
 
69
         if (ia .eq. iamax) then
 
70
c     all finished with a for this column just copy the rest from b
 
71
            do while (ib .lt. ibmax)
 
72
               val = b(ib)
 
73
               rb = rowb(ib)
 
74
               ib = ib+1
 
75
               if (val.ne.0.0) then
 
76
                  if (nnzc.ge.nnzcmax) goto 999
 
77
                  c(nnzc) = val
 
78
                  rowc(nnzc) = rb
 
79
                  ptrc(j+1) = ptrc(j+1) + 1
 
80
                  nnzc = nnzc + 1
 
81
               end if
 
82
            end do
 
83
         else if (ib .eq. ibmax) then
 
84
c     all finished with b for this column just copy the rest from a
 
85
            do while (ia .lt. iamax)
 
86
               val = a(ia)
 
87
               ra = rowa(ia)
 
88
               ia = ia + 1
 
89
               if (val.ne.0.0) then
 
90
                  if (nnzc.ge.nnzcmax) goto 999
 
91
                  c(nnzc) = val
 
92
                  rowc(nnzc) = ra
 
93
                  ptrc(j+1) = ptrc(j+1) + 1
 
94
                  nnzc = nnzc + 1
 
95
               end if
 
96
            end do
 
97
         end if 
 
98
      end do
 
99
 
 
100
c   successful completion (fix the ptr array)
 
101
      cumsum = 0
 
102
      do k = 1, n
 
103
         cumsum = cumsum + ptrc(k)
 
104
         ptrc(k) = cumsum         
 
105
      end do
 
106
      return
 
107
      
 
108
      
 
109
 999  continue
 
110
      ierr = 1
 
111
      return
 
112
      end subroutine <_c>cscadd
 
113
 
 
114
c     element-by-element multiplication
 
115
c     can have at most the minimum of nnzamax and nnzbmax
 
116
 
 
117
      subroutine <_c>cscmul(n,a,rowa,ptra,nnzamax,b,rowb,ptrb,nnzbmax,
 
118
     $     c,rowc,ptrc,nnzcmax,ierr)
 
119
      <_t> a(0:nnzamax-1), b(0:nnzbmax-1), c(0:nnzcmax-1), val
 
120
      integer n, rowa(0:nnzamax-1), rowb(0:nnzbmax-1), rowc(0:nnzcmax-1)
 
121
      integer ptra(0:n), ptrb(0:n), ptrc(0:n), nnzamax, nnzbmax
 
122
      integer ia, ib, ra, rb, j, nnzcmax
 
123
      integer ierr, k, cumsum
 
124
      
 
125
      ierr=0
 
126
      nnzc=0
 
127
      ia=ptra(0)
 
128
      ib=ptrb(0)
 
129
      do j = 0, n-1
 
130
         iamax = ptra(j+1)
 
131
         ibmax = ptrb(j+1)
 
132
         do while ((ia.lt.iamax).and.(ib.lt.ibmax))
 
133
            ra = rowa(ia)
 
134
            rb = rowb(ib)
 
135
            if (ra.eq.rb) then
 
136
c     this is a common element
 
137
               val = a(ia)*b(ib)
 
138
               ia = ia + 1
 
139
               ib = ib + 1
 
140
               if (val.eq.0.0) goto 10
 
141
               if (nnzc.ge.nnzcmax) goto 999
 
142
               c(nnzc) = val
 
143
               rowc(nnzc) = ra
 
144
               ptrc(j+1) = ptrc(j+1) + 1
 
145
               nnzc = nnzc + 1
 
146
            else if (ra .lt. rb) then
 
147
c     a has this but not b
 
148
               ia = ia + 1
 
149
            else
 
150
c     b has this but not a
 
151
               ib = ib + 1
 
152
            end if
 
153
 10         continue
 
154
         end do
 
155
      end do
 
156
      
 
157
c   successful completion (fix the ptr array)
 
158
      cumsum = 0
 
159
      do k = 1, n
 
160
         cumsum = cumsum + ptrc(k)
 
161
         ptrc(k) = cumsum         
 
162
      end do
 
163
 
 
164
      return
 
165
      
 
166
 999  continue
 
167
      ierr = 1
 
168
      return
 
169
      end subroutine <_c>cscmul
 
170
 
 
171
 
 
172
c  matrix-vector multiplication
 
173
 
 
174
 
 
175
      subroutine <_c>cscmux(a,rowa,ptra,nnzamax,ncol,x,mrow,y)
 
176
      <_t> a(0:nnzamax-1),x(0:ncol-1),y(0:mrow-1)
 
177
      integer rowa(0:nnzamax-1), ptra(0:ncol)
 
178
      integer nnzamax, mrow, ncol
 
179
      integer i, j, ia, ra
 
180
      
 
181
      do i = 0, mrow-1
 
182
         y(i) = 0.0
 
183
      end do
 
184
 
 
185
      do j = 0, ncol-1
 
186
         do ia = ptra(j), ptra(j+1)-1
 
187
            ra = rowa(ia)
 
188
            y(ra) = y(ra) + a(ia)*x(j)
 
189
         end do
 
190
      end do
 
191
 
 
192
      return
 
193
      end subroutine <_c>cscmux
 
194
 
 
195
 
 
196
      subroutine <_c>csrmux(a,cola,ptra,nnzamax,ncol,x,mrow,y)
 
197
      <_t> a(0:nnzamax-1),x(0:ncol-1),y(0:mrow-1)
 
198
      integer cola(0:nnzamax-1), ptra(0:mrow)
 
199
      integer nnzamax, mrow, ncol
 
200
      integer i, ja, ca
 
201
      
 
202
      do i = 0, mrow-1
 
203
         y(i) = 0.0
 
204
         do ja = ptra(i), ptra(i+1)-1
 
205
            ca = cola(ja)
 
206
            y(i) = y(i) + a(ja)*x(ca)           
 
207
         end do
 
208
      end do
 
209
 
 
210
      return
 
211
      end subroutine <_c>csrmux
 
212
 
 
213
 
 
214
c     matrix multiplication
 
215
c             c = a * b
 
216
c
 
217
c     where c is mxn csc matrix 
 
218
c           a is mxk csc matrix
 
219
c           b is kxn csr matrix
 
220
c
 
221
c     irow and kcol give position to start 
 
222
c     nnzc gives current last element in output array (initialize to 0)
 
223
c     intended to recall so that calculation can continue 
 
224
c     if memory ran-out at last attempt
 
225
 
 
226
 
 
227
 
 
228
      subroutine <_c>cscmucsr(m,k,n, a,rowa,ptra,nnzamax, b,colb,ptrb,
 
229
     $     nnzbmax, c,rowc,ptrc,nnzcmax, irow,kcol,ierr)
 
230
 
 
231
      <_t> a(0:nnzamax-1), b(0:nnzbmax-1), c(0:nnzcmax-1)
 
232
      integer rowa(0:nnzamax-1), colb(0:nnzbmax-1), rowc(0:nnzcmax-1)
 
233
      integer ptra(0:k), ptrb(0:k), ptrc(0:n)
 
234
      integer nnzamax, nnzbmax, nnzcmax
 
235
      integer m, k, n, irow, kcol, ierr
 
236
 
 
237
      integer kk, ii, jjb, jb, cb, ia, nnzc
 
238
      <_t> val
 
239
 
 
240
      ierr = 0
 
241
      nnzc = ptrc(n)
 
242
      do kk = kcol, n-1
 
243
         do ii = irow, m-1
 
244
            if (nnzc.ge.nnzcmax) goto 999
 
245
            val = 0.0
 
246
c   loop through the column array of b using the ptrb array
 
247
c        so that we know which row we are on.
 
248
            do jjb = 0, k-1
 
249
               do jb = ptrb(jjb), ptrb(jjb+1)-1
 
250
                  cb = colb(jb)
 
251
                  if (cb.eq.kk) then
 
252
c     see if aij is nonzero and if so add to val
 
253
                     do ia = ptra(jjb), ptra(jjb+1)-1
 
254
                        if (rowa(ia).eq.ii) then
 
255
                           val = val + a(ia)*b(jb)
 
256
                        end if
 
257
                     end do
 
258
                  end if
 
259
               end do
 
260
            end do
 
261
            if (val.ne.0.0) then
 
262
c     there is a non-zero value for this value of i and k
 
263
               c(nnzc) = val
 
264
               rowc(nnzc) = ii
 
265
               ptrc(kk+1) = ptrc(kk+1) + 1
 
266
               nnzc = nnzc+1               
 
267
            end if
 
268
         end do         
 
269
      end do
 
270
 
 
271
      return
 
272
 
 
273
 999  continue
 
274
      kcol = kk
 
275
      irow = ii
 
276
      ierr = 1
 
277
 
 
278
      return
 
279
      end subroutine <_c>cscmucsr
 
280
 
 
281
 
 
282
 
 
283
c     matrix multiplication
 
284
c             c = a * b
 
285
c
 
286
c     where c is mxn csc matrix 
 
287
c           a is mxk csr matrix
 
288
c           b is kxn csc matrix
 
289
c
 
290
c     irow and jcol give position to start 
 
291
c     bnum gives current size of output array
 
292
c     intended to recall so that calculation can continue 
 
293
c     where it left off if memory runs-out during calculation
 
294
c
 
295
      subroutine <_c>csrmucsc(m,n,a,cola,ptra,nnzamax,b,rowb,ptrb,
 
296
     $     nnzbmax,c,rowc,ptrc,nnzcmax,irow,kcol,ierr)
 
297
 
 
298
      <_t> a(0:nnzamax-1), b(0:nnzbmax-1), c(0:nnzcmax-1)
 
299
      integer cola(0:nnzamax-1), rowb(0:nnzbmax-1), rowc(0:nnzcmax-1)
 
300
      integer ptra(0:m), ptrb(0:n), ptrc(0:n)
 
301
      integer nnzamax, nnzbmax, nnzcmax
 
302
      integer m, n, irow, kcol, ierr
 
303
 
 
304
      integer kk, ii, jjb, jb, cb, ia, nnzc
 
305
      <_t> val
 
306
 
 
307
      ierr = 0
 
308
      nnzc = ptrc(n)
 
309
      do kk = kcol, n-1
 
310
         do ii = irow, m-1
 
311
            if (nnzc.ge.nnzcmax) goto 999
 
312
            val = 0.0
 
313
c   loop through the non-zero rows of b
 
314
            do ib = ptrb(kk), ptrb(kk+1)-1
 
315
               rb = rowb(ib)
 
316
               do ja = ptra(ii), ptra(ii+1)-1
 
317
                  if (cola(ja).eq.rb) then
 
318
                     val = val + a(ja)*b(ib)
 
319
                  end if
 
320
               end do
 
321
            end do
 
322
            if (val.ne.0.0) then
 
323
c     there is a non-zero value for this value of i and k
 
324
               c(nnzc) = val
 
325
               rowc(nnzc) = ii
 
326
               ptrc(kk+1) = ptrc(kk+1) + 1
 
327
               nnzc = nnzc+1 
 
328
            end if
 
329
         end do
 
330
      end do
 
331
 
 
332
      return
 
333
 
 
334
 999  continue
 
335
      kcol = kk
 
336
      irow = ii
 
337
      ierr = 1
 
338
 
 
339
      return
 
340
      end subroutine <_c>csrmucsc
 
341
 
 
342
 
 
343
c     matrix-matrix multiplication
 
344
c
 
345
c          c = a * b    where  a, b, and c are 
 
346
c                       compressed sparse column matrices
 
347
c
 
348
c     or it computes
 
349
c
 
350
c          c = b * a    where  b, a, and c are 
 
351
c                       compressed sparse row matrices
 
352
c   
 
353
c     intentended for re-entry if nnzcmax is not big enough
 
354
 
 
355
      subroutine <_c>cscmucsc(m,k,n,a,rowa,ptra,nnzamax,b,rowb,ptrb,
 
356
     $     nnzbmax,c,rowc,ptrc,nnzcmax,irow,kcol,ierr)
 
357
 
 
358
      <_t> a(0:nnzamax-1), b(0:nnzbmax-1), c(0:nnzcmax-1)
 
359
      integer rowa(0:nnzamax-1), rowb(0:nnzbmax-1), rowc(0:nnzcmax-1)
 
360
      integer ptra(0:k), ptrb(0:n), ptrc(0:n)
 
361
      integer nnzamax, nnzbmax, nnzcmax
 
362
      integer m, k, n, irow, kcol, ierr
 
363
 
 
364
      integer kk, ii, jjb, jb, cb, ia, nnzc, cumsum
 
365
      <_t> val
 
366
      
 
367
      ierr = 0
 
368
      nnzc = ierr
 
369
      do kk = kcol, n-1
 
370
         do ii = irow, m-1
 
371
            val = 0.0
 
372
c   loop through the row array of b for this column
 
373
            do jb = ptrb(kk), ptrb(kk+1)-1
 
374
               ra = rowb(jb)
 
375
               do ia = ptra(ra), ptra(ra+1)-1
 
376
c     see if aij is nonzero and if so add to val
 
377
                  if (rowa(ia).eq.ii) then
 
378
                        val = val + a(ia)*b(jb)
 
379
                  end if
 
380
               end do
 
381
            end do
 
382
            if (val.ne.0.0) then
 
383
c     there is a non-zero value for this value of i and k
 
384
               if (nnzc.ge.nnzcmax) goto 999
 
385
               c(nnzc) = val
 
386
               rowc(nnzc) = ii
 
387
               ptrc(kk+1) = ptrc(kk+1) + 1
 
388
               nnzc = nnzc+1               
 
389
            end if
 
390
         end do         
 
391
      end do
 
392
 
 
393
c   successful completion (fix the ptr array)
 
394
      cumsum = 0
 
395
      do kk = 1, n
 
396
         cumsum = cumsum + ptrc(kk)
 
397
         ptrc(kk) = cumsum         
 
398
      end do
 
399
      return
 
400
      
 
401
 999  continue
 
402
      kcol = kk
 
403
      irow = ii
 
404
      ierr = nnzc
 
405
      
 
406
      return
 
407
      end subroutine <_c>cscmucsc
 
408
     
 
409