summaryrefslogtreecommitdiffstats
path: root/scilab/modules/polynomials/sci_gateway/fortran/sci_f_roots.f
blob: f74384d1925464c36e38d9e5a9a10d72bab4195b (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
c Copyright (C) ????-2008 - INRIA - Serge Steer
c Copyright (C) 2011 - DIGITEO - Michael Baudin
c
c This file must be used under the terms of the CeCILL.
c This source file is licensed as described in the file COPYING, which
c you should have received as part of this distribution.  The terms
c are also available at
c http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt

      subroutine introots(id)
      INCLUDE 'stack.h'
      integer iadr, sadr
      integer id(nsiz)
      double precision t,sr,si
      integer vol
      integer fail
      logical ref,eigen
      integer algo_eig
      data algo_eig/14/
      integer algo_fast
      data algo_fast/15/
c
      iadr(l)=l+l-1
      sadr(l)=(l/2)+1
c
      lw=lstk(top+1)
      if(lhs.ne.1) then
         call error(41)
         return
      endif
      if(rhs.ne.1.and.rhs.ne.2) then
         call error(42)
         return
      endif
c
c if algo=="e", then eigen = .true.
c if algo=="f", then eigen = .false.
c
      eigen=.true.
c Read the algo option
      if (rhs.eq.2) then
         ilopt=iadr(lstk(top))
         if(istk(ilopt).lt.0) ilopt=iadr(istk(ilopt+1))
         if(istk(ilopt).ne.10) then
            err=2
            call error(55)
            return
         endif
         if(istk(ilopt+1).ne.1.or.istk(ilopt+2).ne.1) then
            err=2
            call error(89)
            return
         endif
         if(istk(ilopt+5).ne.2) then
            err=2
            call error(116)
            return
         else
            if (istk(ilopt+6)==algo_eig) then
              eigen=.true.
            elseif (istk(ilopt+6)==algo_fast) then
              eigen=.false.
            else
              err=2
              call error(116)
              return
            endif
         endif
         rhs=rhs-1
         top=top-1
      endif

      il1=iadr(lstk(top+1-rhs))
      ilr=il1
      if(istk(il1).lt.0) il1=iadr(istk(il1+1))
      ref=il1.ne.ilr
      if(istk(il1).gt.2) then
         fun=-1
         call funnam(ids(1,pt+1),'roots',il1)
         return
      endif
c
      m1=istk(il1+1)
      n1=istk(il1+2)
      mn1=m1*n1
      if(m1*n1.eq.0) then
        return
      endif

      it1=istk(il1+3)
c If "fast" algo was chosen and polynomial is complex, 
c then produce an error.
      if ( .not.eigen  .and. it1 .eq. 1 ) then
        err=2
        call error(116)
        return
      endif

      if(istk(il1).eq.1) then
c     for Matlab compatibility root of the vector of coefficients
         vol=mn1
         l1=sadr(il1+4)

         n=mn1
         if(.not.ref) then
            call dtild(n*(it1+1),stk(l1),1)
            lc=l1
         else
            lc=lw
            lw=lc+n*(it1+1)
            err=lw-lstk(bot)
            if(err.gt.0) then
               call error(17)
               return
            endif
            call dcopy(n,stk(l1),1,stk(lc),-1)
            if(it1.eq.1)  call dcopy(n,stk(l1+n),1,stk(lc+n),-1)
         endif
      else
         if(mn1.ne.1) then
            err=1
            call error(43)
            return
         endif

         if(istk(il1).lt.2) goto 24

         id1=il1+8
         lc=sadr(id1+mn1+1)
         vol=istk(id1+mn1)-1
         call icopy(4,istk(il1+4),1,id,1)
         n=vol
      endif
      l1=sadr(ilr+4)
 21   n=n-1
      if(n.lt.0) then
        goto 24
      endif
      t=abs(stk(lc+n))
      if(it1.eq.1) then
        t=t+abs(stk(lc+n+vol))
      endif
      if(t.eq.0.0d+0) then
        goto 21
      endif

      if ( .not.eigen  .and.n.gt.100) then
c If "fast" algo was chosen and polynomial has degree greater than 100, 
c then produce an error.
        err=2
        call error(116)
        return
      endif

      if (.not.eigen) then
c
c     real polynomial: rpoly algorithm
c     this alg is much more speedy, but it may happens that it gives
C     erroneous results without messages : example
C     roots(%s^31-8*%s^30+9*%s^29+0.995) should have two real roots near
C     1.355 and 6.65 and the other ones inside a circle centered in 0
C     with radius 1
C

         lp=max(lw,l1+2*n)
         err=lp+n+1-lstk(bot)
         if(err.gt.0) then
            call error(17)
            return
         endif
         call unsfdcopy(n+1,stk(lc),1,stk(lp),1)
         call dtild(n+1,stk(lp),1)
         call rpoly(stk(lp),n,stk(l1),stk(l1+n),fail)
         if(fail.eq.1) then
            call error(24)
            return
         elseif(fail.eq.2) then
            call error(74)
            return
         elseif(fail.eq.3) then
            call error(75)
            return

         endif
         istk(ilr)=1
         istk(ilr+1)=n
         istk(ilr+2)=1
         if(n.eq.0) istk(ilr+2)=0
         istk(ilr+3)=1
         lstk(top+1)=l1+2*n
         return
      else
c
c     Companion matrix method
        lw=lw+n*n*(it1+1)
        err=lw+n*(it1+1)-lstk(bot)
        if(err.gt.0) then
           call error(17)
           return
        endif
        sr=stk(lc+n)
        call unsfdcopy(n,stk(lc),-1,stk(lw),1)
        if(it1.eq.0) then
           call dscal(n,-1.0d+0/sr,stk(lw),1)
        else
           si=stk(lc+vol+n)
           t=sr*sr+si*si
           sr=-sr/t
           si=si/t
           call unsfdcopy(n,stk(lc+vol),-1,stk(lw+n),1)
           call wscal(n,sr,si,stk(lw),stk(lw+n),1)
        endif
        call dset(n*n*(it1+1),0.0d+0,stk(l1),1)
        call dset(n-1,1.0d+0,stk(l1+n),n+1)
        call unsfdcopy(n,stk(lw),1,stk(l1),1)
        if(it1.eq.1) then
          call unsfdcopy(n,stk(lw+n),1,stk(l1+n*n),1)
        endif
        lstk(top+1)=l1+n*n*(it1+1)
        istk(ilr)=1
        istk(ilr+1)=n
        istk(ilr+2)=n
        istk(ilr+3)=it1
        fin=3
        fun=2
c     *call* matds(r c)
        return
      endif

c     polynome de degre 0
 24   istk(ilr)=1
      istk(ilr+1)=0
      istk(ilr+2)=0
      istk(ilr+3)=0
      lstk(top+1)=sadr(ilr+4)
	  
      return
      end
c			=======================================