SUBROUTINE SBACC (G,LDG,NB,IR,MT,JT,JTPREV, IERR2) c Copyright (c) 1996 California Institute of Technology, Pasadena, CA. c ALL RIGHTS RESERVED. c Based on Government Sponsored Research NAS7-03001. C>> 1994-10-20 SBACC Krogh Changes to use M77CON C>> 1992-06-16 SBACC CLL C>> 1989-10-20 CLL C>> 1987-11-24 SBACC Lawson Initial code. c Sequential accumulation of data for a banded least-squares c problem. For solution use _BSOL. c ------------------------------------------------------------------ c Subroutine arguments c c G(,) [inout] Array in which user provides data and this subr c builds a packed representation of the orthogonally c transformed equivalent problem. c LDG [in] Leading dimensioning parameter for G(,). c NB [in] Bandwidth of data matrix. Set on first call for a c new problem and not changed thereafter. c IR [inout] User must set IR = 1 on 1st call for a new problem. c User must not thereafter alter IR. This subr will update c IR as IR := JT + min(NB+1, MT+max(IR-JT,0)). c IR identifies the row of G(,) in which user's new block c of data begins. c MT [in] Number of new rows of data being provided on the c current entry. c JT [in] Column index in the user's problem of the data being c provided in column 1 of G(,). c JTPREV [inout] Should not be set or altered by the user. c This subr will update JTPREV as JTPREV := JT. c JTPREV identifies a row in G(,) at which the mode of c packed storage changes. Data in row JTPREV and below c are subject to possible change as future data are c introduced, whereas data above row JTPREV are not. c IERR2 [inout] Error status indicator. Must be set to zero on c initial call to this subr for a new problem. Values on c return have meaning as follows: c 0 means no errors detected. c 1 means MT exceeds LDG - IR + 1 c 2 means JT .lt. JTPREV c 3 means JT .gt. min(JTPREV + NB, IR) which causes the c problem to be structurally singular. The transformations c by this subr can be completed in this case but the c resulting transformed problem cannot be solved simply c by use of _BSOL. c 4 means the condition of Error 3 above applys but computa- c tion must be abandonded due to storage limitations c indicated by MT > LDG - JT + 1. c 5 This subr will STOP because the input value of IERR2 was c nonzero. c ------------------------------------------------------------------ C C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUN 12 C Reference: 'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL, 1974 c Comments in this code refer to Algorithm Step numbers on c pp 213-217 of the book. c ------------------------------------------------------------------ c 1984 July 11, C. L. Lawson, JPL. Adapted code from book c to Fortran 77 for JPL MATH 77 library. c Prefixing subprogram names with S or D for s.p. or d.p. versions. c Using generic names for any intrinsic functions. c July 1987. Using _HTCC in place of _H12. C 1989-10-20 CLL Moved integer declaration earlier to avoid warning c msg from Cray compiler. C 1992-06-16 CLL Using "min" fcn in arg list of call to _HTCC c to avoid "index out of range" when using a bounds checker, even c though no reference would be made to the "out of range" location. c ------------------------------------------------------------------ c--S replaces "?": ?BACC, ?HTCC c Both versions use IERM1, IERV1 c ------------------------------------------------------------------ integer I, IE, IEROLD, IERR2, IG, IG1, IG2, IR, ISTEP, J, JG, JT integer JTPREV, K, KH, KH1, L, LDG, LP1, MH, MT, NB, NBP1 real G(LDG,NB+1), ZERO, UPARAM parameter(ZERO = 0.0E0) c ------------------------------------------------------------------ C C ALG. STEPS 1-4 ARE PERFORMED EXTERNAL TO THIS SUBROUTINE. C if (MT .le. 0) return if (IR .eq. 1) then c Initialize for new problem. JTPREV = 1 IERR2 = 0 elseif(IERR2 .ne. 0) then IEROLD = IERR2 IERR2 = 5 call IERM1('SBACC',IERR2,0, * 'Value of IERR2 nonzero on entry','IERR2 on entry',IEROLD,'.') stop endif c NBP1=NB+1 if( MT .gt. LDG - IR + 1) then IERR2 = 1 call IERM1('SBACC',IERR2,0, * 'MT must not exceed LDG - IR + 1','MT',MT,',') call IERV1('LDG',LDG,',') call IERV1(' IR', IR,'.') return endif C ALG. STEP 5 IF (JT.EQ.JTPREV) GO TO 70 C ALG. STEP 6 if(JT .lt. JTPREV) then IERR2 = 2 call IERM1('SBACC',IERR2,0, * 'Require JT .ge. JTPREV', 'JT',JT,',') call IERV1('JTPREV', JTPREV,'.') return endif C ALG. STEP 7 if (JT .gt. min(JTPREV+NB,IR) ) then IERR2 = 3 call IERM1('SBACC',IERR2,0, * 'With JT .gt. min(JTPREV+NB,IR) problem will be singular.', * 'JT',JT,',') call IERV1('JTPREV', JTPREV,',') call IERV1('NB', NB,',') call IERV1('IR', IR,'.') c c Note: At this point we know the problem will be structurally c singular and thus will not be solvable simply by use of _BSOL. c This does not preclude, however, completion of the transforma- c tion to the band triangular form and so we continue. The c calling program will be informed of this condition by the set- c ting of of IERR2 = 3, and thus can quit or continue as desired. c If this is not due to a usage mistake c then the user should probably introduce more data, modify his c or her mathematical model, or use a stabilizing method such as c is discussed on pp. 218-219 of the L & H book. c if(MT .gt. LDG - JT + 1) then IERR2 = 4 call IERM1('SBACC',IERR2,0, * 'MT must not exceed LDG - JT + 1','MT',MT,',') call IERV1('LDG', LDG,',') call IERV1('JT', JT,'.') return endif C ALG. STEPS 8-9 do 10 I=1,MT IG1=JT+MT-I IG2=IR+MT-I do 10 J=1,NBP1 10 G(IG1,J)=G(IG2,J) C ALG. STEP 10 IE=JT-IR do 20 I=1,IE IG=IR+I-1 do 20 J=1,NBP1 20 G(IG,J)=ZERO C ALG. STEP 11 IR=JT endif c c Here we must shift some rows left by variable amounts. c min(NB-1,IR-JTPREV-1) is the number of rows to be shifted. c K = min(L,JT-JTPREV) is the extent of the left shift c for row JTPREV+L. c C ALG. STEP 12-13 do 50 L=1, min(NB-1,IR-JTPREV-1) C ALG. STEP 14 K=min(L,JT-JTPREV) C ALG. STEP 15 LP1=L+1 IG=JTPREV+L do 40 I=LP1,NB 40 G(IG,I-K)=G(IG,I) C ALG. STEP 16 do 50 JG = NBP1-K, NB 50 G(IG,JG)=ZERO C ALG. STEP 17 JTPREV=JT C ALG. STEPS 18-19 70 continue c c MH = No. of rows in the block to which the orthogonal c triangularization will be applied. c KH1 = Number of pivot columns to be used in the orghogonal c transformations. c MH=IR+MT-JTPREV KH1 = min(NBP1,MH-1) ISTEP = IR-JTPREV+1 C ALG. STEP 20 do 80 I=1,KH1 80 CALL SHTCC (1,I,max(I+1,ISTEP),MH,G(JTPREV,I),UPARAM, * G(JTPREV,min(I+1,NBP1)),LDG,NBP1-I) C ALG. STEP 21 KH=min(NBP1,MH) IR=JTPREV+KH C ALG. STEP 22 if (KH .GE. NBP1) then C ALG. STEP 23 do 90 I=1,NB 90 G(IR-1,I)=ZERO C ALG. STEP 24 endif C ALG. STEP 25 return end