GNU Octave: libcruft/odepack/solsy.f Source File

Go to the documentation of this file.
00001       SUBROUTINE SOLSY (WM, IWM, X, TEM)
00002 CLLL. OPTIMIZE
00003       INTEGER IWM
00004       INTEGER IOWND, IOWNS,
00005      1   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER,
00006      2   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
00007       INTEGER I, MEBAND, ML, MU
00008       DOUBLE PRECISION WM, X, TEM
00009       DOUBLE PRECISION ROWNS, 
00010      1   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
00011       DOUBLE PRECISION DI, HL0, PHL0, R 
00012       DIMENSION WM(*), IWM(*), X(*), TEM(*)
00013       COMMON /LS0001/ ROWNS(209),
00014      2   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND,
00015      3   IOWND(14), IOWNS(6), 
00016      4   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER,
00017      5   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
00018 C-----------------------------------------------------------------------
00019 C THIS ROUTINE MANAGES THE SOLUTION OF THE LINEAR SYSTEM ARISING FROM 
00020 C A CHORD ITERATION.  IT IS CALLED IF MITER .NE. 0.
00021 C IF MITER IS 1 OR 2, IT CALLS DGETRS TO ACCOMPLISH THIS.
00022 C IF MITER = 3 IT UPDATES THE COEFFICIENT H*EL0 IN THE DIAGONAL
00023 C MATRIX, AND THEN COMPUTES THE SOLUTION.
00024 C IF MITER IS 4 OR 5, IT CALLS DGBTRS.
00025 C COMMUNICATION WITH SOLSY USES THE FOLLOWING VARIABLES..
00026 C WM    = REAL WORK SPACE CONTAINING THE INVERSE DIAGONAL MATRIX IF
00027 C         MITER = 3 AND THE LU DECOMPOSITION OF THE MATRIX OTHERWISE. 
00028 C         STORAGE OF MATRIX ELEMENTS STARTS AT WM(3).
00029 C         WM ALSO CONTAINS THE FOLLOWING MATRIX-RELATED DATA..
00030 C         WM(1) = SQRT(UROUND) (NOT USED HERE),
00031 C         WM(2) = HL0, THE PREVIOUS VALUE OF H*EL0, USED IF MITER = 3.
00032 C IWM   = INTEGER WORK SPACE CONTAINING PIVOT INFORMATION, STARTING AT
00033 C         IWM(21), IF MITER IS 1, 2, 4, OR 5.  IWM ALSO CONTAINS BAND 
00034 C         PARAMETERS ML = IWM(1) AND MU = IWM(2) IF MITER IS 4 OR 5.
00035 C X     = THE RIGHT-HAND SIDE VECTOR ON INPUT, AND THE SOLUTION VECTOR
00036 C         ON OUTPUT, OF LENGTH N.
00037 C TEM   = VECTOR OF WORK SPACE OF LENGTH N, NOT USED IN THIS VERSION. 
00038 C IERSL = OUTPUT FLAG (IN COMMON).  IERSL = 0 IF NO TROUBLE OCCURRED. 
00039 C         IERSL = 1 IF A SINGULAR MATRIX AROSE WITH MITER = 3.
00040 C THIS ROUTINE ALSO USES THE COMMON VARIABLES EL0, H, MITER, AND N.
00041 C-----------------------------------------------------------------------
00042       IERSL = 0
00043       GO TO (100, 100, 300, 400, 400), MITER
00044  100  CALL DGETRS ( 'N', N, 1, WM(3), N, IWM(21), X, N, INLPCK)
00045       RETURN
00046 C
00047  300  PHL0 = WM(2)
00048       HL0 = H*EL0
00049       WM(2) = HL0
00050       IF (HL0 .EQ. PHL0) GO TO 330
00051       R = HL0/PHL0
00052       DO 320 I = 1,N
00053         DI = 1.0D0 - R*(1.0D0 - 1.0D0/WM(I+2))
00054         IF (DABS(DI) .EQ. 0.0D0) GO TO 390
00055  320    WM(I+2) = 1.0D0/DI
00056  330  DO 340 I = 1,N
00057  340    X(I) = WM(I+2)*X(I)
00058       RETURN
00059  390  IERSL = 1
00060       RETURN
00061 C
00062  400  ML = IWM(1)
00063       MU = IWM(2)
00064       MEBAND = 2*ML + MU + 1
00065       CALL DGBTRS ( 'N', N, ML, MU, 1, WM(3), MEBAND, IWM(21), X, N, 
00066      * INLPCK)
00067       RETURN
00068 C----------------------- END OF SUBROUTINE SOLSY -----------------------
00069       END