1      SUBROUTINE DBSOL (MODE,G,LDG,NB,IR,JTPREV,X,N,RNORM, IERR3)
2c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
3c ALL RIGHTS RESERVED.
4c Based on Government Sponsored Research NAS7-03001.
5C>> 1994-11-11 DBSOL  Krogh  Declared all vars.
6C>> 1994-10-20 DBSOL  Krogh  Changes to use M77CON
7C>> 1987-11-24 DBSOL  Lawson Initial code.
8c--D replaces "?": ?BSOL
9c     Solution of banded least squares problem following use of _BACC to
10c     accumulate the data.
11c     ------------------------------------------------------------------
12c                            Subroutine arguments
13c
14C     MODE  [in]  Flag selecting operation to be done.
15c          = 1     SOLVE R*X=Y   WHERE R AND Y ARE IN THE G( ) ARRAY
16C                  AND X WILL BE STORED IN THE X( ) ARRAY.
17C            2     SOLVE (R**T)*X=Y   WHERE R IS IN G( ),
18C                  Y IS INITIALLY IN X( ), AND X REPLACES Y IN X( ),
19C            3     SOLVE R*X=Y   WHERE R IS IN G( ).
20C                  Y IS INITIALLY IN X( ), AND X REPLACES Y IN X( ).
21C
22c     G(,)  [in]  Array in which the transformed problem data has
23c           been left by _BACC.
24c     LDG   [in]  Leading dimensioning parameter for G(,).
25c     NB    [in]  Bandwidth of data matrix.
26c     IR    [in]  Must have value set by previous call to _BACC.
27c     JTPREV    [in]  Must have value set by previous call to _BACC.
28c     X()   [inout]  Array of length at least N.
29c           On entry with MODE = 2 or 3 contains the y vector.
30c           In all (non-error) cases contains the solution vector on
31c           return.
32c     N     [in]  Specifies the dimension of the desired solution
33c           vector.  Can be set smaller than the max value the data
34c           would support.
35c     RNORM  [out]  Set to the euclidean norm of the residual vector
36c           when MODE = 1.  Set to zero otherwise.
37c     IERR3  [out]  =0 means no errors detected.
38c                  =1 means a diagonal element is zero.
39c     ------------------------------------------------------------------
40C     C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUN 12
41C     Reference:   'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL, 1974
42c     Comments in this code refer to Algorithm Step numbers on
43c     pp 213-217 of the book.
44c     ------------------------------------------------------------------
45c     1984 July 11, C. L. Lawson, JPL.  Adapted code from book
46c     to Fortran 77 for JPL MATH 77 library.
47c     Prefixing subprogram names with S or D for s.p. or d.p. versions.
48c     Using generic names for any intrinsic functions.
49c     July 1987.  Using _HTCC in place of _H12.
50c     ------------------------------------------------------------------
51c     Subprograms called: IERM1, IERV1
52c     ------------------------------------------------------------------
53      integer MODE, LDG, NB, JTPREV, IR, N, IERR3
54      integer J, NP1, IRM1, II, I, L, IE, JG, IX, I1, I2
55      double precision G(LDG,NB+1),X(N), ZERO, RNORM, RSQ, S
56      parameter( ZERO = 0.0D0)
57c     ------------------------------------------------------------------
58C
59      IERR3 = 0
60      RNORM=ZERO
61      if( MODE .eq. 1) then
62C                                   ********************* MODE = 1
63C                                   ALG. STEP 26
64         DO 20 J=1,N
65   20       X(J)=G(J,NB+1)
66         RSQ=ZERO
67         NP1=N+1
68         IRM1=IR-1
69         if (NP1 .le. IRM1) then
70            DO 30 J=NP1,IRM1
71   30            RSQ=RSQ+G(J,NB+1)**2
72            RNORM=SQRT(RSQ)
73         endif
74      endif
75c
76      if( MODE .eq. 2) go to 90
77C                                   ********************* MODE = 1 or 3
78C                                   ALG. STEP 27
79           DO 80 II=1,N
80           I=N+1-II
81C                                   ALG. STEP 28
82           S=ZERO
83           L=MAX(0,I-JTPREV)
84C                                   ALG. STEP 29
85           if (I .ne. N) then
86C                                   ALG. STEP 30
87              IE=MIN(N+1-I,NB)
88              DO 60 J=2,IE
89                 JG=J+L
90                 IX=I-1+J
91   60            S=S+G(I,JG)*X(IX)
92           endif
93C                                   ALG. STEP 31
94           if (G(I,L+1) .eq. ZERO) go to 130
95   80      X(I)=(X(I)-S)/G(I,L+1)
96C                                   ALG. STEP 32
97      return
98C     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
99C                                   ********************* MODE = 2
100   90      DO 120 J=1,N
101              S=ZERO
102              if (J .ne. 1) then
103                 I1=MAX(1,J-NB+1)
104                 I2=J-1
105                 DO 100 I=I1,I2
106                    L=J-I+1+MAX(0,I-JTPREV)
107  100               S=S+X(I)*G(I,L)
108              endif
109              L=MAX(0,J-JTPREV)
110              if (G(J,L+1) .eq. ZERO) then
111                 I = J
112                 go to 130
113              endif
114  120         X(J)=(X(J)-S)/G(J,L+1)
115      return
116C     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117  130 continue
118      IERR3 = I
119      call IERM1('DBSOL', IERR3, 0,'Singular matrix',
120     * 'MODE', MODE, ',')
121      call IERV1('Index of zero diag elt', I, '.')
122      return
123      end
124