← Back to branch summary

~ubuntu-branches/ubuntu/karmic/scilab/karmic

« back to all changes in this revision

Viewing changes to routines/sparse/dsposp.f

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-21 16:57:43 UTC
  • Revision ID: james.westby@ubuntu.com-20020321165743-e9mv12c1tb1plztg
Tags: upstream-2.6
Import upstream version 2.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
routines/sparse/dsposp.f
 
1
      subroutine dsposp(op,ma,na,a,nela,inda,mb,nb,b,nelb,indb,
 
2
c     Copyright INRIA
 
3
     $     nelc,indc,ierr)
 
4
c!purpose
 
5
c     compare the elements of two  sparse matrices.
 
6
c!parameters 
 
7
c     op   : integer flag which specifies the comparison to perform
 
8
c            see routine dcompa for more precision
 
9
c     ma,na: row and column dimension of the a matrix  
 
10
c     mb,nb: row and column dimension of the b matrix  
 
11
c     a,b  : arrays. 
 
12
c     Contain non zero elements of first,second sparse matrices.
 
13
c     nela :integer: number of non zero elements of a
 
14
c     nelb :integer: number of non zero elements of b
 
15
c     nelc :integer: 
 
16
c                   on entry maximum number  of non zero elements of c
 
17
c                   on return number of non zero elements of c
 
18
c     inda : a matrix control data:
 
19
c            inda(i) 1<=i<=nr contains the number of ith row non zero elements
 
20
c             of a
 
21
c            inda(m+i) 1<=i<=nela column index of each non zero element
 
22
c     indb : b matrix control data:
 
23
c            indb(i) 1<=i<=nr contains the number of ith row non zero elements
 
24
c            of b
 
25
c             indb(m+i) 1<=i<=nelb column index of each non zero element
 
26
c     
 
27
c     indc : on return contains c matrix control data:
 
28
c            indc(i) 1<=i<=nr contains the number of ith row non zero elements
 
29
c            of c
 
30
c            indc(m+i) 1<=i<=nelb column index of each non zero element
 
31
c     ierr : if non zero initial value of nelc is to small
 
32
c     !
 
33
      double precision a(*),b(*),t
 
34
      integer op,nela,inda(*),nelb,indb(*),nelc,indc(*),ierr
 
35
c     
 
36
      integer jc,ka,kb,jb,i,ja,j1,j2
 
37
      logical dcompa,z
 
38
      external dcompa
 
39
c     
 
40
      nr=max(ma,mb)
 
41
      nc=max(na,nb)
 
42
c     
 
43
      nelmx=nelc
 
44
      ierr=0
 
45
 
 
46
c     jc counts elements of c.
 
47
      jc     = 1
 
48
c     ka,kb are numbers in first i rows of a,b.
 
49
      ka     = 1
 
50
      kb     = 1
 
51
      kc     = 1
 
52
c     jb counts elements of b.
 
53
      jb     = 1
 
54
c     i counts rows of a,b,c.
 
55
      if(ma*na.eq.1.and.mb*nb.gt.1) then
 
56
c     compare all element of b with scalar a
 
57
         t=0.0d0
 
58
         if(inda(1).eq.1) t=a(1)
 
59
         z=dcompa(t,0.0d0,op)   
 
60
         do 10 i=1,nr
 
61
            indc(i)=0
 
62
            nirb=indb(i)
 
63
            jb=kb        
 
64
            jc=kc
 
65
            if(nirb.eq.0) then
 
66
               do 03 j=1,nc
 
67
                  if (z) then
 
68
                     if(jc+1.gt.nelmx) goto 99
 
69
                     indc(nr+jc)=j
 
70
                     jc=jc+1
 
71
                  endif
 
72
 03            continue
 
73
            else
 
74
               j2=indb(nr+jb)
 
75
               do 04 j=1,nc
 
76
                  if(j2.eq.j) then
 
77
                     if (dcompa(t,b(jb),op)) then
 
78
                        if(jc+1.gt.nelmx) goto 99
 
79
                        indc(nr+jc)=j
 
80
                        jc=jc+1
 
81
                     endif 
 
82
                     if(jb-kb+1.lt.nirb) jb=jb+1
 
83
                     j2=indb(nr+jb)
 
84
                  else
 
85
                     if (z) then
 
86
                        if(jc+1.gt.nelmx) goto 99
 
87
                        indc(nr+jc)=j
 
88
                        jc=jc+1
 
89
                     endif 
 
90
                  endif
 
91
 04            continue
 
92
            endif
 
93
            kb=kb+indb(i)
 
94
            indc(i)=jc-kc
 
95
            kc=jc
 
96
 10      continue
 
97
 
 
98
      elseif(ma*na.gt.1.and.mb*nb.eq.1) then
 
99
c     compare all elements of a with scalar b  
 
100
         t=0.0d0
 
101
         if(indb(1).eq.1) t=b(1)
 
102
         z=dcompa(0.0d0,t,op)
 
103
         do 20 i=1,nr
 
104
            indc(i)=0
 
105
            nira=inda(i)
 
106
            ja=ka       
 
107
            jc=kc
 
108
            if(nira.eq.0) then
 
109
               if(z) then
 
110
                  if(kc+nc.gt.nelmx) goto 99
 
111
                  indc(i)=nc
 
112
                  do 11 j=1,nc
 
113
                     indc(nr+kc-1+j)=j
 
114
 11               continue
 
115
                  jc=kc+nc
 
116
               endif
 
117
            else
 
118
               j1=inda(nr+ja)
 
119
               do 12 j=1,nc
 
120
                  if(j1.eq.j) then
 
121
                     if (dcompa(a(ja),t,op)) then
 
122
                        if(jc+1.gt.nelmx) goto 99
 
123
                        indc(nr+jc)=j
 
124
                        jc=jc+1
 
125
                     endif
 
126
                     if(ja-ka+1.lt.nira) ja=ja+1
 
127
                     j1=inda(nr+ja)
 
128
                  elseif(z) then
 
129
                     if(jc+1.gt.nelmx) goto 99
 
130
                     indc(nr+jc)=j
 
131
                     jc=jc+1
 
132
                  endif
 
133
 12            continue
 
134
            endif
 
135
            indc(i)=jc-kc
 
136
            ka=ka+nira
 
137
            kc=jc
 
138
 20      continue 
 
139
      else
 
140
         z=dcompa(0.0d0,0.0d0,op)   
 
141
         do 30 i=1,nr
 
142
            indc(i)=0
 
143
            nira=inda(i)
 
144
            nirb=indb(i)
 
145
            ja=ka
 
146
            jb=kb        
 
147
            jc=kc
 
148
            if(nira.eq.0) then
 
149
               if(nirb.eq.0) then
 
150
                  if(z) then
 
151
                     if(kc+nc.gt.nelmx) goto 99
 
152
                     indc(i)=nc
 
153
                     do 21 j=1,nc
 
154
                        indc(nr+kc-1+j)=j
 
155
 21                  continue
 
156
                     jc=kc+nc
 
157
                  endif
 
158
               else
 
159
                  j2=indb(nr+jb)
 
160
                  do 22 j=1,nc
 
161
                     if(j2.eq.j) then
 
162
                        if (dcompa(0.0d0,b(jb),op)) then
 
163
                           if(jc+1.gt.nelmx) goto 99
 
164
                           indc(nr+jc)=j
 
165
                           jc=jc+1
 
166
                        endif
 
167
                        if(jb-kb+1.lt.nirb) jb=jb+1
 
168
                        j2=indb(nr+jb)
 
169
                     elseif(z) then
 
170
                        if(jc+1.gt.nelmx) goto 99
 
171
                        indc(nr+jc)=j
 
172
                        jc=jc+1
 
173
                     endif
 
174
 22               continue
 
175
               endif
 
176
            else
 
177
               if(nirb.eq.0) then
 
178
                  j1=inda(nr+ja)
 
179
                  do 23 j=1,nc
 
180
                     if(j1.eq.j) then
 
181
                        if (dcompa(a(ja),0.0d0,op)) then
 
182
                           if(jc+1.gt.nelmx) goto 99
 
183
                           indc(nr+jc)=j
 
184
                           jc=jc+1
 
185
                        endif
 
186
                        if(ja-ka+1.lt.nira) ja=ja+1
 
187
                        j1=inda(nr+ja)
 
188
                     elseif(z) then
 
189
                        if(jc+1.gt.nelmx) goto 99
 
190
                        indc(nr+jc)=j
 
191
                        jc=jc+1
 
192
                     endif
 
193
 23               continue
 
194
               else
 
195
                  j1=inda(nr+ja)
 
196
                  j2=indb(nr+jb)
 
197
                  do 24 j=1,nc
 
198
                     if(j1.eq.j) then
 
199
                        if(j2.eq.j) then
 
200
                           if (dcompa(a(ja),b(jb),op)) then
 
201
                              if(jc+1.gt.nelmx) goto 99
 
202
                              indc(nr+jc)=j
 
203
                              jc=jc+1
 
204
                           endif 
 
205
                           if(ja-ka+1.lt.nira) ja=ja+1
 
206
                           if(jb-kb+1.lt.nirb) jb=jb+1
 
207
                           j1=inda(nr+ja)
 
208
                           j2=indb(nr+jb)
 
209
                        else
 
210
                           if (dcompa(a(ja),0.0d0,op)) then
 
211
                              if(jc+1.gt.nelmx) goto 99
 
212
                              indc(nr+jc)=j
 
213
                              jc=jc+1
 
214
                           endif 
 
215
                           if(ja-ka+1.lt.nira) ja=ja+1
 
216
                           j1=inda(nr+ja)
 
217
                        endif
 
218
                     else
 
219
                        if(j2.eq.j) then
 
220
                           if (dcompa(0.0d0,b(jb),op)) then
 
221
                              if(jc+1.gt.nelmx) goto 99
 
222
                              indc(nr+jc)=j
 
223
                              jc=jc+1
 
224
                           endif
 
225
                           if(jb-kb+1.lt.nirb) jb=jb+1
 
226
                           j2=indb(nr+jb)
 
227
                        elseif(z) then
 
228
                           if(jc+1.gt.nelmx) goto 99
 
229
                           indc(nr+jc)=j
 
230
                           jc=jc+1
 
231
                        endif 
 
232
                     endif
 
233
 24               continue
 
234
               endif
 
235
            endif
 
236
            ka=ka+inda(i)
 
237
            kb=kb+indb(i)
 
238
            indc(i)=jc-kc
 
239
            kc=jc
 
240
 30      continue
 
241
      endif
 
242
      nelc  = jc-1
 
243
      return
 
244
c     error messages.
 
245
 99   ierr=1
 
246
c     no more place for c
 
247
 
 
248
      return 
 
249
      end
 
250
 
 
251