bf009f1bc4851731622689aac0093fc8b2322b60
[scilab.git] / scilab / modules / polynomials / sci_gateway / fortran / sci_f_psimp.f
1 c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 c Copyright (C) ????-2008 - INRIA
3 c
4 c This file must be used under the terms of the CeCILL.
5 c This source file is licensed as described in the file COPYING, which
6 c you should have received as part of this distribution.  The terms
7 c are also available at
8 c http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
9       subroutine intpsimp(id)
10       INCLUDE 'stack.h'
11       integer iadr, sadr
12       integer id(nsiz)
13       integer vola,volb
14       logical chkvar
15 c
16       iadr(l)=l+l-1
17       sadr(l)=(l/2)+1
18 c
19
20 c     simp(num,den)
21 c
22       if(lhs.ne.2) then
23          call error(41)
24          return
25       endif
26 c
27       ila=iadr(lstk(top+1-rhs))
28
29       il=ila
30       if(istk(il).lt.0) il=iadr(istk(il+1))
31       if(istk(il).gt.2) then
32          fun=-1
33          call funnam(ids(1,pt+1),'simp',il)
34          return
35       endif
36 c
37       ilb=iadr(lstk(top))
38
39       il=ilb
40       if(istk(il).lt.0) il=iadr(istk(il+1))
41       if(istk(il).gt.2) then
42          fun=-1
43          call funnam(ids(1,pt+1),'simp',il)
44          return
45       endif
46 c
47       if(istk(ila).lt.0) then
48 c     move b and copy value of a in place of its reference
49          k=istk(ila+2)
50          vola=lstk(k+1)-lstk(k)
51          volb=lstk(top+1)-lstk(top)
52          err=lstk(top)+vola+volb-lstk(bot)
53          if(err.gt.0) then
54             call error(17)
55             return
56          endif
57          call unsfdcopy(volb,stk(lstk(top)),-1,stk(lstk(top-1)+vola),-1)
58          call unsfdcopy(vola,stk(lstk(k)),1,stk(lstk(top-1)),1)
59          lstk(top)=lstk(top-1)+vola
60          ilb=iadr(lstk(top))
61          lstk(top+1)=lstk(top)+volb
62       endif
63
64       if(istk(ilb).lt.0) then
65 c      copy value of b in place of its reference
66          k=istk(ilb+2)
67          volb=lstk(k+1)-lstk(k)
68          err=lstk(top)+volb-lstk(bot)
69          if(err.gt.0) then
70             call error(17)
71             return
72          endif
73          call unsfdcopy(volb,stk(lstk(k)),1,stk(lstk(top)),1)
74          lstk(top+1)=lstk(top)+volb
75       endif
76       lw=lstk(top+1)
77
78 c
79       if(istk(ila+3).ne.0.or.istk(ilb+3).ne.0) then
80          fun=-1
81          call funnam(ids(1,pt+1),'simp',ilb)
82          return
83       endif
84 c
85       mna=istk(ila+1)*istk(ila+2)
86       id(1)=0
87       if(istk(ila).eq.2) then
88          ida=ila+8
89          la=sadr(ida+mna+1)
90          call icopy(4,istk(ila+4),1,id,1)
91       else
92          la=sadr(ila+4)
93          ida=iadr(lw)
94          lw=sadr(ida+mna+1)
95          err=lw-lstk(bot)
96          if(err.gt.0) then
97             call error(17)
98             return
99          endif
100          do 74 i=1,mna+1
101             istk(ida+i-1)=i
102  74      continue
103       endif
104 c
105       mnb=istk(ilb+1)*istk(ilb+2)
106       if(istk(ilb).eq.2) then
107          idb=ilb+8
108          lb=sadr(idb+mnb+1)
109          if(id(1).eq.0) then
110             call icopy(4,istk(ilb+4),1,id,1)
111          else
112             if(.not.chkvar(id,istk(ilb+4))) then
113                call error(43)
114                return
115             endif
116          endif
117       else
118          lb=sadr(ilb+4)
119          idb=iadr(lw)
120          lw=sadr(idb+mna+1)
121          err=lw-lstk(bot)
122          if(err.gt.0) then
123             call error(17)
124             return
125          endif
126          do 75 i=1,mna+1
127             istk(idb+i-1)=i
128  75      continue
129       endif
130
131       if(mnb.ne.mna)then
132          call error(60)
133          return
134       endif
135 c
136 c     determine max of the degrees
137       na=0
138       nb=0
139       do 76 i=1,mna
140          na=max(na,istk(ida+i)-istk(ida-1+i))
141          nb=max(nb,istk(idb+i)-istk(idb-1+i))
142  76   continue
143 c     preserve adress of the beginning of a and b coefficients
144
145       lar=la
146       numx0 = stk(la)
147       numx1 = stk(la+1)
148       numx2 = stk(la+2)
149       numx3 = stk(la+3)
150
151       lbr=lb
152       denx0 = stk(lb)
153       denx1 = stk(lb+1)
154       denx2 = stk(lb+2)
155       denx3 = stk(lb+3)
156 c     allocate memory for intermediate results
157       law=lw
158       lbw=law+na+1
159       lw=lbw+nb+1
160
161 c     simplify
162       la1=la
163       lb1=lb
164       do 79 i=1,mna
165          na=istk(ida+i)-istk(ida-1+i)-1
166          nb=istk(idb+i)-istk(idb-1+i)-1
167          ierr=lstk(bot)-lw
168          call  dpsimp(stk(la),na,stk(lb),nb,stk(law),nnum,
169      $        stk(lbw),nden,stk(lw),ierr)
170          if(ierr.eq.1) then
171             call error(27)
172             return
173          elseif(ierr.eq.2) then
174             call msgs(43,i)
175          endif
176 c     .  copy overwrite initial polynomials with simplified ones
177          call dcopy(nnum,stk(law),1,stk(la1),1)
178          call dcopy(nden,stk(lbw),1,stk(lb1),1)
179
180          numx0 = stk(la1)
181          numx1 = stk(la1+1)
182          numx2 = stk(la1+2)
183          numx3 = stk(la1+3)
184
185          denx0 = stk(lb1)
186          denx1 = stk(lb1+1)
187          denx2 = stk(lb1+2)
188          denx3 = stk(lb1+3)
189
190          la=la+na+1
191          lb=lb+nb+1
192          la1=la1+nnum
193          lb1=lb1+nden
194          istk(ida-1+i)=nnum
195          istk(idb-1+i)=nden
196  79   continue
197 c
198 c     form vector of pointers from vector of degrees+1
199       ma=1
200       mb=1
201       do 80 i=1,mna+1
202          na=istk(ida-1+i)
203          nb=istk(idb-1+i)
204          istk(ida-1+i)=ma
205          istk(idb-1+i)=mb
206          ma=ma+na
207          mb=mb+nb
208  80   continue
209 c
210 c     compute position of the a and b simplified in the result
211       lstk(top)=lar+istk(ida+mna)-1
212       il=iadr(lstk(top))
213 c
214 c     put new b variable in place
215       if(istk(ilb).eq.2) then
216 c     b matrice de polynome
217          l=sadr(il+9+mna)
218 c     .  move b data up
219          call icopy(9+mna,istk(ilb),1,istk(il),1)
220          call unsfdcopy(istk(il+8+mna),stk(lbr),1,stk(l),1)
221          l=l+istk(il+8+mna)-1
222       else
223 c     b matrice de scalaires
224          call icopy(4,istk(ilb),1,istk(il),1)
225          l=sadr(il+4)
226          call unsfdcopy(mna,stk(lbr),1,stk(l),1)
227          l=l+mna
228       endif
229       lstk(top+1)=l
230       end
231
232 c                       =======================================