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