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