home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / libcruft / odepack / solsy.f < prev    next >
Text File  |  1996-09-28  |  3KB  |  69 lines
  1.       SUBROUTINE SOLSY (WM, IWM, X, TEM)
  2. CLLL. OPTIMIZE
  3.       INTEGER IWM
  4.       INTEGER IOWND, IOWNS,
  5.      1   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER,
  6.      2   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
  7.       INTEGER I, MEBAND, ML, MU
  8.       DOUBLE PRECISION WM, X, TEM
  9.       DOUBLE PRECISION ROWNS, 
  10.      1   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
  11.       DOUBLE PRECISION DI, HL0, PHL0, R 
  12.       DIMENSION WM(1), IWM(1), X(1), TEM(1)
  13.       COMMON /LS0001/ ROWNS(209),
  14.      2   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND,
  15.      3   IOWND(14), IOWNS(6), 
  16.      4   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER,
  17.      5   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
  18. C-----------------------------------------------------------------------
  19. C THIS ROUTINE MANAGES THE SOLUTION OF THE LINEAR SYSTEM ARISING FROM 
  20. C A CHORD ITERATION.  IT IS CALLED IF MITER .NE. 0.
  21. C IF MITER IS 1 OR 2, IT CALLS DGESL TO ACCOMPLISH THIS.
  22. C IF MITER = 3 IT UPDATES THE COEFFICIENT H*EL0 IN THE DIAGONAL
  23. C MATRIX, AND THEN COMPUTES THE SOLUTION.
  24. C IF MITER IS 4 OR 5, IT CALLS DGBSL.
  25. C COMMUNICATION WITH SOLSY USES THE FOLLOWING VARIABLES..
  26. C WM    = REAL WORK SPACE CONTAINING THE INVERSE DIAGONAL MATRIX IF
  27. C         MITER = 3 AND THE LU DECOMPOSITION OF THE MATRIX OTHERWISE. 
  28. C         STORAGE OF MATRIX ELEMENTS STARTS AT WM(3).
  29. C         WM ALSO CONTAINS THE FOLLOWING MATRIX-RELATED DATA..
  30. C         WM(1) = SQRT(UROUND) (NOT USED HERE),
  31. C         WM(2) = HL0, THE PREVIOUS VALUE OF H*EL0, USED IF MITER = 3.
  32. C IWM   = INTEGER WORK SPACE CONTAINING PIVOT INFORMATION, STARTING AT
  33. C         IWM(21), IF MITER IS 1, 2, 4, OR 5.  IWM ALSO CONTAINS BAND 
  34. C         PARAMETERS ML = IWM(1) AND MU = IWM(2) IF MITER IS 4 OR 5.
  35. C X     = THE RIGHT-HAND SIDE VECTOR ON INPUT, AND THE SOLUTION VECTOR
  36. C         ON OUTPUT, OF LENGTH N.
  37. C TEM   = VECTOR OF WORK SPACE OF LENGTH N, NOT USED IN THIS VERSION. 
  38. C IERSL = OUTPUT FLAG (IN COMMON).  IERSL = 0 IF NO TROUBLE OCCURRED. 
  39. C         IERSL = 1 IF A SINGULAR MATRIX AROSE WITH MITER = 3.
  40. C THIS ROUTINE ALSO USES THE COMMON VARIABLES EL0, H, MITER, AND N.
  41. C-----------------------------------------------------------------------
  42.       IERSL = 0
  43.       GO TO (100, 100, 300, 400, 400), MITER
  44.  100  CALL DGESL (WM(3), N, N, IWM(21), X, 0)
  45.       RETURN
  46. C
  47.  300  PHL0 = WM(2)
  48.       HL0 = H*EL0
  49.       WM(2) = HL0
  50.       IF (HL0 .EQ. PHL0) GO TO 330
  51.       R = HL0/PHL0
  52.       DO 320 I = 1,N
  53.         DI = 1.0D0 - R*(1.0D0 - 1.0D0/WM(I+2))
  54.         IF (DABS(DI) .EQ. 0.0D0) GO TO 390
  55.  320    WM(I+2) = 1.0D0/DI
  56.  330  DO 340 I = 1,N
  57.  340    X(I) = WM(I+2)*X(I)
  58.       RETURN
  59.  390  IERSL = 1
  60.       RETURN
  61. C
  62.  400  ML = IWM(1)
  63.       MU = IWM(2)
  64.       MEBAND = 2*ML + MU + 1
  65.       CALL DGBSL (WM(3), MEBAND, N, ML, MU, IWM(21), X, 0)
  66.       RETURN
  67. C----------------------- END OF SUBROUTINE SOLSY -----------------------
  68.       END 
  69.