diff libcruft/odepack/solsy.f @ 1644:395bb6d6c096

[project @ 1995-12-14 08:32:12 by jwe] Initial revision
author jwe
date Thu, 14 Dec 1995 08:32:21 +0000
parents
children 44ed237bdc1e
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libcruft/odepack/solsy.f	Thu Dec 14 08:32:21 1995 +0000
@@ -0,0 +1,68 @@
+      SUBROUTINE SOLSY (WM, IWM, X, TEM)
+CLLL. OPTIMIZE
+      INTEGER IWM
+      INTEGER IOWND, IOWNS,
+     1   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER,
+     2   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
+      INTEGER I, MEBAND, ML, MU
+      DOUBLE PRECISION WM, X, TEM
+      DOUBLE PRECISION ROWNS, 
+     1   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
+      DOUBLE PRECISION DI, HL0, PHL0, R 
+      DIMENSION WM(1), IWM(1), X(1), TEM(1)
+      COMMON /LS0001/ ROWNS(209),
+     2   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND,
+     3   IOWND(14), IOWNS(6), 
+     4   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER,
+     5   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
+C-----------------------------------------------------------------------
+C THIS ROUTINE MANAGES THE SOLUTION OF THE LINEAR SYSTEM ARISING FROM 
+C A CHORD ITERATION.  IT IS CALLED IF MITER .NE. 0.
+C IF MITER IS 1 OR 2, IT CALLS DGESL TO ACCOMPLISH THIS.
+C IF MITER = 3 IT UPDATES THE COEFFICIENT H*EL0 IN THE DIAGONAL
+C MATRIX, AND THEN COMPUTES THE SOLUTION.
+C IF MITER IS 4 OR 5, IT CALLS DGBSL.
+C COMMUNICATION WITH SOLSY USES THE FOLLOWING VARIABLES..
+C WM    = REAL WORK SPACE CONTAINING THE INVERSE DIAGONAL MATRIX IF
+C         MITER = 3 AND THE LU DECOMPOSITION OF THE MATRIX OTHERWISE. 
+C         STORAGE OF MATRIX ELEMENTS STARTS AT WM(3).
+C         WM ALSO CONTAINS THE FOLLOWING MATRIX-RELATED DATA..
+C         WM(1) = SQRT(UROUND) (NOT USED HERE),
+C         WM(2) = HL0, THE PREVIOUS VALUE OF H*EL0, USED IF MITER = 3.
+C IWM   = INTEGER WORK SPACE CONTAINING PIVOT INFORMATION, STARTING AT
+C         IWM(21), IF MITER IS 1, 2, 4, OR 5.  IWM ALSO CONTAINS BAND 
+C         PARAMETERS ML = IWM(1) AND MU = IWM(2) IF MITER IS 4 OR 5.
+C X     = THE RIGHT-HAND SIDE VECTOR ON INPUT, AND THE SOLUTION VECTOR
+C         ON OUTPUT, OF LENGTH N.
+C TEM   = VECTOR OF WORK SPACE OF LENGTH N, NOT USED IN THIS VERSION. 
+C IERSL = OUTPUT FLAG (IN COMMON).  IERSL = 0 IF NO TROUBLE OCCURRED. 
+C         IERSL = 1 IF A SINGULAR MATRIX AROSE WITH MITER = 3.
+C THIS ROUTINE ALSO USES THE COMMON VARIABLES EL0, H, MITER, AND N.
+C-----------------------------------------------------------------------
+      IERSL = 0
+      GO TO (100, 100, 300, 400, 400), MITER
+ 100  CALL DGESL (WM(3), N, N, IWM(21), X, 0)
+      RETURN
+C
+ 300  PHL0 = WM(2)
+      HL0 = H*EL0
+      WM(2) = HL0
+      IF (HL0 .EQ. PHL0) GO TO 330
+      R = HL0/PHL0
+      DO 320 I = 1,N
+        DI = 1.0D0 - R*(1.0D0 - 1.0D0/WM(I+2))
+        IF (DABS(DI) .EQ. 0.0D0) GO TO 390
+ 320    WM(I+2) = 1.0D0/DI
+ 330  DO 340 I = 1,N
+ 340    X(I) = WM(I+2)*X(I)
+      RETURN
+ 390  IERSL = 1
+      RETURN
+C
+ 400  ML = IWM(1)
+      MU = IWM(2)
+      MEBAND = 2*ML + MU + 1
+      CALL DGBSL (WM(3), MEBAND, N, ML, MU, IWM(21), X, 0)
+      RETURN
+C----------------------- END OF SUBROUTINE SOLSY -----------------------
+      END