fix import/export of ierode common on Windows
[scilab.git] / scilab / modules / differential_equations / src / fortran / ainvg.f
1 C/MEMBR ADD NAME=AINVG,SSI=0
2       subroutine ainvg (res, adda, neq, t, y, ydot, miter,
3      1                   ml, mu, pw, ipvt, ier )
4 clll. optimize
5       
6       include 'stack.h'
7       external res, adda
8       integer neq, miter, ml, mu, ipvt, ier
9       integer i, lenpw, mlp1, nrowpw
10       double precision t, y, ydot, pw
11       dimension y(*), ydot(*), pw(*), ipvt(*)
12 c-----------------------------------------------------------------------
13 c%purpose
14 c this subroutine computes the initial value
15 c of the vector ydot satisfying
16 c     a * ydot = g(t,y)
17 c when a is nonsingular.  it is called by lsodi for
18 c initialization only, when istate = 0 .
19 c ainvg returns an error flag ier..
20 c   ier  =  0  means ainvg was successful.
21 c   ier .ge. 2 means res returned an error flag ires = ier.
22 c   ier .lt. 0 means the a-matrix was found to be singular.
23 c!
24 c-----------------------------------------------------------------------
25 c
26       if (miter .ge. 4)  go to 100
27 c
28 c full matrix case -----------------------------------------------------
29 c
30       lenpw = neq*neq
31       do 10  i = 1, lenpw
32    10    pw(i) = 0.0d+0
33 c
34       ier = 1
35       call res ( neq, t, y, pw, ydot, ier )
36       if(ierror.gt.0) return
37       if (ier .gt. 1) return
38 c
39       call adda ( neq, t, y, 0, 0, pw, neq )
40       if(ierror.gt.0) return
41       call dgefa ( pw, neq, neq, ipvt, ier )
42       if (ier .eq. 0) go to 20
43          ier = -ier
44          return
45    20 call dgesl ( pw, neq, neq, ipvt, ydot, 0 )
46       return
47 c
48 c band matrix case -----------------------------------------------------
49 c
50   100 continue
51       nrowpw = 2*ml + mu + 1
52       lenpw = neq * nrowpw
53       do 110  i = 1, lenpw
54   110    pw(i) = 0.0d+0
55 c
56       ier = 1
57       call res ( neq, t, y, pw, ydot, ier )
58       if(ierror.gt.0) return
59       if (ier .gt. 1) return
60 c
61       mlp1 = ml + 1
62       call adda ( neq, t, y, ml, mu, pw(mlp1), nrowpw )
63       if(ierror.gt.0) return
64       call dgbfa ( pw, nrowpw, neq, ml, mu, ipvt, ier )
65       if (ier .eq. 0) go to 120
66          ier = -ier
67          return
68   120 call dgbsl ( pw, nrowpw, neq, ml, mu, ipvt, ydot, 0 )
69       return
70 c-------------------- end of subroutine ainvg --------------------------
71       end