2 c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 c Copyright (C) Scilab Enterprises - 2013 - Paul Bignier
5 c This file must be used under the terms of the CeCILL.
6 c This source file is licensed as described in the file COPYING, which
7 c you should have received as part of this distribution. The terms
8 c are also available at
9 c http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
11 subroutine daskri(fname)
12 c ====================================================================
14 c ====================================================================
16 c This code was inspired by sci_f_dasrt.f from the same folder.
17 c daskr() has more optional parameters than dasrt(), so the gateway
18 c has to adapt to the new options
22 character*(nlgh+1) namjac
24 integer iadr,sadr,gettype
26 double precision atol,rtol,t0
27 integer info(20),topk,topw
28 integer LENWP, LENIWP, maxl, kmp, nrmax, mxnit, mxnj, mxnh, lsoff
29 logical hotstart,type,getexternal,getrvect
30 logical checkrhs,checklhs,getrmat,cremat,getscalar
31 double precision tout,tstop,maxstep,stepin
32 double precision epli, stptol, epinit
33 character*(nlgh+1) namer,namej,names,namep,namepj
34 common /dassln/ namer,namej,names,namep,namepj
35 external bresd,bjacd,bsurfd,bpsold,bpjacd
36 external setfresd,setfjacd,setfsurfd,setfpsold,setfpjacd
38 data atol/1.d-7/,rtol/1.d-9/
43 c SCILAB function : daskr
44 c --------------------------
45 c [y0,nvs[,hotdata]]=daskr(y0,t0,t1[,atol[,rtol]],res[,jac],nh,h[,info
46 c [,psol][,pjac]][,hotdata])
54 if (.not.checkrhs(fname,6,13)) return
55 if (.not.checklhs(fname,2,3)) return
56 c Checking variable y0 (number 1)
57 c -------------------------------
59 if(.not.getrmat(fname,topk,ky,n1,m1,l1))return
64 if (.not.cremat(fname,topw,0,n1,1,lydot,lc)) return
67 call dset(n1,0.0d0,stk(lydot),1)
73 il1 = iadr(lstk(top-rhs+1))
76 c Checking variable t0 (number 2)
77 c ----------------------------
79 if(.not.getscalar(fname,topk,kt0,lr2))return
81 c Checking variable t1 (number 3)
82 c -------------------------------
83 if(.not.getrmat(fname,topk,top-rhs+3,m3,n3,l3))return
86 c Checking variable atol (number 4)
87 c --------------------------------
89 itype = gettype(top-rhs+4)
90 if ( itype .ne. 1) then
91 if (.not.cremat(fname,topw,0,1,1,latol,lc)) return
93 if (.not.cremat(fname,topw,0,1,1,lrtol,lc)) return
101 if(.not.getrvect(fname,topk,top-rhs+4,m4,n4,latol))return
103 c Checking variable rtol (number 5)
104 c --------------------------------
105 itype = gettype(top-rhs+5)
106 if (itype .ne. 1) then
107 if (.not.cremat(fname,topw,0,1,1,lrtol,lc)) return
114 if(.not.getrvect(fname,topk,top-rhs+5,m5,n5,lrtol))return
116 if(m5.ne.m4.or.(m5.ne.1.and.m5.ne.neq)) then
126 c Checking variable res (number 6)
127 c --------------------------------
128 1105 kres=top-rhs+6-iskip
129 if (.not.getexternal(fname,topk,kres,namer,type,
132 c Checking variable number 7
133 c -----------------------------
141 c . info or jac ? get first element type to decide
143 if (istk(il7).lt.0) il7=istk(il7+1)
146 if (abs(istk(iadr(l71))).eq.11.or.
147 $ abs(istk(iadr(l71))).eq.13) then
152 if((is7.ne.10).and.(is7.ne.11).and.(is7.ne.13)) then
157 if (.not.getexternal(fname,topk,kjac,namej,type,
162 c Checking variable number 8
163 c -----------------------------
164 if(.not.getscalar(fname,topk,top-rhs+8-iskip,lr8))return
166 c Checking variable number 9
167 ksurf=top-rhs+9-iskip
168 if (.not.getexternal(fname,topk,ksurf,names,type,
171 c Checking variable info (number 10)
172 c ------------------------------------
173 kinfo = top-rhs+10-iskip
174 if (kinfo.gt.top) then
192 il10 = iadr(lstk(top-rhs+10-iskip))
193 if (istk(il10) .ne. 15) then
194 c Default info values
215 c -- subvariable tstop(info) --
216 il10e1=iadr(l10+istk(il10+1+1)-1)
217 l10e1 = sadr(il10e1+4)
218 m10e1 = istk(il10e1+1)*istk(il10e1+2)
227 c -- subvariable imode(info) --
228 il10e2=iadr(l10+istk(il10+1+2)-1)
229 l10e2 = sadr(il10e2+4)
233 c -- subvariable band(info) --
234 il10e3=iadr(l10+istk(il10+1+3)-1)
235 m10e3 =istk(il10e3+2)*istk(il10e3+2)
236 l10e3 = sadr(il10e3+4)
239 elseif(m10e3.eq.2) then
250 c -- subvariable maxstep(info) --
251 il10e4=iadr(l10+istk(il10+1+4)-1)
252 m10e4 =istk(il10e4+2)*istk(il10e4+2)
253 l10e4 = sadr(il10e4+4)
262 c -- subvariable stepin(info) --
263 il10e5=iadr(l10+istk(il10+1+5)-1)
264 m10e5 =istk(il10e5+2)*istk(il10e5+2)
265 l10e5 = sadr(il10e5+4)
274 c -- subvariable nonneg(info) --
275 il10e6=iadr(l10+istk(il10+1+6)-1)
276 l10e6 = sadr(il10e6+4)
280 c -- subvariable consistent(info) --
281 il10e7=iadr(l10+istk(il10+1+7)-1)
282 m10e7 =istk(il10e7+2)*istk(il10e7+2)
283 l10e7 = sadr(il10e7+4)
284 if(m10e7.eq.0.or.(m10e7.eq.1.and.stk(l10e7).eq.0)) then
285 c info is then [] or [0]
288 c info then looks like list(..., [+-1 +-1 ... +-1 +-1],...)
290 if (info(10).eq.0.or.info(10).eq.2) then
295 c Make no further changes here, but copy the [+-1] tab
296 c in 'iwork' once its size 'liw' is defined.
300 c -- subvariable iteration(info) --
301 il10e8=iadr(l10+istk(il10+1+8)-1)
302 l10e8 = sadr(il10e8+4)
311 c -- subvariable defaultKrylov(info) --
312 il10e9=iadr(l10+istk(il10+1+9)-1)
313 m10e9 =istk(il10e9+2)*istk(il10e9+2)
314 l10e9 = sadr(il10e9+4)
317 c maxl and kmp need default values
321 c info then looks like list(..., [maxl kmp nrmax epli],...)
330 c -- subvariable justConsistentComp(info) --
331 il10e10=iadr(l10+istk(il10+1+10)-1)
332 l10e10 = sadr(il10e10+4)
334 if(proceed.eq.0) then
337 c Check that info(11) = 1, meaning that the provided initial values
339 if (info(11).eq.1) info(14) = 1
343 c -- subvariable psolJac(info) --
344 il10e11=iadr(l10+istk(il10+1+11)-1)
345 l10e11 = sadr(il10e11+4)
347 if (psolJac.eq.0) then
354 c -- subvariable excludeAlgebraic(info) --
355 il10e12=iadr(l10+istk(il10+1+12)-1)
356 m10e12 =istk(il10e12+2)*istk(il10e12+2)
357 l10e12 = sadr(il10e12+4)
358 if (m10e12.eq.0) then
361 c info then looks like list(..., [+-1 +-1 ... +-1 +-1],...)
363 if (info(10).eq.0.or.info(10).eq.2) then
368 c Make no further changes here, but copy the [+-1] tab
369 c in 'iwork' once its size 'liw' is defined.
373 c -- subvariable defaultHeuristic(info) --
374 il10e13=iadr(l10+istk(il10+1+13)-1)
375 m10e13 =istk(il10e13+2)*istk(il10e13+2)
376 l10e13 = sadr(il10e13+4)
377 if (m10e13.eq.0) then
380 c info then looks like list(..., [mxnit mxnj lsoff stptol epinit],...)
386 stptol = stk(l10e9+4)
387 epinit = stk(l10e9+5)
391 c -- subvariable verbosity(info) --
392 il10e14=iadr(l10+istk(il10+1+14)-1)
393 l10e14 = sadr(il10e14+4)
394 verbosity=stk(l10e14)
395 if (verbosity.eq.1) then
397 elseif(verbosity.eq.2) then
403 c Checking variable psol (number 11)
404 c --------------------------------------
405 8 kpsol=top-rhs+11-iskip
406 if (kpsol.gt.top) then
410 il11 = iadr(lstk(top-rhs+11-iskip))
412 if((pType.ne.10).and.(pType.ne.11).and.(pType.ne.13)) then
416 if (.not.getexternal(fname,topk,kpsol,namep,type,
419 c Checking variable pjac (number 12)
420 c --------------------------------------
421 9 kpjac=top-rhs+12-iskip
422 if (kpjac.gt.top) then
426 il12 = iadr(lstk(top-rhs+12-iskip))
428 if((pType.ne.10).and.(pType.ne.11).and.(pType.ne.13)) then
432 if (.not.getexternal(fname,topk,kpjac,namepj,type,
436 if(rhs.eq.13-iskip) then
439 c Checking variable hotdata (number 13)
440 c --------------------------------------
441 il13 = iadr(lstk(top-rhs+13-iskip))
442 if (istk(il13) .ne. 1) then
447 n13 = istk(il13+1)*istk(il13+2)
449 elseif(rhs.ne.12-iskip) then
454 c --------------------Work Tables
455 if (.not.cremat(fname,topw,0,1,2,lw15,lc)) return
457 if (.not.cremat(fname,topw,0,1,30,lw17,lc)) return
460 c Set ipar and rpar to 0 by default
462 istk(iadr(lw15)+i) = 0
468 if (.not.cremat(fname,topw,0,1,nh,lgr,lc)) return
472 c base = lrw = 60 then augment size
473 c according to the case (full dense, banded, ...)
475 if (info(12).eq.0) then
476 lrw = base + max(maxord+4,7)*neq + 3*nh
477 if (info(6).eq.0) then
478 c For the full (dense) JACOBIAN case
480 elseif(info(5).eq.1) then
481 c For the banded user-defined JACOBIAN case
482 lrw = lrw + (2*ml+mu+1)*neq
483 elseif(info(5).eq.0) then
484 c For the banded finite-difference-generated JACOBIAN case
485 lrw = lrw + (2*ml+mu+1)*neq + 2*(neq/(ml+mu+1)+1)
487 elseif(info(12).eq.1) then
488 c LENWP is the length ot the rwork segment containing
489 c the matrix elements of the preconditioner P
491 lrw = base + (maxord+5)*neq + 3*nh
492 $ + (maxl + 3 + min(1,maxl-kmp))*neq
493 $ + (maxl+3)*maxl + 1 + LENWP
495 if(info(16).eq.1) lrw = lrw + neq
498 c base = liw = 40, then augment size according to the case
500 if (info(12).eq.0) then
502 elseif(info(12).eq.1) then
503 c LENIWP is the length ot the iwork segment containing
504 c the matrix indexes of the preconditioner P
505 c (compressed sparse row format)
509 if(info(10).eq.1.or.info(10).eq.3) liw = liw + neq
510 if(info(11).eq.1.or.info(16).eq.1) liw = liw + neq
511 if(.not.hotstart) then
512 if (.not.cremat(fname,topw,0,1,lrw,lrwork,lc)) return
514 if (.not.cremat(fname,topw,0,1,sadr(liw)+1,liwork,lc)) return
517 if(lrw+liw.gt.n13) then
518 c Test if output (/input) array 'hotdata' is large enough
519 c to hold at least rwork and iwork
526 call entier(liw,stk(liwork),istk(iadr(liwork)))
529 if(info(4).eq.1) then
532 if(info(7).eq.1) then
533 stk(lrwork+1)=maxstep
535 if(info(8).eq.1) then
538 if(info(6).eq.1) then
539 istk(iadr(liwork))=ml
540 istk(iadr(liwork)+1)=mu
542 if(info(11).eq.1) then
544 istk(iadr(liwork)+LID+i) = stk(l10e7+i)
547 if(info(12).eq.1) then
548 istk(iadr(liwork)+26) = LENWP
549 istk(iadr(liwork)+27) = LENIWP
551 if(info(13).eq.1) then
552 istk(iadr(liwork)+23) = maxl
553 istk(iadr(liwork)+24) = kmp
554 istk(iadr(liwork)+25) = nrmax
555 stk(lrwork+10) = epli
557 if (info(15).eq.1) then
559 istk(iadr(lw15)) = neq
560 istk(iadr(lw15)+1) = neq
561 istk(iadr(lw15)+2) = 2
562 istk(iadr(lw15)+3) = 2
563 istk(iadr(lw15)+4) = 1
567 if(info(16).eq.1) then
569 istk(iadr(liwork)+LID+i) = stk(l10e12+i)
572 if(info(17).eq.1) then
573 istk(iadr(liwork)+31) = mxnit
574 istk(iadr(liwork)+32) = mxnj
575 istk(iadr(liwork)+33) = mxnh
576 istk(iadr(liwork)+34) = lsoff
577 stk(lrwork+13) = stptol
578 stk(lrwork+14) = epinit
580 istk(iadr(liwork)+16) = lrw
581 istk(iadr(liwork)+17) = liw
583 c create header of wp
584 if (.not.cremat(fname,topw,0,neq*neq,1,lr,lc)) return
587 c create header of wp and iwp
588 if (.not.cremat(fname,topw,0,neq*neq,2,lr,lc)) return
592 c Structure d'info pour les externals
597 istk(ilext+1)=ilext+6
598 istk(ilext+2)=ilext+10
599 istk(ilext+3)=ilext+14
600 istk(ilext+4)=ilext+17
601 istk(ilext+5)=ilext+21
616 c header used to create wp
617 istk(ilext+18)=iheaderwp
618 c header used to create iwp
619 istk(ilext+19)=iheaderiwp
620 c header of y used to create b
625 c header of t0 used to create neq, t0, h, cj
627 c header of y used to create y, ydot, rewt and savr
640 if(hotstart) info(1)=1
649 margin=(k-1)*(2*n1+1)+4
651 if(lhs.eq.3) lw1=lw1+4+lrw+liw
652 if(lw1-lstk(bot).gt.0) then
654 call msgstxt('Not enough memory to go further')
658 if (tout .eq. t0) then
660 call unsfdcopy(n1,stk(l1),1,stk(lyri+1),1)
661 call unsfdcopy(n1,stk(lydot),1,stk(lyri+n1+1),1)
668 call unsfdcopy(n1,stk(l1),1,stk(lyri+1),1)
669 call unsfdcopy(n1,stk(lydot),1,stk(lyri+n1+1),1)
672 if (info(15).eq.1) then
673 call ddaskr(bresd, n1, t0, stk(l1), stk(lydot),
674 & stk(lyri), info, stk(lrtol), stk(latol), idid,
675 & stk(lrwork), lrw, istk(iadr(liwork)), liw, stk(lw15),
676 & istk(il17), bpjacd, bpsold, bsurfd, nh, istk(lgroot))
678 call ddaskr(bresd, n1, t0, stk(l1), stk(lydot),
679 & stk(lyri), info, stk(lrtol), stk(latol), idid,
680 & stk(lrwork), lrw, istk(iadr(liwork)), liw, stk(lw15),
681 & istk(il17), bjacd, bpsold, bsurfd, nh, istk(lgroot))
683 C SUBROUTINE DDASKR(RES, NEQ, T, Y, YPRIME, TOUT, INFO, RTOL, ATOL,
684 C * IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, JAC, PSOL,
687 if(err.gt.0.or.err1.gt.0) return
689 C A step was successfully taken in the intermediate-output mode.
690 C The code has not yet reached TOUT.
695 elseif(idid.eq.2.or.idid.eq.4) then
696 C The integration to TSTOP was successfully completed (T=TSTOP)
699 elseif(idid.eq.3) then
700 C The integration to TOUT was successfully completed (T=TOUT) by
701 C stepping past TOUT. Y and Ydot are obtained by interpolation.
705 elseif(idid.eq.5) then
706 C One or more root found
710 elseif(idid.eq.-1) then
711 C A large amount of work has been expended (About 500 steps)
712 call msgstxt('Too many steps necessary to reach next '//
713 & 'required time discretization point.')
714 call msgstxt('Change discretization of time vector t '//
715 & 'or decrease accuracy.')
718 elseif(idid.eq.-2) then
719 C The error tolerances are too stringent.
723 c buf='The error tolerances are too stringent'
726 elseif(idid.eq.-3) then
727 C The local error test cannot be satisfied because you specified
728 C a zero component in ATOL and the corresponding computed solution
729 C component is zero. Thus, a pure relative error test is impossible
730 C for this component.
731 buf='atol and computed test value are zero.'
734 elseif(idid.eq.-5) then
735 C There were repeated failures in the evaluation
736 C or processing of the preconditioner (in JAC).
737 buf='Cannot evaluate the preconditioner.'
740 elseif(idid.eq.-6) then
741 C Repeated error test failures on the last attempted step.
742 call msgstxt('A singularity in the solution '//
745 elseif(idid.eq.-7) then
746 C The corrector could not converge.
747 call msgstxt('May be inaccurate or ill-conditioned '//
750 elseif(idid.eq.-8) then
751 C The matrix of partial derivatives is singular.
752 buf='The matrix of partial derivatives is singular'//
753 & 'Some of your equations may be redundant.'
756 elseif(idid.eq.-9) then
757 C The corrector could not converge. there were repeated error
758 c test failures in this step.
759 call msgstxt('Either ill-posed problem or '//
760 & 'discontinuity or singularity encountered.')
762 elseif(idid.eq.-10) then
763 call msgstxt('External ''res'' returned many times'//
766 elseif(idid.eq.-11) then
767 C IRES equal to -2 was encountered and control is being returned to the
769 buf='Error in external ''res''.'
772 elseif(idid.eq.-12) then
773 C DDASKR failed to compute the initial Yprime.
774 buf='daskr failed to compute the initial Ydot.'
777 elseif(idid.eq.-13) then
778 C An unrecoverable error was encountered inside the user's PSOL routine,
779 C and control is being returned to the calling program.
780 buf='Error in external ''Psol''.'
783 elseif(idid.eq.-14) then
784 C The Krylov linear system solver could not achieve convergence.
785 buf='The Krylov linear system solver did not converge.'
788 elseif(idid.eq.-33) then
789 C The code has encountered trouble from which
790 C it cannot recover. A message is printed
791 C explaining the trouble and control is returned
792 C to the calling program. For example, this occurs
793 C when invalid input is detected.
794 call msgstxt('daskr encountered trouble.')
804 c Variable de sortie: y0
807 if(k.eq.0) istk(ilyr+1)=0
812 c Variable de sortie: roots
816 err=lw+4+nh+1-lstk(bot)
829 if(istk(lgroot+i).ne.0) then
832 istk(ilw+2)=istk(ilw+2)+1
838 if(lhs.eq.2) goto 1150
840 c Variable de sortie: rwork
844 err=lw+4+lrw+liw-lstk(bot)
854 call unsfdcopy(lrw,stk(lrwork),1,stk(lw),1)
855 call int2db(liw,istk(iadr(liwork)),1,stk(lw+lrw),1)
859 c Remise en place de la pile
860 1150 call unsfdcopy(lw-lw0,stk(lw0),1,stk(l0),1)