1      SUBROUTINE SACCUM (A,LDA,N,B,LDB,NB,IR1, NROWS, NCOUNT)
2c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
3c ALL RIGHTS RESERVED.
4c Based on Government Sponsored Research NAS7-03001.
5c>> 1996-03-30 SACCUM Krogh   Added external statement.
6C>> 1994-11-11 SACCUM Krogh   Declared all vars.
7C>> 1994-10-20 SACCUM Krogh  Changes to use M77CON
8C>> 1987-11-24 SACCUM Lawson  Initial code.
9c--S replaces "?": ?ACCUM, ?ROTG, ?HTCC, ?NRM2
10c     Sequential accumulation of rows of data for a linear least
11c     squares problem.  Using Givens orthogonal transformations when
12c     NROWS is small and Householder orthogonal transformations when
13c     NROWS is larger.
14c
15c        On first call user must set IR1 = 1.  On each return this subr
16c     will update IR1 to min(IR1+NROWS, N+2).  The user should not alter
17c     IR1 while processing data for the same case.
18c
19c        On all calls, including the first, the user sets NROWS to
20c     indicate the number of rows of new data being provided.  The user
21c     must put the new data in rows IR1 through IR1 + NROWS - 1 of the
22c     arrays A() and B().  On return, with IR1 updated to
23c     min(IR1+NROWS, N+2), the first IR1-1 rows of A() will contain
24c     transformed data on and to the right of the diagonal, and zeros
25c     to the left of the diagonal.  The first IR1-1 rows of B() will
26c     also contain transformed data.
27c     ------------------------------------------------------------------
28c     Reference: C. L. Lawson & R. J. Hanson,
29c                Solving Least Squares Problems, Prentice-Hall, 1974.
30c     Original code by R. J. Hanson, JPL, Sept 11, 1968.
31c     Adapted to Fortran 77 for the JPL MATH77 library
32c     by Lawson, 6/2/87.
33c     ------------------------------------------------------------------
34c                      Subroutine arguments
35c
36c     A(,)  [inout]
37c     LDA   [in]
38c     N     [in]
39c     B(,)  [inout]
40c     LDB   [in]
41c     NB    [in]
42c     IR1   [inout]
43c     NROWS  [in]
44c     NCOUNT  [out]
45c     ------------------------------------------------------------------
46c     Subprograms called directly: SHTCC, SNRM2, SROTG, IERM1, IERV1
47c     ------------------------------------------------------------------
48      external SNRM2
49      integer    IR1, NB, L1, LAST, LDA, LDB, N, NROWS, MARK
50      integer    M, IP, IROW, J, I, NCOUNT
51      real              SNRM2
52      real              A(LDA,N), B(LDB,NB), UPARAM, ZERO
53      real             CFAC, SFAC, TEMP
54      parameter(ZERO = 0.0E0, MARK = 2)
55c     ------------------------------------------------------------------
56C
57      if (NROWS .le. 0) return
58C
59      M=IR1+NROWS-1
60      if(M .gt. LDA .or. M .gt. LDB) then
61         call IERM1('SACCUM',1,0,
62     *   'Require LDA .ge. M and LDB .ge M where M = IR1+NROWS-1',
63     *   'IR1',IR1, ',')
64         call IERV1('NROWS',NROWS, ',')
65         call IERV1('LDA',LDA, ',')
66         call IERV1('LDB',LDB, '.')
67         return
68      endif
69      if(M .le. N) then
70         LAST = M-1
71      else
72         LAST = N
73      endif
74c
75      if(NROWS .le. MARK) then
76c      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
77c           Use Givens transformation for better efficiency when the
78c           number of rows being introduced is .le. MARK.
79c
80         do 50 IP = 1, LAST
81            L1 = max(IR1, IP+1)
82            do 30 IROW = L1, M
83               call SROTG(A(IP,IP), A(IROW,IP), CFAC, SFAC)
84               A(IROW,IP) = ZERO
85               do 10 J = IP+1,N
86                  TEMP = A(IP,J)
87                  A(IP,J)   =   CFAC * TEMP + SFAC * A(IROW,J)
88                  A(IROW,J) =  -SFAC * TEMP + CFAC * A(IROW,J)
89   10          continue
90               do 20 J = 1,NB
91                  TEMP = B(IP,J)
92                  B(IP,J)   =   CFAC * TEMP + SFAC * B(IROW,J)
93                  B(IROW,J) =  -SFAC * TEMP + CFAC * B(IROW,J)
94   20          continue
95   30       continue
96   50    continue
97c        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
98      else
99c        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100c                 Use Householder transformations
101c                 when the number of rows being introduced exceeds MARK.
102c
103         do 80 IP = 1, LAST
104            L1 = max(IR1, IP+1)
105            call SHTCC (1,IP,L1,M,A(1,IP),UPARAM,
106     *                           A(1,min(IP+1,N)),LDA,N-IP)
107            call SHTCC (2, IP, L1, M, A(1,IP),UPARAM,
108     *                           B(1,1),LDB,NB)
109C
110C                      Clear elements just made implicitly zero.
111C
112            do 60 I = L1, min(M,N+1)
113               A(I,IP)=ZERO
114   60       continue
115   80    continue
116c        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117      endif
118c
119      if(IR1 .eq. 1) then
120         NCOUNT = NROWS
121      else
122         NCOUNT = NCOUNT + NROWS
123      endif
124      IR1=min(IR1 + NROWS, N+2)
125      if (M .le. N+1) return
126C
127C     Pack lengths of B() column vectors below row N
128C     to single locations each.
129C
130      do 90 J=1,NB
131         B(N+1,J) = SNRM2 (M-N, B(N+1,J), 1)
132   90 continue
133      return
134      end
135