← Back to branch summary

~ubuntu-branches/ubuntu/hoary/scilab/hoary

« back to all changes in this revision

Viewing changes to routines/interf/matops.f

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2005-01-09 22:58:21 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20050109225821-473xr8vhgugxxx5j
Tags: 3.0-12
http://bugs.debian.org/255869
changed configure.in to build scilab's own malloc.o, closes: #255869

Show diffs side-by-side

added added

removed removed

Lines of Context:
routines/interf/matops.f
26
26
c     
27
27
      fun = 0
28
28
c     
29
 
c        cconc  extrac insert rconc
 
29
c        cconc   insert extrac rconc
30
30
      goto(75  ,  95  ,  78   ,76) op
31
31
c     
32
32
c           :  +  -  * /  \  =          '
238
238
         go to 41
239
239
      elseif (mn2.eq.1) then
240
240
c     .   vector+const
 
241
         err=l1+mn1*(itr+1)-lstk(bot)
 
242
         if(err.gt.0) then
 
243
            call error(17)
 
244
            return
 
245
         endif
241
246
         call dadd(mn1,stk(l2),0,stk(l1),1)
242
247
         if(it2+2*it1.eq.1) call unsfdcopy(mn1,stk(l2+mn2),0,
243
248
     $        stk(l1+mn1),1)
246
251
         istk(il1+3)=itr
247
252
      elseif (mn1.eq.1) then
248
253
c     .  cst+vector
 
254
         err=l1+mn2*(itr+1)-lstk(bot)
 
255
         if(err.gt.0) then
 
256
            call error(17)
 
257
            return
 
258
         endif
249
259
         cstr=stk(l1)
250
260
         csti=stk(l1+1)
251
261
         call unsfdcopy((it2+1)*mn2,stk(l2),1,stk(l1),1)
262
272
            call error(8)
263
273
            return
264
274
         endif
 
275
         err=l1+mn1*(itr+1)-lstk(bot)
 
276
         if(err.gt.0) then
 
277
            call error(17)
 
278
            return
 
279
         endif
265
280
         call dadd(mn1,stk(l2),1,stk(l1),1)
266
281
         if(it2+2*it1.eq.1) then
267
282
            call unsfdcopy(mn1,stk(l2+mn1),1,stk(l1+mn1),1)
604
619
            return
605
620
         endif
606
621
         lr=l2+mn2*(it2+1)
607
 
         err=lr+m1*n2*(itr+1)-lstk(bot)
608
 
         if(err.gt.0) then
 
622
c     .  m1*n2 may overflow
 
623
         temp=float(lr)+float(m1)*n2*(itr+1)-lstk(bot)
 
624
         if(temp.gt.0.0d0) then
 
625
            err=int(temp)
609
626
            call error(17)
610
627
            return
611
628
         endif
612
629
         if(it1*it2.ne.1) then
613
 
            call dmmul(stk(l1),m1,stk(l2),m2,stk(lr),m1,m1,n1,n2)
614
 
            if(it1.eq.1) call dmmul(stk(l1+mn1),m1,stk(l2),m2
615
 
     $           ,stk(lr+m1*n2),m1,m1,n1,n2)
616
 
            if(it2.eq.1) call dmmul(stk(l1),m1,stk(l2+mn2),m2
617
 
     $           ,stk(lr+m1*n2),m1,m1,n1,n2)
 
630
*           remplacement de dmmul par dgemm (Bruno le 31/10/2001)
 
631
            call dgemm('n','n',m1,n2,n1,1.d0,stk(l1),m1,stk(l2),m2,
 
632
     $           0.d0,stk(lr),m1)
 
633
            if(it1.eq.1) call dgemm('n','n',m1,n2,n1,1.d0,stk(l1+mn1),
 
634
     $           m1,stk(l2),m2,0.d0,stk(lr+m1*n2),m1)
 
635
            if(it2.eq.1) call dgemm('n','n',m1,n2,n1,1.d0,stk(l1),m1,
 
636
     $           stk(l2+mn2),m2,0.d0,stk(lr+m1*n2),m1)      
618
637
         else
619
638
c     .     a et a2 sont complexes
620
639
            call wmmul(stk(l1),stk(l1+mn1),m1,stk(l2),stk(l2
675
694
            call error(14)
676
695
            return
677
696
         endif
678
 
         if (m2 .eq. n2) fun = 1
679
 
         if (m2 .ne. n2) fun = 4
680
 
         fin = -1
681
697
         top = top+1
682
 
         rhs = 2
 
698
         rhs = 2        
 
699
         call intslash('slash')
 
700
         if (fin.ge.0) call putlhsvar()
 
701
c         if (m2 .eq. n2) fun = 1
 
702
c         if (m2 .ne. n2) fun = 4
 
703
c         fin = -1
683
704
      else
684
705
c     .  vector / cst
685
706
         if(m2.lt.0.and.mn1.ne.1) then 
773
794
            call error(14)
774
795
            return
775
796
         endif
776
 
         if (m1 .eq. n1) fun = 1
777
 
         if (m1 .ne. n1) fun = 4
778
797
         top = top+1
779
 
         fin = -2
780
 
         rhs = 2
 
798
         rhs = 2 
 
799
         call intbackslash('backslash')
 
800
         if (fin.ge.0) call putlhsvar()
 
801
c         if (m1 .eq. n1) fun = 1
 
802
c         if (m1 .ne. n1) fun = 4
 
803
c         fin = -2
781
804
      else
782
805
c     .  cst \ vector
783
806
         if(m1.lt.0.and.mn2.ne.1) then 
900
923
      if(it2.eq.0) then
901
924
         if(it1.eq.0) then
902
925
            call ddpow1(mn,stk(l1),inc1,stk(l2),inc2,
903
 
     $           stk(lw),stk(lw+mn),1,err,itr)
 
926
     $           stk(lw),stk(lw+mn),1,ierr,itr)
904
927
         else
905
928
            call wdpow1(mn,stk(l1),stk(l1+mn1),inc1,stk(l2),inc2,
906
 
     $           stk(lw),stk(lw+mn),1,err)
 
929
     $           stk(lw),stk(lw+mn),1,ierr)
907
930
         endif
908
931
      else
909
932
         if(it1.eq.0) then
910
933
            call dwpow1(mn,stk(l1),inc1,stk(l2),stk(l2+mn2),inc2,
911
 
     &           stk(lw),stk(lw+mn),1,err)
 
934
     &           stk(lw),stk(lw+mn),1,ierr)
912
935
         else
913
936
            call wwpow1(mn,stk(l1),stk(l1+mn1),inc1,stk(l2),stk(l2+mn2),
914
 
     &           inc2,stk(lw),stk(lw+mn),1,err)
 
937
     &           inc2,stk(lw),stk(lw+mn),1,ierr)
915
938
         endif
916
939
      endif
917
 
      if(err.eq.1) then
 
940
      if(ierr.eq.1) then
918
941
         call error(30)
919
942
         return
920
943
      endif
921
 
      if(err.eq.2) then
 
944
      if(ierr.eq.2) then
922
945
         if(ieee.eq.0) then
923
946
            call error(27)
924
947
            return
1100
1123
      endif
1101
1124
c     
1102
1125
      if (nexp.le.0) then
1103
 
         fun=10
1104
 
         fin=1
1105
1126
         rhs=1
1106
 
         call matlu
1107
 
         if(err.gt.0) return
 
1127
         call intinv('pow')
 
1128
         call putlhsvar()
 
1129
c         call matlu
 
1130
         if(err.gt.0.or.err1.gt.0) return
1108
1131
         nexp=-nexp
1109
 
         fun=0
1110
1132
      endif
1111
1133
      l2=l1+mn1*(it1+1)
1112
1134
c     
1248
1270
         if(istk(pstk(pt)).eq.16) goto  54
1249
1271
      endif
1250
1272
c
 
1273
      if (st .eq. 0.0d+0) then
 
1274
         istk(il1+1)=1
 
1275
         istk(il1+1)=0
 
1276
         istk(il1+2)=0
 
1277
         istk(il1+3)=0
 
1278
         lstk(top+1)=l1
 
1279
         return
 
1280
      endif
 
1281
 
 
1282
      if(isanan(e1).eq.1.or.isanan(st).eq.1.or.isanan(e2).eq.1) then
 
1283
         stk(l1)=e1+st+e2
 
1284
         istk(il1+1)=1
 
1285
         istk(il1+2)=1
 
1286
         istk(il1+3)=0
 
1287
         lstk(top+1)=l1+1
 
1288
         return
 
1289
      endif
 
1290
 
1251
1291
c     floating point used to avoid integer overflow
1252
1292
      e1r=dble(l1) + max(3.0d0,(e2-e1)/st) - dble(lstk(bot))
 
1293
      if(isanan(e1r).eq.1) then
 
1294
        stk(l1)=eiR
 
1295
         istk(il1+1)=1
 
1296
         istk(il1+2)=1
 
1297
         istk(il1+3)=0
 
1298
         lstk(top+1)=l1+1
 
1299
       return
 
1300
      endif
 
1301
 
1253
1302
      if (e1r .gt. 0.0d0) then
1254
1303
         err=e1r
1255
1304
         call error(17)
1276
1325
      return
1277
1326
c     
1278
1327
c     for clause
1279
 
 54   stk(l1) = e1
1280
 
      stk(l1+1) = st
1281
 
      stk(l1+2) = e2
1282
 
      istk(il1+1)=-3
1283
 
      istk(il1+2)=-1
1284
 
      lstk(top+1)=l1+3
 
1328
 54   if(isanan(e1).eq.1.or.isanan(st).eq.1.or.isanan(e2).eq.1) then
 
1329
         stk(l1) = e1+st+e2
 
1330
         istk(il1+1)=1
 
1331
         istk(il1+2)=1
 
1332
         lstk(top+1)=l1+3
 
1333
      else
 
1334
         stk(l1) = e1
 
1335
         stk(l1+1) = st
 
1336
         stk(l1+2) = e2
 
1337
         istk(il1+1)=-3
 
1338
         istk(il1+2)=-1
 
1339
         lstk(top+1)=l1+3
 
1340
      endif
1285
1341
      return
1286
1342
      end
1287
1343
 
2086
2142
         istk(il1+1)=mn2
2087
2143
         istk(il1+2)=1
2088
2144
         istk(il1+3)=istk(il2+3)
 
2145
         l1=sadr(il1+4)
2089
2146
         call unsfdcopy(mn2*(it2+1),stk(l2),1,stk(l1),1)
2090
 
         lstk(top+1)=sadr(il1+4)+mn2*(it2+1)
 
2147
         lstk(top+1)=l1+mn2*(it2+1)
2091
2148
         return
2092
2149
      endif
2093
2150
c     check and convert indices variable
2274
2331
c     
2275
2332
c     Copyright INRIA
2276
2333
      include '../stack.h'
 
2334
      common /mtlbc/ mmode
2277
2335
c     
2278
2336
      logical isany
2279
2337
      integer top0
2296
2354
c     
2297
2355
      il2=iadr(lstk(top))
2298
2356
      if(istk(il2).lt.0) il2=iadr(istk(il2+1))
 
2357
  
2299
2358
      m2=istk(il2+1)
2300
2359
      n2=istk(il2+2)
2301
2360
      it2=istk(il2+3)
2306
2365
      il1=iadr(lstk(top))
2307
2366
      ilrs=il1
2308
2367
      if(istk(il1).lt.0) il1=iadr(istk(il1+1))
 
2368
      if (istk(il1).eq.10.or.istk(il1).eq.15) then
 
2369
         top=top0
 
2370
         fin=-fin
 
2371
         return
 
2372
      endif
2309
2373
      m1=istk(il1+1)
2310
2374
      n1=istk(il1+2)
2311
2375
      it1=istk(il1+3)
2312
2376
      l1=sadr(il1+4)
2313
2377
      mn1=m1*n1
2314
 
 
2315
 
 
2316
2378
c     arg3(arg1)=arg2
2317
2379
c     
2318
2380
      if (istk(il2)*istk(il1).eq.0) then
2497
2559
         l1=sadr(ilrs+4)
2498
2560
         call unsfdcopy(mnr*(itr+1),stk(lr),1,stk(l1),1)
2499
2561
         istk(ilrs)=1
 
2562
         if(mmode.eq.1.and.nr.eq.1.and.m3.eq.0) then
 
2563
         istk(ilrs+1)=nr
 
2564
         istk(ilrs+2)=mr
 
2565
         else
2500
2566
         istk(ilrs+1)=mr
2501
2567
         istk(ilrs+2)=nr
 
2568
         endif
2502
2569
         istk(ilrs+3)=itr
2503
2570
         lstk(top+1)=l1+mnr*(itr+1)
2504
2571
      else
3072
3139
            endif
3073
3140
            return
3074
3141
         endif
 
3142
         if (.true.) then 
 
3143
c           add explicit nan tests when requested 
 
3144
            call idcmp(stk(l1),stk(l2),mn1,istk(il1+3),op)
 
3145
         else
3075
3146
         do 132 i=0,mn1-1
3076
3147
            e1=stk(l1+i)
3077
3148
            e2=stk(l2+i)
3093
3164
               istk(il1+3+i)=0
3094
3165
            endif
3095
3166
 132     continue
 
3167
         endif
3096
3168
         lstk(top+1)=sadr(il1+3+mn1)
3097
3169
      endif
3098
3170
      return
3158
3230
 10         continue
3159
3231
            lstk(top+1)=sadr(il1+3+mn1)
3160
3232
         endif
3161
 
      elseif(fin.eq.et) then
3162
 
         if(mn1.eq.0.or.mn2.eq.0) then
3163
 
            istk(il1)=1
3164
 
            istk(il1+1)=0
3165
 
            istk(il1+2)=0
3166
 
            istk(il1+3)=0
3167
 
            lstk(top+1)=sadr(il1+4)
3168
 
            return
3169
 
         endif
3170
 
         if(mn1.eq.1) then
3171
 
            i1=0
3172
 
            mn1=mn2
3173
 
         else
3174
 
            i1=1
3175
 
         endif
3176
 
         if(mn2.eq.1) then
3177
 
            i2=0
3178
 
            mn2=mn1
3179
 
         else
3180
 
            i2=1
3181
 
         endif
3182
 
         istk(il1)=4
3183
 
         do 20 i=0,mn1-1
3184
 
            e1=stk(l1+i)
3185
 
            e2=stk(l2+i)
3186
 
            if(e1.ne.0.0d0.and.e2.ne.0.0d0) then
3187
 
               istk(il1+3+i)=1
3188
 
            else
3189
 
               istk(il1+3+i)=0
3190
 
            endif
3191
 
 20      continue
3192
 
         lstk(top+1)=sadr(il1+3+mn1)
3193
 
         return
3194
 
      elseif(fin.eq.ou) then
3195
 
         if(mn1.eq.0.or.mn2.eq.0) then
3196
 
            istk(il1)=1
3197
 
            istk(il1+1)=0
3198
 
            istk(il1+2)=0
3199
 
            istk(il1+3)=0
3200
 
            lstk(top+1)=sadr(il1+4)
3201
 
            return
3202
 
         endif
3203
 
         if(mn1.eq.1) then
3204
 
            i1=0
3205
 
            mn1=mn2
3206
 
         else
3207
 
            i1=1
3208
 
         endif
3209
 
         if(mn2.eq.1) then
3210
 
            i2=0
3211
 
            mn2=mn1
3212
 
         else
3213
 
            i2=1
3214
 
         endif
3215
 
         do 30 i=0,mn1-1
3216
 
            e1=stk(l1+i)
3217
 
            e2=stk(l2+i)
3218
 
            if(e1.ne.0.0d0.or.e2.ne.0.0d0) then
3219
 
               istk(il1+3+i)=1
3220
 
            else
3221
 
               istk(il1+3+i)=0
3222
 
            endif
3223
 
 30      continue
 
3233
      else
 
3234
         if(mn1.eq.0.or.mn2.eq.0) then
 
3235
            istk(il1)=1
 
3236
            istk(il1+1)=0
 
3237
            istk(il1+2)=0
 
3238
            istk(il1+3)=0
 
3239
            lstk(top+1)=sadr(il1+4)
 
3240
            return
 
3241
         endif
 
3242
         if(mn1.eq.1) then
 
3243
            i1=0
 
3244
            mn1=mn2
 
3245
         else
 
3246
            i1=1
 
3247
         endif
 
3248
         if(mn2.eq.1) then
 
3249
            i2=0
 
3250
            mn2=mn1
 
3251
         else
 
3252
            i2=1
 
3253
         endif
 
3254
         if(mn1.ne.mn2) then
 
3255
            call error(60)
 
3256
            return
 
3257
         endif
 
3258
         if(fin.eq.ou) then
 
3259
            do 30 i=0,mn1-1
 
3260
               e1=stk(l1+i*i1)
 
3261
               e2=stk(l2+i*i2)
 
3262
               if(e1.ne.0.0d0.or.e2.ne.0.0d0) then
 
3263
                  istk(il1+3+i)=1
 
3264
               else
 
3265
                  istk(il1+3+i)=0
 
3266
               endif
 
3267
 30         continue
 
3268
         else
 
3269
            do 31 i=0,mn1-1
 
3270
               e1=stk(l1+i*i1)
 
3271
               e2=stk(l2+i*i2)
 
3272
               if(e1.ne.0.0d0.and.e2.ne.0.0d0) then
 
3273
                  istk(il1+3+i)=1
 
3274
               else
 
3275
                  istk(il1+3+i)=0
 
3276
               endif
 
3277
 31         continue
 
3278
         endif
3224
3279
         istk(il1)=4
3225
3280
         istk(il1+1)=max(m1,m2)
3226
3281
         istk(il1+2)=max(n1,n2)