scilab/modules/differential_equations/sci_gateway/fortran/bpjacd.f

   1 c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
   2 c Copyright (C) Scilab Enterprises - 2013
   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 c
  10       subroutine bpjacd (res, ires, neq, t, y, ydot, rewt, savr,
  11      *                     wk, h, cj, wp, iwp, ier, rpar, ipar)
  12 c
  13 c ======================================================================
  14 c     Management of external dealing with preconditioning of linear systems.
  15 c ======================================================================
  16 c
  17       INCLUDE 'stack.h'
  18       integer iadr,sadr
  19 c
  20       logical allowptr
  21 c
  22       character tmpbuf * (bsiz)
  23       double precision res(*), t, y(*), ydot(*), rewt(*), savr(*),
  24      *                  wk(*), h, cj, wp(*), rpar(*)
  25       integer ires, neq, iwp(*), ier, ipar(*)
  26       double precision dneq
  27       integer vol,tops,nordre,hsize
  28       data nordre/5/,mlhs/3/
  29 c
  30       iadr(l)=l+l-1
  31       sadr(l)=(l/2)+1
  32 c
  33 c
  34
  35       if (ddt .eq. 4) then
  36          write(tmpbuf(1:12),'(3i4)') top,r,sym
  37          call basout(io,wte,' bpjacd  top:'//tmpbuf(1:4))
  38       endif
  39 c     nordre is the external's order number in the data structure,
  40 c     mlhs (mrhs) is the number of output (input) parameters
  41 c     of the simulator
  42 c
  43       mrhs=8
  44       ierror=0
  45 c
  46       ilp=iadr(lstk(top))
  47       il=istk(ilp+nordre)
  48 c
  49       tops=istk(il)
  50       ils=iadr(lstk(tops))
  51 c
  52       if(istk(ils).eq.10) then
  53 c     Fortran simulator case
  54          call fpjacd (res, ires, neq, t, y, ydot, rewt, savr,
  55      *                wk, h, cj, wp, iwp, ier, rpar, ipar)
  56          return
  57       endif
  58 c     external is a Scilab function
  59
  60 c     On return ierror=1 is used to notify to the ode solver that
  61 c     scilab was not able to evaluate the external
  62       ierror=1
  63
  64 c     Putting Fortran arguments on Scilab stack
  65 c
  66 c     Minimum entry arguments for the simulator. The value of these arguments
  67 c     comes from the Fortran context (call list)
  68 c+
  69 c     function res
  70       dneq = neq
  71       call ftob(dneq,1,istk(il+1))
  72       if(err.gt.0.or.err1.gt.0) return
  73       call ftob(t,1,istk(il+1))
  74       if(err.gt.0.or.err1.gt.0) return
  75       call ftob(y,neq,istk(il+2))
  76       if(err.gt.0.or.err1.gt.0) return
  77       call ftob(ydot,neq,istk(il+2))
  78       if(err.gt.0.or.err1.gt.0) return
  79       call ftob(h,1,istk(il+1))
  80       if(err.gt.0.or.err1.gt.0) return
  81       call ftob(cj,1,istk(il+1))
  82       if(err.gt.0.or.err1.gt.0) return
  83       call ftob(rewt,neq,istk(il+2))
  84       if(err.gt.0.or.err1.gt.0) return
  85       call ftob(savr,neq,istk(il+2))
  86       if(err.gt.0.or.err1.gt.0) return
  87 c+
  88 c
  89       if(istk(ils).eq.15) goto 10
  90 c
  91 c     Retrieving the simulator's address
  92       fin=lstk(tops)
  93 c
  94       goto 40
  95 c     If the simulator is defined by a list
  96  10   nelt=istk(ils+1)
  97       l=sadr(ils+3+nelt)
  98       ils=ils+2
  99 c
 100 c     Retrieving the simulator's address
 101       fin=l
 102 c
 103 c     Managing the additional simulator parameters coming from
 104 c     the context (elements of the list describing the simulator)
 105 c
 106       nelt=nelt-1
 107       if(nelt.ne.0) then
 108          l=l+istk(ils+1)-istk(ils)
 109          vol=istk(ils+nelt+1)-istk(ils+1)
 110          if(top+1+nelt.ge.bot) then
 111             call error(18)
 112             return
 113          endif
 114          err=lstk(top+1)+vol-lstk(bot)
 115          if(err.gt.0) then
 116             call error(17)
 117             return
 118          endif
 119          call unsfdcopy(vol,stk(l),1,stk(lstk(top+1)),1)
 120          do 11 i=1,nelt
 121             top=top+1
 122             lstk(top+1)=lstk(top)+istk(ils+i+1)-istk(ils+i)
 123  11      continue
 124          mrhs=mrhs+nelt
 125       endif
 126  40   continue
 127 c
 128 c     Executing the macro defining the simulator
 129 c
 130       pt=pt+1
 131       if(pt.gt.psiz) then
 132          call  error(26)
 133          return
 134       endif
 135
 136       ids(1,pt)=lhs
 137       ids(2,pt)=rhs
 138       rstk(pt)=1001
 139       lhs=mlhs
 140       rhs=mrhs
 141       niv=niv+1
 142       fun=0
 143 c
 144       icall=5
 145       include 'callinter.h.f'
 146 c
 147  200  lhs=ids(1,pt)
 148       rhs=ids(2,pt)
 149       pt=pt-1
 150       niv=niv-1
 151 c+
 152 c     Transferring the output to Fortran
 153       call btof(ier,1)
 154       if(err.gt.0.or.err1.gt.0) return
 155 c      call btof(iwp,2*neq*neq)
 156       il2=iadr(lstk(top))
 157       hsize=4
 158       n=istk(il2+1)*istk(il2+2)*(istk(il2+3)+1)
 159       liwp = 2*neq*neq
 160 c     Test if the variable on the stack has same type and size as the theoretical iwp
 161       if (istk(il2).ne.1.or.n.ne.liwp) then
 162          call error(98)
 163          return
 164       endif
 165       lx=sadr(il2+hsize)
 166       do 900 i=1,liwp
 167          iwp(i) = stk(lx+i-1)
 168 900   continue
 169       top = top-1
 170       if(err.gt.0.or.err1.gt.0) return
 171       lwp = neq*neq
 172       call btof(wp,lwp)
 173       if(err.gt.0.or.err1.gt.0) return
 174 c+
 175
 176 c     Normal return ierror set to 0
 177       ierror=0
 178       return
 179 c
 180  9999 continue
 181       niv=niv-1
 182       if(err1.gt.0) then
 183          lhs=ids(1,pt)
 184          rhs=ids(2,pt)
 185          pt=pt-1
 186          fun=0
 187       endif
 188       return
 189       end
 190