← Back to branch summary

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

« back to all changes in this revision

Viewing changes to routines/interf/matelm.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/matelm.f
26
26
c      13   14   15   16   17   18  19-21 22   23   24   25   26
27
27
c     sqrt log   ^  sign clean floor ceil expm cumsum  cumprod testmatrix
28
28
c      27   28   29  30   31     32   33   34    35      36      37
29
 
c     isreal frexp
30
 
c       38     39
 
29
c     isreal frexp zeros tan  log1p imult  asin acos number_properties
 
30
c       38     39    40   41     42    43   44   45      46
 
31
c     nearfloat dsearch isequal
 
32
c       47        48      49
31
33
c!
32
34
c
33
35
      goto (10 ,15 ,20 ,25 ,30 ,35 ,40 ,45 ,50 ,60,
34
36
     1      61 ,62 ,70 ,72 ,71 ,90 ,91 ,105,110,110,
35
37
     2      110,130,140,150,160,170,180,190,200,210,
36
 
     3      220,37 ,39 ,173,46 ,47, 230,240, 250),fin
 
38
     3      220,37 ,39 ,173,46 ,47, 230,240,250,260,
 
39
     4      165,195,196,152,154,300,310,320,330     ),fin
37
40
 
38
41
 10   continue
39
42
      call intabs(id)
166
169
 150  continue
167
170
      call intcos(id)
168
171
      goto 900
 
172
c
 
173
c     asin
 
174
c
 
175
 152  continue
 
176
      call intasin(id)
 
177
      goto 900
 
178
c
 
179
c     acos
 
180
c
 
181
 154  continue
 
182
      call intacos(id)
 
183
      goto 900  
169
184
c     
170
185
c     atan
171
186
c     
172
187
 160  continue
173
188
      call intatan(id)
174
189
      goto 900
 
190
c
 
191
c     tan
 
192
c
 
193
 165  continue
 
194
      call inttan(id)
 
195
      goto 900
175
196
 
176
197
c     exp element wise
177
198
 170  continue
194
215
 190  continue
195
216
      call intlog(id)
196
217
      goto 900
 
218
c     
 
219
c     log1p
 
220
c     
 
221
 195  continue
 
222
      call intlog1p(id)
 
223
      goto 900
 
224
c     
 
225
c     imult
 
226
c     
 
227
 196  continue
 
228
      call intimult(id)
 
229
      goto 900
197
230
 
198
231
c     
199
232
c     ** non integer power of square matrices or  scalar^matrix
232
265
      call intfrexp(id)
233
266
      goto 900
234
267
c
 
268
c     zeros
 
269
c
 
270
 
 
271
 260  continue
 
272
      call intzeros(id)
 
273
      goto 900
 
274
 
 
275
c
 
276
c     number_properties
 
277
c
 
278
 300  call  intnbprop(id)
 
279
      go to 900
 
280
c
 
281
c     number_properties
 
282
c
 
283
 310  call  intnearfl(id)
 
284
      go to 900
 
285
 
 
286
c
 
287
c     dsearch
 
288
c
 
289
 320  call intdsearch(id)
 
290
      goto 900
 
291
c
 
292
c     isequal
 
293
c
 
294
 330  call intisequal(id)
 
295
      goto 900
 
296
c
235
297
 900  return
236
298
      end
237
299
 
238
300
      subroutine intrand(fname,id)
239
301
c     -------------------------------
240
302
      character*(*) fname
241
 
c     Interface for rand function 
 
303
c     Interface for rand function
242
304
      INCLUDE '../stack.h'
243
 
      double precision s,sr,si
 
305
      double precision s,sr,si,r
 
306
      save             si, r
244
307
      integer id(nsiz),tops,topk,name(nlgh)
245
308
      double precision urand
246
309
      logical checkrhs,checklhs,getsmat,getscalar,cremat,getmat
247
310
      logical cresmat2
248
311
      integer gettype
249
312
      character*(20) randtype
 
313
      logical phase
 
314
      save    phase
250
315
      integer iadr,sadr
 
316
      data    phase /.true./
251
317
c
252
318
      iadr(l)=l+l-1
253
319
c      sadr(l)=(l/2)+1
391
457
      if (.not.cremat(fname,top,itres,m,n,lr,lc)) return
392
458
      if ( m*n .ne. 0 ) then 
393
459
         if ( ran(2).eq.0 ) then 
394
 
            do 76 j = 0, (itres+1)*m*n-1
 
460
*           U(0,1) random numbers
 
461
            do j = 0, (itres+1)*m*n-1
395
462
               stk(lr+j) = urand(ran(1))
396
 
 76         continue
 
463
            enddo
397
464
         elseif (ran(2).eq.1) then 
398
 
            do 77 j = 0, m*n-1
399
 
 75            sr=2.0d+0*urand(ran(1))-1.0d+0
400
 
               si=2.0d+0*urand(ran(1))-1.0d+0
401
 
               t = sr*sr+si*si
402
 
               if (t .gt. 1.0d+0) go to 75
403
 
               stk(lr+j) = sr*sqrt(-2.0d+0*log(t)/t)
404
 
 77         continue
 
465
*           N(0,1) random numbers (modified by Bruno 11/10/2001
 
466
*                  to use si*r and to correct the bug in the
 
467
*                  complex case)
 
468
            do j = 0, (itres+1)*m*n-1
 
469
               if (phase) then
 
470
 75               sr=2.d0*urand(ran(1)) - 1.d0
 
471
                  si=2.d0*urand(ran(1)) - 1.d0
 
472
                  t = sr*sr + si*si
 
473
                  if (t .gt. 1.d0) go to 75
 
474
                  r = sqrt(-2.d0*log(t)/t)
 
475
                  stk(lr+j) = sr*r
 
476
               else
 
477
                  stk(lr+j) = si*r
 
478
               endif
 
479
               phase = .not. phase
 
480
            enddo
405
481
         endif
406
482
      endif
407
483
C     switching back to the default randvalue
531
607
c     ------------max of each column of a 
532
608
         if (.not.cremat(fname,topk,0,1,n,lr,lir)) return
533
609
         if (.not.cremat(fname,topk+1,0,1,n,lkr,lkir)) return
534
 
         do 20 j=0,n-1
535
 
            stk(lr+j)=stk(lr1+m*j)
536
 
            stk(lkr+j)=1
537
 
            if(fin.eq.17) then 
538
 
               do 15 i=1,m-1
539
 
                  x1=stk(lr1+i+m*j)
540
 
                  if(x1.lt.stk(lr+j).or.isanan(x1).eq.1) then 
541
 
                     stk(lr+j)=x1
542
 
                     stk(lkr+j)=i+1
543
 
                  endif
544
 
 15           continue
545
 
            else
546
 
               do 16 i=1,m-1
547
 
                  x1=stk(lr1+i+m*j)
548
 
                  if(x1.gt.stk(lr+j).or.isanan(x1).eq.1) then 
549
 
                     stk(lr+j)=x1
550
 
                     stk(lkr+j)=i+1
551
 
                  endif
552
 
 16            continue
553
 
            endif
554
 
 20     continue
 
610
         if(fin.eq.17) then
 
611
c     .    min
 
612
            do 15 j=0,n-1
 
613
               k=idmin(m,stk(lr1+m*j),1)
 
614
               stk(lr+j)=stk(lr1+m*j+k-1)
 
615
               stk(lkr+j)=k
 
616
 15         continue
 
617
         else
 
618
c     .    max
 
619
            do 16 j=0,n-1
 
620
               k=idmax(m,stk(lr1+m*j),1)
 
621
               stk(lr+j)=stk(lr1+m*j+k-1)
 
622
               stk(lkr+j)=k
 
623
 16         continue
 
624
         endif
555
625
         call copyobj(fname,topk,topk-rhs+1)
556
626
         if (lhs.eq.2) then 
557
627
            call copyobj(fname,topk+1,topk-rhs+2)
561
631
      else if ( type(1:1).eq.'c') then       
562
632
         if (.not.cremat(fname,topk,0,m,1,lr,lir)) return
563
633
         if (.not.cremat(fname,topk+1,0,m,1,lkr,lkir)) return
564
 
         do 30 j=0,m-1
565
 
            stk(lr+j)=stk(lr1+j)
566
 
            stk(lkr+j)=1
567
 
            if(fin.eq.17) then 
568
 
               do 25 i=1,n-1
569
 
                  x1=stk(lr1+j+m*i)
570
 
                  if (x1.lt.stk(lr+j).or.isanan(x1).eq.1) then 
571
 
                     stk(lr+j)=x1
572
 
                     stk(lkr+j)=i+1
573
 
                  endif
574
 
 25            continue
575
 
            else
576
 
               do 26 i=1,n-1
577
 
                  x1=stk(lr1+j+m*i)
578
 
                  if (x1.gt.stk(lr+j).or.isanan(x1).eq.1) then 
579
 
                     stk(lr+j)=x1
580
 
                     stk(lkr+j)=i+1
581
 
                  endif
582
 
 26           continue
583
 
            endif
584
 
 30      continue
 
634
         if(fin.eq.17) then
 
635
c     .    min
 
636
            do 25 j=0,m-1
 
637
               k=idmin(n,stk(lr1+j),m)
 
638
               stk(lr+j)=stk(lr1+j+(k-1)*m)
 
639
               stk(lkr+j)=k
 
640
 25         continue
 
641
         else
 
642
c     .    max
 
643
            do 26 j=0,m-1
 
644
               k=idmax(n,stk(lr1+j),m)
 
645
               stk(lr+j)=stk(lr1+j+(k-1)*m)
 
646
               stk(lkr+j)=k
 
647
 26         continue
 
648
         endif
585
649
         call copyobj(fname,topk,topk-rhs+1)
586
650
         if (lhs.eq.2) then 
587
651
            call copyobj(fname,topk+1,topk-rhs+2)
588
652
         endif
589
653
         top=topk-rhs+lhs            
590
654
c     ----- general maxi or mini 
591
 
      else if ( type(1:1).eq.'g') then 
592
 
         k=lr1
593
 
         ki=lr1
594
 
         x1=stk(k)
 
655
      else if ( type(1:1).eq.'g'.or.type(1:1).eq.'*') then 
 
656
         if (rhs.eq.2) topk=top
 
657
 
 
658
         x1=stk(lr1)
 
659
         k=1
595
660
         if(fin.eq.17) then 
596
 
c     mini
597
 
            do 41 i=2,m*n
598
 
               lr1=lr1+1
599
 
               if(stk(lr1).lt.x1.or.isanan(stk(lr1)).eq.1) then 
600
 
                  k=lr1
601
 
                  x1=stk(k)
602
 
               endif
603
 
 41         continue
604
 
c     maxi
 
661
c     .     mini
 
662
            k=idmin(m*n,stk(lr1),1)
605
663
         else
606
 
            do 42 i=2,m*n
607
 
               lr1=lr1+1
608
 
               if(stk(lr1).gt.x1.or.isanan(stk(lr1)).eq.1) then 
609
 
                  k=lr1
610
 
                  x1=stk(k)
611
 
               endif
612
 
 42         continue
 
664
c     .     maxi
 
665
            k=idmax(m*n,stk(lr1),1)
613
666
         endif
 
667
         x1=stk(lr1-1+k)
614
668
C     return the max or min 
615
669
         if (.not.cremat(fname,topk,0,1,1,l1,li1)) return
616
670
         stk(l1)=x1
618
672
c     for matrices 
619
673
         if(lhs.eq.2) then 
620
674
            top=topk+1
621
 
            k=k-ki+1
622
675
            if(m.eq.1.or.n.eq.1) then 
623
676
               if (.not.cremat(fname,top,0,1,1,lr1,lc1)) return
624
677
               stk(lr1)=dble(k)
798
851
c     WARNING : argument of this interface may be passed by reference
799
852
      INCLUDE '../stack.h'
800
853
      integer id(nsiz)
801
 
      double precision pythag
 
854
      double precision dlapy2
802
855
      logical ref
803
856
      integer head
804
857
      integer iadr,sadr
864
917
      else
865
918
         k1=l1+mn
866
919
         do 13 i=1,mn
867
 
            stk(lr1+i)=pythag(stk(l1+i),stk(k1+i))
 
920
            stk(lr1+i)=dlapy2(stk(l1+i),stk(k1+i))
868
921
 13      continue
869
922
         istk(ilr+3)=0
870
923
      endif
1085
1138
c     WARNING : argument of this interface may be passed by reference
1086
1139
      INCLUDE '../stack.h'
1087
1140
      integer id(nsiz)
1088
 
      double precision round
1089
1141
      logical ref
1090
1142
      integer iadr,sadr
1091
1143
c     
1113
1165
      n=istk(il+2)
1114
1166
      it=istk(il+3)
1115
1167
 
1116
 
      if(istk(il).eq.1) then
 
1168
      if(istk(il).eq.8) then
 
1169
         return
 
1170
      elseif(istk(il).eq.1) then
1117
1171
         mn=m*n
1118
1172
         l=sadr(il+4)
1119
1173
         lr=sadr(ilr+4)
1146
1200
      if(it.eq.1) mn=2*mn
1147
1201
      do 10 i=1,mn
1148
1202
       i1=i-1
1149
 
       stk(lr+i1)=round(stk(l+i1))
 
1203
       stk(lr+i1)=anint(stk(l+i1))
1150
1204
 10   continue
1151
1205
      lstk(top+1)=lr+mn
1152
1206
      return
1182
1236
      m=istk(il+1)
1183
1237
      n=istk(il+2)
1184
1238
      it=istk(il+3)
1185
 
 
1186
 
      if(istk(il).eq.1) then
 
1239
      if(istk(il).eq.8) then
 
1240
         return
 
1241
      elseif(istk(il).eq.1) then
1187
1242
         mn=m*n
1188
1243
         l=sadr(il+4)
1189
1244
         lr=sadr(ilr+4)
1226
1281
c     WARNING : argument of this interface may be passed by reference
1227
1282
      INCLUDE '../stack.h'
1228
1283
      integer id(nsiz)
1229
 
      double precision t,t1,round
 
1284
      double precision t,t1
1230
1285
      logical ref
1231
1286
      integer iadr,sadr
1232
1287
c     
1253
1308
      m=istk(il+1)
1254
1309
      n=istk(il+2)
1255
1310
      it=istk(il+3)
1256
 
 
1257
 
      if(istk(il).eq.1) then
 
1311
      if(istk(il).eq.8) then
 
1312
         return
 
1313
      elseif(istk(il).eq.1) then
1258
1314
         mn=m*n
1259
1315
         l=sadr(il+4)
1260
1316
         lr=sadr(ilr+4)
1293
1349
c     WARNING : argument of this interface may be passed by reference
1294
1350
      INCLUDE '../stack.h'
1295
1351
      integer id(nsiz)
1296
 
      double precision t,t1,round
 
1352
      double precision t,t1
1297
1353
      logical ref
1298
1354
      integer iadr,sadr
1299
1355
c     
1320
1376
      m=istk(il+1)
1321
1377
      n=istk(il+2)
1322
1378
      it=istk(il+3)
1323
 
 
1324
 
      if(istk(il).eq.1) then
 
1379
      if(istk(il).eq.8) then
 
1380
         return
 
1381
      elseif(istk(il).eq.1) then
1325
1382
         mn=m*n
1326
1383
         l=sadr(il+4)
1327
1384
         lr=sadr(ilr+4)
1457
1514
               istk(ilr+3)=0
1458
1515
               stk(lr) = n
1459
1516
               lstk(top+1)=lr+1
1460
 
            elseif(sel.eq.3) then
 
1517
            elseif(sel.eq.0) then
1461
1518
               istk(ilr+2)=1
1462
1519
               istk(ilr+3)=0
1463
1520
               stk(lr) = m*n
1809
1866
      lstk(top+1)=l1+mn*(it+1)
1810
1867
      if(sel.eq.0) then
1811
1868
c     op(a) <=> op(a,'*')
1812
 
         call cupro(mn,stk(l1))
1813
 
         if(it.eq.1) call cupro(mn,stk(l1+mn))
 
1869
         if(it.eq.0) then
 
1870
            call cupro(mn,stk(l1))
 
1871
         else
 
1872
            call cuproi(mn,stk(l1),stk(l1+mn))
 
1873
         endif
1814
1874
      elseif(sel.eq.1) then
1815
1875
c     op(a,'r')  <=>  op(a,1)
1816
1876
         if(it.eq.0) then
2052
2112
         elseif(istk(il+6).eq.col) then
2053
2113
            sel=2
2054
2114
         elseif(istk(il+6).eq.star) then
2055
 
            sel=3
 
2115
            sel=0
2056
2116
         else
2057
2117
            err=2
2058
2118
            call error(44)
2628
2688
      return
2629
2689
      end
2630
2690
 
 
2691
      subroutine intzeros(id)
 
2692
      INCLUDE '../stack.h'
 
2693
      integer id(nsiz),name(nlgh)
 
2694
 
 
2695
      integer tops
 
2696
      double precision s
 
2697
      integer iadr,sadr
 
2698
c
 
2699
      iadr(l)=l+l-1
 
2700
      sadr(l)=(l/2)+1
 
2701
c
 
2702
      if (lhs .ne. 1) then
 
2703
         call error(41)
 
2704
         return
 
2705
      endif
 
2706
      if(rhs.gt.2) then
 
2707
         call setfunnam(ids(1,pt+1),'%hm_zeros',9)
 
2708
         fun=-1
 
2709
c         call error(42)
 
2710
         return
 
2711
      endif
 
2712
 
 
2713
      tops=top
 
2714
c
 
2715
 
 
2716
      if(rhs.le.0) then
 
2717
c     ones sans argument
 
2718
         top=top+1
 
2719
         m=1
 
2720
         n=1
 
2721
      elseif(rhs.eq.1) then
 
2722
         il=iadr(lstk(top))
 
2723
         if(abs(istk(il)).gt.10) then
 
2724
            call funnam(ids(1,pt+1),'zeros',il)
 
2725
            fun=-1
 
2726
            return
 
2727
         endif
 
2728
         if(istk(il).lt.0) il=iadr(istk(il+1))
 
2729
         m=istk(il+1)
 
2730
         n=istk(il+2)
 
2731
c     ones(matrice)
 
2732
      elseif(rhs.eq.2) then
 
2733
c     ones(m,n)
 
2734
         il=iadr(lstk(top))
 
2735
         if(istk(il).lt.0) il=iadr(istk(il+1))
 
2736
         if(istk(il).ne.1) then
 
2737
            err=1
 
2738
            call  error(53)
 
2739
            return
 
2740
         endif
 
2741
         if(istk(il+3).ne.0) then
 
2742
            err=1
 
2743
            call  error(52)
 
2744
            return
 
2745
         endif
 
2746
         if(istk(il+1)*istk(il+2).ne.1) then
 
2747
            err=1
 
2748
            call error(89)
 
2749
            return
 
2750
         endif
 
2751
         n=max(int(stk(sadr(il+4))),0)
 
2752
c
 
2753
         top=top-1
 
2754
         il=iadr(lstk(top))
 
2755
         if(istk(il).lt.0) il=iadr(istk(il+1))
 
2756
         if(istk(il).ne.1) then
 
2757
            err=1
 
2758
            call  error(53)
 
2759
            return
 
2760
         endif
 
2761
         if(istk(il+3).ne.0) then
 
2762
            err=1
 
2763
            call  error(52)
 
2764
            return
 
2765
         endif
 
2766
         if(istk(il+1)*istk(il+2).ne.1) then
 
2767
            err=1
 
2768
            call error(89)
 
2769
            return
 
2770
         endif
 
2771
         m=max(int(stk(sadr(il+4))),0)
 
2772
      endif
 
2773
c
 
2774
      mn=m*n
 
2775
      if(m.eq.0) n=0
 
2776
      if(n.eq.0) m=0
 
2777
 
 
2778
      il=iadr(lstk(top))
 
2779
      l=sadr(il+4)
 
2780
 
 
2781
c     to avoid integer overflow
 
2782
      s=l+dble(m)*dble(n)- lstk(bot)
 
2783
      if(s.gt.0.0d0) then
 
2784
         err=s
 
2785
         call error(17)
 
2786
         return
 
2787
      endif
 
2788
      istk(il)=1
 
2789
      istk(il+1)=m
 
2790
      istk(il+2)=n
 
2791
      istk(il+3)=0
 
2792
      lstk(top+1)=l+mn
 
2793
      if(mn.eq.0) return
 
2794
      call dset(mn,0.0d+0,stk(l),1)
 
2795
      return
 
2796
      end
 
2797
 
2631
2798
      subroutine intsort(id)
2632
2799
      INCLUDE '../stack.h'
2633
2800
      integer id(nsiz)
2858
3025
      la=sadr(il+4)
2859
3026
      mna=ma*na
2860
3027
c
 
3028
      l=sadr(ilr+4)
 
3029
      l1=l+mnb*mna*(max(itb,ita)+1)
 
3030
      lstk(top+1)=l1
 
3031
      lw=l1
 
3032
c
2861
3033
      if(fin.eq.19) goto 115
2862
 
      if(fin.eq.20) goto 111
2863
 
      l=la
2864
 
      mn=mna
2865
 
      it=ita
 
3034
      if(fin.eq.20) then
 
3035
         if(refb) then
 
3036
            err=lw+mnb*(itb+1)-lstk(bot)
 
3037
            if(err.gt.0) then
 
3038
               call error(17)
 
3039
               return
 
3040
            endif
 
3041
            call dcopy(mnb*(itb+1),stk(lb),1,stk(lw),1)
 
3042
         endif
 
3043
         lb=lw
 
3044
         lw=lb+mnb*(itb+1)
 
3045
         l=lb
 
3046
         mn=mnb
 
3047
         it=itb
 
3048
      elseif(fin.eq.21) then
 
3049
         if(refa) then
 
3050
            err=lw+mna*(ita+1)-lstk(bot)
 
3051
            if(err.gt.0) then
 
3052
               call error(17)
 
3053
               return
 
3054
            endif
 
3055
            call dcopy(mna*(ita+1),stk(la),1,stk(lw),1)
 
3056
         endif
 
3057
         la=lw
 
3058
         lw=la+mna*(ita+1)
 
3059
         l=la
 
3060
         mn=mna
 
3061
         it=ita
 
3062
      endif
 
3063
 
2866
3064
  111 if(it.eq.1) goto 113
2867
3065
      do 112 k=1,mn
2868
3066
         lk=l+k-1
2887
3085
 114  continue
2888
3086
c
2889
3087
 115  continue
2890
 
      l=sadr(ilr+4)
2891
 
      l1=l+mnb*mna*(max(itb,ita)+1)
2892
 
      lstk(top+1)=l1
 
3088
 
2893
3089
c
2894
3090
c move a and b if necessary
2895
 
      lw=l1
2896
3091
      if(.not.refb) then
2897
3092
         lw=lw+mnb*(itb+1)
2898
3093
         err=lw-lstk(bot)
2976
3171
 
2977
3172
c changement de dimension d'une matrice
2978
3173
      il=iadr(lstk(top+1-rhs))
2979
 
 
2980
3174
      if(abs(istk(il)).eq.5.or.abs(istk(il)).eq.6) then
2981
3175
         top=tops
2982
3176
         call ref2val
2988
3182
 
2989
3183
      ityp=abs(istk(il))
2990
3184
      if(ityp.ne.1.and.ityp.ne.2
2991
 
     $     .and.ityp.ne.4.and.ityp.ne.10) then
 
3185
     $     .and.ityp.ne.4.and.ityp.ne.9.and.ityp.ne.10) then
2992
3186
         top=tops
2993
3187
         call ref2val
2994
3188
         call funnam(ids(1,pt+1),'matrix',iadr(lstk(top-rhs+1)))
3089
3283
         call error(42)
3090
3284
         return
3091
3285
      endif
 
3286
      if(m.eq.0.or.n.eq.0) then
 
3287
         if(mn.ne.0) then
 
3288
            call error(60)
 
3289
            return
 
3290
         endif
 
3291
         return
 
3292
      endif
3092
3293
      if(m.eq.-1) m=mn/n
3093
3294
      if(n.eq.-1) n=mn/m
3094
3295
 
3303
3504
            do 11 i=0,mn-1
3304
3505
               sr=stk(l+i)
3305
3506
               si=stk(l+mn+i)
3306
 
               if(abs(sinh(sr)).eq.1.0d+0.and.si.eq.0.0d+0) then
 
3507
c correction (Bruno) : les points singuliers de atan sont +- i
 
3508
               if ( sr .eq. 0.d0 .and. abs(si) .eq. 1.d0 ) then
3307
3509
                  if(ieee.eq.0) then
3308
3510
                     call error(32)
3309
3511
                     return
3669
3871
         do 10 i=0,mn-1
3670
3872
            if(stk(l+i).lt.0.0d+0) then
3671
3873
               itr=1
3672
 
               goto 20
 
3874
c               goto 20
3673
3875
            elseif(stk(l+i).eq.0.0d+0) then
3674
3876
               if(ieee.eq.0) then
3675
3877
                  call error(32)
3677
3879
               elseif(ieee.eq.1) then
3678
3880
                  call msgs(64)
3679
3881
               endif
3680
 
               goto 20
 
3882
c               goto 20
3681
3883
            endif
3682
3884
 10      continue
3683
3885
 
3684
3886
 20      if(itr.eq.0) then
3685
 
c     .     argument is a real positive matrix
 
3887
c     .     argument is a real positive matrix with entries >= 0
3686
3888
            do 193 i=0,mn-1
 
3889
               if(stk(l+i).eq.0.0d+0) then
 
3890
                  if(ieee.eq.0) then
 
3891
                     call error(32)
 
3892
                     return
 
3893
                  elseif(ieee.eq.1) then
 
3894
                     call msgs(64)
 
3895
                  endif
 
3896
               endif
3687
3897
               stk(lr+i)=log(stk(l+i))
3688
3898
 193        continue
3689
3899
         else
3690
 
c     .     argument is a real matrix with negative entries
 
3900
c     .     argument is a real matrix with  at least one entry < 0     
3691
3901
            err=lr+2*mn-lstk(bot)
3692
3902
            if(err.gt.0) then
3693
3903
               call error(17)
3695
3905
            endif
3696
3906
            lstk(top+1)=lr+2*mn
3697
3907
            do 194 i=0,mn-1
 
3908
               if(stk(l+i).eq.0.0d+0) then
 
3909
                  if(ieee.eq.0) then
 
3910
                     call error(32)
 
3911
                     return
 
3912
                  elseif(ieee.eq.1) then
 
3913
                     call msgs(64)
 
3914
                  endif
 
3915
               endif
3698
3916
               call wlog(stk(l+i),0.0d+0,stk(lr+i),stk(lr+mn+i))
3699
3917
 194        continue
3700
3918
            istk(ilr+3)=itr
3704
3922
         do 195 i=0,mn-1
3705
3923
            sr=stk(l+i)
3706
3924
            si=stk(l+mn+i)
3707
 
            if(sr*sr+si*si.eq.0.0d+0) then
 
3925
            if(sr.eq.0d0 .and. si.eq.0.0d+0) then  
3708
3926
               if(ieee.eq.0) then
3709
3927
                  call error(32)
3710
3928
                  return
3723
3941
      INCLUDE '../stack.h'
3724
3942
      integer id(nsiz)
3725
3943
 
3726
 
      double precision sr,si,pythag
 
3944
      double precision sr,si,dlapy2
3727
3945
      integer iadr,sadr
3728
3946
c
3729
3947
      iadr(l)=l+l-1
3783
4001
 10      continue
3784
4002
      else
3785
4003
         do 20 i=0,mn-1
3786
 
            sr=pythag(stk(l+i),stk(l+mn+i))
 
4004
            sr=dlapy2(stk(l+i),stk(l+mn+i))
3787
4005
            if(sr.eq.0) then
3788
4006
               stk(lr+i)=0.0d0
3789
4007
               stk(lr+mn+i)=0.0d0
3800
4018
      INCLUDE '../stack.h'
3801
4019
      integer id(nsiz)
3802
4020
 
3803
 
      double precision sr,si,pythag,epsa,epsr,norm,eps
 
4021
      double precision sr,si,dlapy2,epsa,epsr,norm,eps
3804
4022
      double precision dasum,wasum
3805
4023
      integer iadr,sadr
3806
4024
c
4093
4311
      endif
4094
4312
      return
4095
4313
      end
 
4314
 
4096
4315
      subroutine intfrexp(id)
4097
4316
c     WARNING : argument of this interface may be passed by reference
4098
4317
      INCLUDE '../stack.h'
4099
4318
      integer id(nsiz)
4100
 
      double precision t,t1,round
 
4319
      double precision t,t1
4101
4320
      logical ref
4102
4321
      integer iadr,sadr
4103
4322
c     
4159
4378
      call vfrexp(mn,stk(l),1,stk(lr),1,stk(lr1),1)
4160
4379
      return
4161
4380
      end
 
4381
      subroutine inttan(id)
 
4382
      INCLUDE '../stack.h'
 
4383
c     
 
4384
c     interface for tan
 
4385
c     
 
4386
      integer id(nsiz)
 
4387
 
 
4388
      integer iadr,sadr
 
4389
     
 
4390
      iadr(l)=l+l-1
 
4391
      sadr(l)=(l/2)+1
 
4392
 
 
4393
      if (lhs .ne. 1) then
 
4394
         call error(41)
 
4395
         return
 
4396
      endif
 
4397
      if(rhs .ne. 1) then
 
4398
         call error(42)
 
4399
         return
 
4400
      endif
 
4401
 
 
4402
      il=iadr(lstk(top))
 
4403
 
 
4404
      if(abs(istk(il)).ne.1) then
 
4405
         call funnam(ids(1,pt+1),'tan',iadr(lstk(top)))
 
4406
         fun=-1
 
4407
         return
 
4408
      endif
 
4409
 
 
4410
      if(istk(il).lt.0) then
 
4411
c        argument is passed by reference
 
4412
         ilr=il
 
4413
         il=iadr(istk(il+1))
 
4414
         mn=istk(il+1)*istk(il+2)
 
4415
         it=istk(il+3)
 
4416
         l=sadr(il+4)
 
4417
         lr=sadr(ilr+4)
 
4418
         err=lr+mn*(it+1)-lstk(bot)
 
4419
         if(err.gt.0) then
 
4420
            call error(17)
 
4421
            return
 
4422
         endif
 
4423
         call icopy(4,istk(il),1,istk(ilr),1)
 
4424
         lstk(top+1)=lr+mn*(it+1)
 
4425
      else
 
4426
         mn=istk(il+1)*istk(il+2)
 
4427
         it=istk(il+3)
 
4428
         l=sadr(il+4)
 
4429
         lr=l
 
4430
      endif
 
4431
 
 
4432
      if(mn.eq.0) return
 
4433
 
 
4434
      if(it.eq.0) then
 
4435
         do i=0,mn-1
 
4436
            stk(lr+i)=tan(stk(l+i))
 
4437
         enddo
 
4438
      else
 
4439
         do i=0,mn-1
 
4440
            call wtan(stk(l+i),stk(l+mn+i),stk(lr+i),stk(lr+mn+i))
 
4441
         enddo
 
4442
      endif
 
4443
      return
 
4444
      end
 
4445
 
 
4446
 
 
4447
      subroutine intimult(id)
 
4448
      INCLUDE '../stack.h'
 
4449
      integer id(nsiz)
 
4450
c
 
4451
c     interface for imult : multiplication by i
 
4452
c
 
4453
      double precision sr,si
 
4454
      integer iadr,sadr
 
4455
 
 
4456
      iadr(l)=l+l-1
 
4457
      sadr(l)=(l/2)+1
 
4458
 
 
4459
      if (lhs .ne. 1) then
 
4460
         call error(41)
 
4461
         return
 
4462
      endif
 
4463
      if(rhs .ne. 1) then
 
4464
         call error(42)
 
4465
         return
 
4466
      endif
 
4467
 
 
4468
      il=iadr(lstk(top))
 
4469
      
 
4470
      if(abs(istk(il)).ne.1) then
 
4471
         call funnam(ids(1,pt+1),'imult',iadr(lstk(top)))
 
4472
         fun=-1
 
4473
         return
 
4474
      endif
 
4475
 
 
4476
      if(istk(il).lt.0) then
 
4477
c        argument is passed by reference
 
4478
         ilr=il
 
4479
         il=iadr(istk(il+1))
 
4480
         mn=istk(il+1)*istk(il+2)
 
4481
         it=istk(il+3)
 
4482
         l=sadr(il+4)
 
4483
         lr=sadr(ilr+4)
 
4484
         err=lr+mn*(it+1)-lstk(bot)
 
4485
         if(err.gt.0) then
 
4486
            call error(17)
 
4487
            return
 
4488
         endif
 
4489
         call icopy(4,istk(il),1,istk(ilr),1)
 
4490
         lstk(top+1)=lr+mn*(it+1)
 
4491
      else
 
4492
         mn=istk(il+1)*istk(il+2)
 
4493
         it=istk(il+3)
 
4494
         l=sadr(il+4)
 
4495
         lr=l
 
4496
         ilr=il
 
4497
      endif
 
4498
      
 
4499
      if(mn.eq.0) return
 
4500
 
 
4501
      if(it.eq.0) then
 
4502
c        argument is real but result is complex
 
4503
         err=lr+2*mn-lstk(bot)
 
4504
         if(err.gt.0) then
 
4505
            call error(17)
 
4506
            return
 
4507
         endif
 
4508
         lstk(top+1)=lr+2*mn
 
4509
         do i=0,mn-1
 
4510
            stk(lr+mn+i) = stk(l+i)
 
4511
            stk(lr+i) = 0.d0
 
4512
         enddo
 
4513
         istk(ilr+3)=1
 
4514
      else
 
4515
c        argument is complex
 
4516
         do i=0,mn-1
 
4517
            sr=stk(l+i)
 
4518
            si=stk(l+mn+i)
 
4519
            stk(lr+i) = -si
 
4520
            stk(lr+i+mn) = sr
 
4521
         enddo
 
4522
      endif
 
4523
      return
 
4524
      end
 
4525
 
 
4526
 
 
4527
      subroutine intlog1p(id)
 
4528
c
 
4529
c     interface for log1p function (log(1+x))
 
4530
c     rmk : don't work in complex
 
4531
c
 
4532
      INCLUDE '../stack.h'
 
4533
      integer id(nsiz)
 
4534
 
 
4535
      double precision logp1
 
4536
      external         logp1
 
4537
 
 
4538
      integer iadr,sadr
 
4539
c
 
4540
      iadr(l)=l+l-1
 
4541
      sadr(l)=(l/2)+1
 
4542
c
 
4543
      if (lhs .ne. 1) then
 
4544
         call error(41)
 
4545
         return
 
4546
      endif
 
4547
      if(rhs .ne. 1) then
 
4548
         call error(42)
 
4549
         return
 
4550
      endif
 
4551
 
 
4552
      il=iadr(lstk(top))
 
4553
      
 
4554
      if(abs(istk(il)).ne.1) then
 
4555
         call funnam(ids(1,pt+1),'log1p',iadr(lstk(top)))
 
4556
         fun=-1
 
4557
         return
 
4558
      endif
 
4559
 
 
4560
      if(istk(il).lt.0) then
 
4561
c        argument is passed by reference
 
4562
         ilr=il
 
4563
         il=iadr(istk(il+1))
 
4564
         mn=istk(il+1)*istk(il+2)
 
4565
         it=istk(il+3)
 
4566
         l=sadr(il+4)
 
4567
         lr=sadr(ilr+4)
 
4568
         err=lr+mn*(it+1)-lstk(bot)
 
4569
         if(err.gt.0) then
 
4570
            call error(17)
 
4571
            return
 
4572
         endif
 
4573
         call icopy(4,istk(il),1,istk(ilr),1)
 
4574
         lstk(top+1)=lr+mn*(it+1)
 
4575
      else
 
4576
         mn=istk(il+1)*istk(il+2)
 
4577
         it=istk(il+3)
 
4578
         l=sadr(il+4)
 
4579
         lr=l
 
4580
         ilr=il
 
4581
      endif
 
4582
      
 
4583
      if(mn.eq.0) return
 
4584
 
 
4585
      if(it.eq.0) then
 
4586
         itr=0
 
4587
         do i=0,mn-1
 
4588
            if(stk(l+i).le.-1.0d+0) then
 
4589
               if(ieee.eq.0) then
 
4590
                  call error(32)
 
4591
                  return
 
4592
               elseif(ieee.eq.1) then
 
4593
                  call msgs(64)
 
4594
               endif
 
4595
               goto 20
 
4596
            endif
 
4597
         enddo
 
4598
 
 
4599
 20      continue
 
4600
 
 
4601
         do i=0,mn-1
 
4602
            stk(lr+i)=logp1(stk(l+i))
 
4603
         enddo
 
4604
      else
 
4605
c        complex case : message "not implemented in scilab ..."
 
4606
         call error(43)
 
4607
         return
 
4608
      endif
 
4609
      end
 
4610
 
 
4611
 
 
4612
      subroutine intasin(id)
 
4613
      INCLUDE '../stack.h'
 
4614
      integer id(nsiz)
 
4615
c
 
4616
c     interface for the arcsin function
 
4617
c
 
4618
      double precision sr,si
 
4619
      integer iadr,sadr
 
4620
c
 
4621
      iadr(l)=l+l-1
 
4622
      sadr(l)=(l/2)+1
 
4623
c
 
4624
      if (lhs .ne. 1) then
 
4625
         call error(41)
 
4626
         return
 
4627
      endif
 
4628
      if(rhs .ne. 1) then
 
4629
         call error(42)
 
4630
         return
 
4631
      endif
 
4632
 
 
4633
      il=iadr(lstk(top))
 
4634
      
 
4635
      if(abs(istk(il)).ne.1) then
 
4636
         call funnam(ids(1,pt+1),'asin',iadr(lstk(top)))
 
4637
         fun=-1
 
4638
         return
 
4639
      endif
 
4640
 
 
4641
      if(istk(il).lt.0) then
 
4642
c        argument is passed by reference
 
4643
         ilr=il
 
4644
         il=iadr(istk(il+1))
 
4645
         mn=istk(il+1)*istk(il+2)
 
4646
         it=istk(il+3)
 
4647
         l=sadr(il+4)
 
4648
         lr=sadr(ilr+4)
 
4649
         err=lr+mn*(it+1)-lstk(bot)
 
4650
         if(err.gt.0) then
 
4651
            call error(17)
 
4652
            return
 
4653
         endif
 
4654
         call icopy(4,istk(il),1,istk(ilr),1)
 
4655
         lstk(top+1)=lr+mn*(it+1)
 
4656
      else
 
4657
         mn=istk(il+1)*istk(il+2)
 
4658
         it=istk(il+3)
 
4659
         l=sadr(il+4)
 
4660
         lr=l
 
4661
         ilr=il
 
4662
      endif
 
4663
      
 
4664
      if(mn.eq.0) return
 
4665
 
 
4666
      if(it.eq.0) then
 
4667
         itr=0
 
4668
         do i=0,mn-1
 
4669
            if(abs(stk(l+i)).gt.1.0d+0) then
 
4670
               itr=1
 
4671
               goto 20
 
4672
            endif
 
4673
         enddo
 
4674
 
 
4675
 20      continue
 
4676
         if(itr.eq.0) then
 
4677
c           argument is a real positive matrix with entries in [-1,1]
 
4678
            do i=0,mn-1
 
4679
               stk(lr+i)=asin(stk(l+i))
 
4680
            enddo
 
4681
         else
 
4682
c           argument is a real matrix with some entries outside [-1,1]
 
4683
            err=lr+2*mn-lstk(bot)
 
4684
            if(err.gt.0) then
 
4685
               call error(17)
 
4686
               return
 
4687
            endif
 
4688
            lstk(top+1)=lr+2*mn
 
4689
            do i=0,mn-1
 
4690
               call wasin(stk(l+i),0.0d+0,stk(lr+i),stk(lr+mn+i))
 
4691
            enddo
 
4692
            istk(ilr+3)=itr
 
4693
         endif
 
4694
      else
 
4695
c        argument is a complex matrix
 
4696
         do i=0,mn-1
 
4697
            call wasin(stk(l+i),stk(l+mn+i),stk(lr+i),stk(lr+i+mn))
 
4698
         enddo
 
4699
      endif
 
4700
      return
 
4701
      end
 
4702
 
 
4703
      subroutine intacos(id)
 
4704
      INCLUDE '../stack.h'
 
4705
      integer id(nsiz)
 
4706
 
 
4707
      double precision sr,si
 
4708
      integer iadr,sadr
 
4709
c
 
4710
      iadr(l)=l+l-1
 
4711
      sadr(l)=(l/2)+1
 
4712
c
 
4713
      if (lhs .ne. 1) then
 
4714
         call error(41)
 
4715
         return
 
4716
      endif
 
4717
      if(rhs .ne. 1) then
 
4718
         call error(42)
 
4719
         return
 
4720
      endif
 
4721
 
 
4722
      il=iadr(lstk(top))
 
4723
      
 
4724
      if(abs(istk(il)).ne.1) then
 
4725
         call funnam(ids(1,pt+1),'acos',iadr(lstk(top)))
 
4726
         fun=-1
 
4727
         return
 
4728
      endif
 
4729
 
 
4730
      if(istk(il).lt.0) then
 
4731
c        argument is passed by reference
 
4732
         ilr=il
 
4733
         il=iadr(istk(il+1))
 
4734
         mn=istk(il+1)*istk(il+2)
 
4735
         it=istk(il+3)
 
4736
         l=sadr(il+4)
 
4737
         lr=sadr(ilr+4)
 
4738
         err=lr+mn*(it+1)-lstk(bot)
 
4739
         if(err.gt.0) then
 
4740
            call error(17)
 
4741
            return
 
4742
         endif
 
4743
         call icopy(4,istk(il),1,istk(ilr),1)
 
4744
         lstk(top+1)=lr+mn*(it+1)
 
4745
      else
 
4746
         mn=istk(il+1)*istk(il+2)
 
4747
         it=istk(il+3)
 
4748
         l=sadr(il+4)
 
4749
         lr=l
 
4750
         ilr=il
 
4751
      endif
 
4752
      
 
4753
      if(mn.eq.0) return
 
4754
 
 
4755
      if(it.eq.0) then
 
4756
         itr=0
 
4757
         do i=0,mn-1
 
4758
            if(abs(stk(l+i)).gt.1.0d+0) then
 
4759
               itr=1
 
4760
               goto 20
 
4761
            endif
 
4762
         enddo
 
4763
 
 
4764
 20      continue
 
4765
         if(itr.eq.0) then
 
4766
c           argument is a real positive matrix with entries in [-1,1]
 
4767
            do i=0,mn-1
 
4768
               stk(lr+i)=acos(stk(l+i))
 
4769
            enddo
 
4770
         else
 
4771
c           argument is a real matrix with some entries outside [-1,1]
 
4772
            err=lr+2*mn-lstk(bot)
 
4773
            if(err.gt.0) then
 
4774
               call error(17)
 
4775
               return
 
4776
            endif
 
4777
            lstk(top+1)=lr+2*mn
 
4778
            do i=0,mn-1
 
4779
               call wacos(stk(l+i),0.0d+0,stk(lr+i),stk(lr+mn+i))
 
4780
            enddo
 
4781
            istk(ilr+3)=itr
 
4782
         endif
 
4783
      else
 
4784
c        argument is a complex matrix
 
4785
         do i=0,mn-1
 
4786
            call wacos(stk(l+i),stk(l+mn+i),stk(lr+i),stk(lr+i+mn))
 
4787
         enddo
 
4788
      endif
 
4789
      return
 
4790
      end
 
4791
 
 
4792
 
 
4793
      subroutine intnbprop(id)
 
4794
 
 
4795
c     Interface for number_properties (an interface on dlamch) :
 
4796
c    
 
4797
c         number_properties("eps")    -> machine epsilon
 
4798
c         number_properties("radix")  -> base
 
4799
c         number_properties("digits") -> number of digits for the mantissa
 
4800
c         number_properties("minexp") -> emin
 
4801
c         number_properties("maxexp") -> emax
 
4802
c         number_properties("huge")   -> max positive float
 
4803
c         number_properties("tiny")   -> min positive normalised float 
 
4804
c         number_properties("denorm") -> (boolean) true if denormalised number are used
 
4805
c         number_properties("tiniest")-> min positive denormalised float 
 
4806
c         
 
4807
c
 
4808
      implicit none
 
4809
 
 
4810
      INCLUDE '../stack.h'
 
4811
 
 
4812
      integer id(nsiz)
 
4813
 
 
4814
c     EXTERNAL FUNCTIONS
 
4815
      double precision dlamch
 
4816
      external         dlamch
 
4817
 
 
4818
c     EXTERNAL API FUNCTIONS
 
4819
      logical  checkrhs, checklhs, getsmat, getrmat, cremat, crebmat
 
4820
      external checkrhs, checklhs, getsmat, getrmat, cremat, crebmat
 
4821
 
 
4822
c     LOCAL VAR
 
4823
      integer topk
 
4824
      integer n, m, idxmat, mt, nt, lstr, nlstr, lm, lr, lc, i
 
4825
      integer   lmax
 
4826
      parameter(lmax = 10)
 
4827
      character*(lmax) inputstring
 
4828
      character*17 fname
 
4829
      double precision tiniest, b
 
4830
 
 
4831
c     TEXT
 
4832
      fname = 'number_properties'
 
4833
      topk=top
 
4834
      rhs=max(0,rhs)
 
4835
 
 
4836
      if (.not.checkrhs(fname,1,1)) return
 
4837
      if (.not.checklhs(fname,1,1)) return
 
4838
 
 
4839
c     1/ get the string
 
4840
      if( .not. getsmat(fname,topk,top,mt,nt,1,1,lstr,nlstr)) return
 
4841
c     rmq : pas de verif qu'il s'agit bien d'une matrice (1,1) ...
 
4842
c     on recupere la chaine dans la variable inputstring
 
4843
      lm = min(nlstr,lmax)
 
4844
      call cvstr(lm,istk(lstr),inputstring,1)
 
4845
c     complete (eventualy) the string with some blanks
 
4846
      inputstring(lm+1:lmax) = '          '
 
4847
 
 
4848
c     2/ go on
 
4849
      if     (inputstring(1:9) .eq. 'eps      ') then
 
4850
         if (.not.cremat(fname,top,0,1,1,lr,lc)) return
 
4851
         stk(lr) = dlamch('e')
 
4852
      elseif (inputstring(1:9) .eq. 'huge     ') then
 
4853
         if (.not.cremat(fname,top,0,1,1,lr,lc)) return
 
4854
         stk(lr) = dlamch('o')
 
4855
      elseif (inputstring(1:9) .eq. 'tiny     ') then
 
4856
         if (.not.cremat(fname,top,0,1,1,lr,lc)) return
 
4857
         stk(lr) = dlamch('u')
 
4858
      elseif (inputstring(1:9) .eq. 'radix    ') then
 
4859
         if (.not.cremat(fname,top,0,1,1,lr,lc)) return
 
4860
         stk(lr) = dlamch('b')
 
4861
      elseif (inputstring(1:9) .eq. 'digits   ') then
 
4862
         if (.not.cremat(fname,top,0,1,1,lr,lc)) return
 
4863
         stk(lr) = dlamch('n')
 
4864
      elseif (inputstring(1:9) .eq. 'minexp   ') then
 
4865
         if (.not.cremat(fname,top,0,1,1,lr,lc)) return
 
4866
         stk(lr) = dlamch('m')
 
4867
      elseif (inputstring(1:9) .eq. 'maxexp   ') then
 
4868
         if (.not.cremat(fname,top,0,1,1,lr,lc)) return
 
4869
         stk(lr) = dlamch('l')
 
4870
      elseif (inputstring(1:9) .eq. 'denorm   ') then
 
4871
         if (.not.crebmat(fname,top,1,1,lr)) return
 
4872
         if (dlamch('u') / dlamch('b') .gt. 0.d0) then
 
4873
            istk(lr) = 1
 
4874
         else
 
4875
            istk(lr) = 0
 
4876
         endif
 
4877
      elseif (inputstring(1:9) .eq. 'tiniest  ') then
 
4878
         if (.not.cremat(fname,top,0,1,1,lr,lc)) return
 
4879
         b = dlamch('b')
 
4880
         tiniest = dlamch('u')
 
4881
         if ( tiniest/b .ne. 0.d0 ) then
 
4882
c     denormalised number are used
 
4883
            do i = 1, dlamch('n') - 1
 
4884
               tiniest = tiniest / b
 
4885
            enddo
 
4886
         endif
 
4887
         stk(lr) = tiniest
 
4888
      else
 
4889
         buf=fname//' : unknown property kind'
 
4890
         call error(999)
 
4891
         return
 
4892
      endif
 
4893
      
 
4894
      end
 
4895
 
 
4896
 
 
4897
      subroutine intnearfl(id)
 
4898
 
 
4899
c     Interface for nearfloat :
 
4900
c    
 
4901
c         nearfloat("succ",x) -> succ of x
 
4902
c         nearfloat("pred",x) -> pred of x
 
4903
c
 
4904
      implicit none
 
4905
 
 
4906
      INCLUDE '../stack.h'
 
4907
 
 
4908
      integer id(nsiz)
 
4909
 
 
4910
c     EXTERNAL FUNCTIONS
 
4911
      double precision nearfloat
 
4912
      external         nearfloat
 
4913
 
 
4914
c     EXTERNAL API FUNCTIONS
 
4915
      logical  checkrhs, checklhs, getsmat, getrmat, cremat
 
4916
      external checkrhs, checklhs, getsmat, getrmat, cremat
 
4917
 
 
4918
c     LOCAL VAR
 
4919
      integer topk
 
4920
      integer n, m, idxmat, mt, nt, lstr, nlstr, lm, lr, lc, i
 
4921
      integer   lmax
 
4922
      parameter(lmax = 4)
 
4923
      character*(lmax) inputstring
 
4924
      character*9 fname
 
4925
 
 
4926
c     TEXT
 
4927
      fname = 'nearfloat'
 
4928
      topk=top
 
4929
      rhs=max(0,rhs)
 
4930
 
 
4931
      if (.not.checkrhs(fname,2,2)) return
 
4932
      if (.not.checklhs(fname,1,1)) return
 
4933
 
 
4934
c     1/ get the adress of the matrix
 
4935
      if( .not. getrmat(fname, topk, top, m, n, idxmat) ) return
 
4936
      top = top - 1
 
4937
c     2/ get the string
 
4938
      if( .not. getsmat(fname,topk,top,mt,nt,1,1,lstr,nlstr)) return
 
4939
c     pas de verif qu'il s'agit bien d'une matrice (1,1) ...
 
4940
c     on recupere la chaine dans la variable inputstring
 
4941
      lm = min(nlstr,lmax)
 
4942
      call cvstr(lm,istk(lstr),inputstring,1)
 
4943
 
 
4944
c     3/ go on
 
4945
      if (inputstring .eq. 'succ') then
 
4946
         if (.not.cremat(fname,top,0,m,n,lr,lc)) return
 
4947
         do i=0, m*n-1
 
4948
            stk(lr+i) = nearfloat(stk(idxmat+i),1.d0)
 
4949
         enddo
 
4950
 
 
4951
      elseif (inputstring .eq. 'pred') then
 
4952
         if (.not.cremat(fname,top,0,m,n,lr,lc)) return
 
4953
         do i=0, m*n-1
 
4954
            stk(lr+i) = nearfloat(stk(idxmat+i),-1.d0)
 
4955
         enddo
 
4956
 
 
4957
      else
 
4958
         buf=fname//' : unknown string specifier (must be pred or succ)'
 
4959
         call error(999)
 
4960
         return
 
4961
      endif
 
4962
 
 
4963
      end
 
4964
 
 
4965
      subroutine intdsearch(id)
 
4966
*
 
4967
*     interface for dsearch (Bruno le 10/12/2001)
 
4968
*
 
4969
*       [ind , occ, info] = dsearch(X, val [, ch])
 
4970
*        
 
4971
*       X and val must be real vectors (says of length m for X and n for val ), 
 
4972
*       if ch is not present then ch = 'c'  (dsearch on "intervals")
 
4973
*       ch must be 'd' or 'c'
 
4974
*
 
4975
*       ind is a vector with the same format than X
 
4976
*       occ is a vector with the same format than val (but with n-1
 
4977
*           components in the case ch='c')
 
4978
*       info is a scalar
 
4979
*
 
4980
      implicit none
 
4981
 
 
4982
      INCLUDE '../stack.h'
 
4983
 
 
4984
      integer id(nsiz)
 
4985
 
 
4986
c     EXTERNAL SUBROUTINES
 
4987
      external  dsearchc, dsearchd
 
4988
 
 
4989
c     EXTERNAL API FUNCTIONS
 
4990
      logical  checkrhs, checklhs, getsmat, getrvect, cremat, getrmat
 
4991
      external checkrhs, checklhs, getsmat, getrvect, cremat, getrmat
 
4992
 
 
4993
c     LOCAL VAR
 
4994
      integer topk, topl
 
4995
      integer mX, nX, lX, mval, nval, lval, mch, nch, lch, nlch
 
4996
      integer lind, mocc, nocc, locc, linfo, lc, j
 
4997
      character*1 ch
 
4998
      character*9 fname
 
4999
 
 
5000
c     STATEMENT FUNC
 
5001
      integer l, iadr,sadr
 
5002
      iadr(l)=l+l-1
 
5003
      sadr(l)=(l/2)+1
 
5004
 
 
5005
c     TEXT
 
5006
      fname = 'dsearch'
 
5007
      topk=top
 
5008
      rhs=max(0,rhs)
 
5009
 
 
5010
      if (.not.checkrhs(fname,2,3)) return
 
5011
      if (.not.checklhs(fname,1,3)) return
 
5012
 
 
5013
*     get ch
 
5014
      if (rhs .eq. 3) then
 
5015
         if( .not. getsmat(fname,topk,top,mch,nch,1,1,lch,nlch)) return
 
5016
         top = top - 1
 
5017
         call cvstr(1,istk(lch),ch,1)
 
5018
      else
 
5019
         ch = 'c'
 
5020
      endif
 
5021
      if (ch.ne.'c' .and. ch.ne.'d') then
 
5022
         buf=fname//' : unknown char specifier (must be ''c'' or ''d'')'
 
5023
         call error(999)
 
5024
         return
 
5025
      endif
 
5026
 
 
5027
c     get val 
 
5028
      if( .not. getrvect(fname, topk, top, mval, nval, lval) ) return
 
5029
      if (ch.eq.'d') then
 
5030
         if (mval*nval.lt.1) then
 
5031
            buf=fname//' : argument 2 must not be an empty vector'
 
5032
            call error(999)
 
5033
            return
 
5034
         endif
 
5035
         mocc = mval
 
5036
         nocc = nval
 
5037
      else    ! case ch='c'
 
5038
         if (mval*nval.lt.2) then
 
5039
            buf=fname//' : in the interval case, argument 2 must be'
 
5040
     $               //' a vector with length > 1'
 
5041
            call error(999)
 
5042
            return
 
5043
         endif
 
5044
         if (mval .eq. 1) then 
 
5045
            mocc = 1
 
5046
            nocc = nval - 1
 
5047
         else
 
5048
            mocc = mval - 1
 
5049
            nocc = nval
 
5050
         endif
 
5051
      endif
 
5052
*     verif that val is in strict increasing order
 
5053
      do j = 1, mval*nval-1
 
5054
         if (.not. stk(lval+j-1) .lt. stk(lval+j)) then  ! cette forme permet de detecter les nans
 
5055
            buf=fname//' : the array val (arg 2) is not well ordered'
 
5056
            call error(999)
 
5057
            return
 
5058
         endif
 
5059
      enddo
 
5060
      top = top - 1
 
5061
      
 
5062
c     get X
 
5063
      if( .not. getrmat(fname, topk, top, mX, nX, lX) ) return
 
5064
 
 
5065
c     reserve space for ind
 
5066
      if (.not.cremat(fname, topk+1, 0, mX, nX, lind, lc)) return
 
5067
 
 
5068
c     reserve space for occ
 
5069
      if (.not.cremat(fname, topk+2, 0, mocc, nocc, locc, lc)) return
 
5070
 
 
5071
c     reserve space for info
 
5072
      if (.not.cremat(fname, topk+3, 0, 1, 1, linfo, lc)) return
 
5073
 
 
5074
      if (mX.eq.0.or.nX.eq.0) then
 
5075
         stk(linfo)=0
 
5076
         call dset(mocc*nocc,0.0D0,stk(locc),1)
 
5077
      else
 
5078
 
 
5079
c     go on for the computation
 
5080
         if ( ch .eq. 'c') then
 
5081
            call dsearchc(stk(lX), mX*nX, stk(lval), mval*nval-1,
 
5082
     $           stk(lind), stk(locc), stk(linfo))
 
5083
         else 
 
5084
            call dsearchd(stk(lX), mX*nX, stk(lval), mval*nval, stk(lind
 
5085
     $           ),stk(locc), stk(linfo))
 
5086
         endif
 
5087
 
 
5088
c     int2db ... (normalement ca doit passer avec -1 sans copie
 
5089
C     supplementaire)
 
5090
         call int2db(mX*nX,     istk(iadr(lind)), -1, stk(lind), -1) 
 
5091
         call int2db(mocc*nocc, istk(iadr(locc)), -1, stk(locc), -1) 
 
5092
         call int2db(1,     istk(iadr(linfo)),-1, stk(linfo),-1) 
 
5093
      endif
 
5094
*     copie en "haut" 
 
5095
      topl = topk - rhs
 
5096
      if(lhs .ge. 1) then
 
5097
         call copyobj(fname,topk+1,topl+1)
 
5098
      endif
 
5099
      if(lhs .ge. 2) then
 
5100
         call copyobj(fname,topk+2,topl+2)
 
5101
      endif
 
5102
      if(lhs .ge. 3) then
 
5103
         call copyobj(fname,topk+3,topl+3)
 
5104
      endif
 
5105
      top=topl+lhs
 
5106
      return 
 
5107
      end
 
5108
 
 
5109
      subroutine intisequal(id)
 
5110
 
 
5111
c     Interface for isequal:
 
5112
 
 
5113
      implicit none
 
5114
 
 
5115
      INCLUDE '../stack.h'
 
5116
 
 
5117
      integer id(nsiz)
 
5118
      integer typ,m,n,l,il,il1,ilk,k,topk,top1,srhs,k1
 
5119
 
 
5120
c     EXTERNAL API FUNCTIONS
 
5121
      logical  checkrhs, checklhs
 
5122
      external checkrhs, checklhs
 
5123
      character*7 fname
 
5124
      integer iadr,sadr
 
5125
      integer equal
 
5126
      data equal/50/
 
5127
c     TEXT
 
5128
      iadr(l)=l+l-1
 
5129
      sadr(l)=(l/2)+1
 
5130
c
 
5131
      fname = 'isequal'
 
5132
 
 
5133
      topk=top
 
5134
      top1=top-rhs+1
 
5135
      rhs=max(0,rhs)
 
5136
 
 
5137
      if (.not.checkrhs(fname,2,2000000)) return
 
5138
      if (.not.checklhs(fname,1,1)) return
 
5139
c first check the types
 
5140
      typ=abs(istk(iadr(lstk(top1))))
 
5141
      do 10 k=1,rhs-1
 
5142
        if (abs(istk(iadr(lstk(top1+k)))).ne.typ) goto 60
 
5143
 10   continue
 
5144
c
 
5145
      if (typ.ge.15.and.typ.le.17) then
 
5146
         call setfunnam(ids(1,pt+1),'%l_isequal',10)
 
5147
         fun=-1
 
5148
         return
 
5149
      endif
 
5150
 
 
5151
      if(typ.gt.14) then
 
5152
         call funnam(ids(1,pt+1),'isequal',iadr(lstk(top1)))
 
5153
         fun=-1
 
5154
         return
 
5155
      endif
 
5156
 
 
5157
 
 
5158
c first check the dimensions
 
5159
      il=iadr(lstk(top-rhs+1))
 
5160
      if(istk(il).lt.0) il=iadr(istk(il+1))
 
5161
      m=istk(il+1)
 
5162
      n=istk(il+2)
 
5163
      do 20 k=2,rhs
 
5164
         il=iadr(lstk(top-rhs+k))
 
5165
         if(istk(il).lt.0) il=iadr(istk(il+1))
 
5166
         if(m.ne.istk(il+1).or.n.ne.istk(il+2)) goto 60
 
5167
 20   continue
 
5168
 
 
5169
      srhs=rhs
 
5170
 
 
5171
 30   do 40 k=2,srhs
 
5172
         call createref(iadr(lstk(top1)),top1,lstk(top1+1)-lstk(top1))
 
5173
         topk=top1+k-1
 
5174
         call createref(iadr(lstk(topk)),topk,lstk(topk+1)-lstk(topk))
 
5175
         fin=equal
 
5176
         rhs=2
 
5177
         call allops()
 
5178
         if(err.gt.0.or.err1.gt.0) return
 
5179
         if(icall.ne.0) then
 
5180
c     should not happen
 
5181
            rhs=srhs
 
5182
            top=top1-1+rhs
 
5183
            call funnam(ids(1,pt+1),'isequal',iadr(lstk(top-rhs+1)))
 
5184
            fun=-1
 
5185
            return
 
5186
         endif
 
5187
c     
 
5188
         il=iadr(lstk(top))
 
5189
         do 35  k1=1,istk(il+1)*istk(il+2)
 
5190
            if(istk(il+2+k1).eq.0) goto 60
 
5191
 35      continue
 
5192
 
 
5193
         top=top-1
 
5194
 40   continue
 
5195
c variables are equal
 
5196
      top=top1
 
5197
      il=iadr(lstk(top))
 
5198
      istk(il)=4
 
5199
      istk(il+1)=1
 
5200
      istk(il+2)=1
 
5201
      istk(il+3)=1
 
5202
      lstk(top+1)=sadr(il+4)
 
5203
      return
 
5204
c variables are different
 
5205
 60   continue
 
5206
      top=top1
 
5207
      il=iadr(lstk(top))
 
5208
      istk(il)=4
 
5209
      istk(il+1)=1
 
5210
      istk(il+2)=1
 
5211
      istk(il+3)=0
 
5212
      lstk(top+1)=sadr(il+4)
 
5213
      return
 
5214
 
 
5215
 
 
5216
      end