1 SUBROUTINE SQRSL(X,LDX,N,K,QRAUX,Y,QY,QTY,B,RSD,XB,JOB,INFO) 2C***BEGIN PROLOGUE SQRSL 3C***DATE WRITTEN 780814 (YYMMDD) 4C***REVISION DATE 820801 (YYMMDD) 5C***REVISION HISTORY (YYMMDD) 6C 000330 Modified array declarations. (JEC) 7C***CATEGORY NO. D9,D2A1 8C***KEYWORDS LINEAR ALGEBRA,LINPACK,MATRIX,ORTHOGONAL TRIANGULAR,SOLVE 9C***AUTHOR STEWART, G. W., (U. OF MARYLAND) 10C***PURPOSE Applies the output of SQRDC to compute coordinate trans- 11C formations projections, and least squares solutions. 12C***DESCRIPTION 13C 14C SQRSL applies the output of SQRDC to compute coordinate 15C transformations, projections, and least squares solutions. 16C For K .LE. MIN(N,P), let XK be the matrix 17C 18C XK = (X(JPVT(1)),X(JPVT(2)), ... ,X(JPVT(K))) 19C 20C formed from columnns JPVT(1), ... ,JPVT(K) of the original 21C N x P matrix X that was input to SQRDC (if no pivoting was 22C done, XK consists of the first K columns of X in their 23C original order). SQRDC produces a factored orthogonal matrix Q 24C and an upper triangular matrix R such that 25C 26C XK = Q * (R) 27C (0) 28C 29C This information is contained in coded form in the arrays 30C X and QRAUX. 31C 32C On Entry 33C 34C X REAL(LDX,P) 35C X contains the output of SQRDC. 36C 37C LDX INTEGER 38C LDX is the leading dimension of the array X. 39C 40C N INTEGER 41C N is the number of rows of the matrix XK. It must 42C have the same value as N in SQRDC. 43C 44C K INTEGER 45C K is the number of columns of the matrix XK. K 46C must not be greater than MIN(N,P), where P is the 47C same as in the calling sequence to SQRDC. 48C 49C QRAUX REAL(P) 50C QRAUX contains the auxiliary output from SQRDC. 51C 52C Y REAL(N) 53C Y contains an N-vector that is to be manipulated 54C by SQRSL. 55C 56C JOB INTEGER 57C JOB specifies what is to be computed. JOB has 58C the decimal expansion ABCDE, with the following 59C meaning. 60C 61C If A .NE. 0, compute QY. 62C If B,C,D, or E .NE. 0, compute QTY. 63C If C .NE. 0, compute B. 64C If D .NE. 0, compute RSD. 65C If E .NE. 0, compute XB. 66C 67C Note that a request to compute B, RSD, or XB 68C automatically triggers the computation of QTY, for 69C which an array must be provided in the calling 70C sequence. 71C 72C On Return 73C 74C QY REAL(N). 75C QY contains Q*Y, if its computation has been 76C requested. 77C 78C QTY REAL(N). 79C QTY contains TRANS(Q)*Y, if its computation has 80C been requested. Here TRANS(Q) is the 81C transpose of the matrix Q. 82C 83C B REAL(K) 84C B contains the solution of the least squares problem 85C 86C minimize norm2(Y - XK*B), 87C 88C if its computation has been requested. (Note that 89C if pivoting was requested in SQRDC, the J-th 90C component of B will be associated with column JPVT(J) 91C of the original matrix X that was input into SQRDC.) 92C 93C RSD REAL(N). 94C RSD contains the least squares residual Y - XK*B, 95C if its computation has been requested. RSD is 96C also the orthogonal projection of Y onto the 97C orthogonal complement of the column space of XK. 98C 99C XB REAL(N). 100C XB contains the least squares approximation XK*B, 101C if its computation has been requested. XB is also 102C the orthogonal projection of Y onto the column space 103C of X. 104C 105C INFO INTEGER. 106C INFO is zero unless the computation of B has 107C been requested and R is exactly singular. In 108C this case, INFO is the index of the first zero 109C diagonal element of R and B is left unaltered. 110C 111C The parameters QY, QTY, B, RSD, and XB are not referenced 112C if their computation is not requested and in this case 113C can be replaced by dummy variables in the calling program. 114C To save storage, the user may in some cases use the same 115C array for different parameters in the calling sequence. A 116C frequently occuring example is when one wishes to compute 117C any of B, RSD, or XB and does not need Y or QTY. In this 118C case one may identify Y, QTY, and one of B, RSD, or XB, while 119C providing separate arrays for anything else that is to be 120C computed. Thus the calling sequence 121C 122C CALL SQRSL(X,LDX,N,K,QRAUX,Y,DUM,Y,B,Y,DUM,110,INFO) 123C 124C will result in the computation of B and RSD, with RSD 125C overwriting Y. More generally, each item in the following 126C list contains groups of permissible identifications for 127C a single callinng sequence. 128C 129C 1. (Y,QTY,B) (RSD) (XB) (QY) 130C 131C 2. (Y,QTY,RSD) (B) (XB) (QY) 132C 133C 3. (Y,QTY,XB) (B) (RSD) (QY) 134C 135C 4. (Y,QY) (QTY,B) (RSD) (XB) 136C 137C 5. (Y,QY) (QTY,RSD) (B) (XB) 138C 139C 6. (Y,QY) (QTY,XB) (B) (RSD) 140C 141C In any group the value returned in the array allocated to 142C the group corresponds to the last member of the group. 143C 144C LINPACK. This version dated 08/14/78 . 145C G. W. Stewart, University of Maryland, Argonne National Lab. 146C 147C SQRSL uses the following functions and subprograms. 148C 149C BLAS SAXPY,SCOPY,SDOT 150C Fortran ABS,MIN0,MOD 151C***REFERENCES DONGARRA J.J., BUNCH J.R., MOLER C.B., STEWART G.W., 152C *LINPACK USERS GUIDE*, SIAM, 1979. 153C***ROUTINES CALLED SAXPY,SCOPY,SDOT 154C***END PROLOGUE SQRSL 155 INTEGER LDX,N,K,JOB,INFO 156 REAL X(LDX,*),QRAUX(*),Y(*),QY(*),QTY(*),B(*),RSD(*),XB(*) 157C 158 INTEGER I,J,JJ,JU,KP1 159 REAL SDOT,T,TEMP 160 LOGICAL CB,CQY,CQTY,CR,CXB 161C 162C SET INFO FLAG. 163C 164C***FIRST EXECUTABLE STATEMENT SQRSL 165 INFO = 0 166C 167C DETERMINE WHAT IS TO BE COMPUTED. 168C 169 CQY = JOB/10000 .NE. 0 170 CQTY = MOD(JOB,10000) .NE. 0 171 CB = MOD(JOB,1000)/100 .NE. 0 172 CR = MOD(JOB,100)/10 .NE. 0 173 CXB = MOD(JOB,10) .NE. 0 174 JU = MIN0(K,N-1) 175C 176C SPECIAL ACTION WHEN N=1. 177C 178 IF (JU .NE. 0) GO TO 40 179 IF (CQY) QY(1) = Y(1) 180 IF (CQTY) QTY(1) = Y(1) 181 IF (CXB) XB(1) = Y(1) 182 IF (.NOT.CB) GO TO 30 183 IF (X(1,1) .NE. 0.0E0) GO TO 10 184 INFO = 1 185 GO TO 20 186 10 CONTINUE 187 B(1) = Y(1)/X(1,1) 188 20 CONTINUE 189 30 CONTINUE 190 IF (CR) RSD(1) = 0.0E0 191 GO TO 250 192 40 CONTINUE 193C 194C SET UP TO COMPUTE QY OR QTY. 195C 196 IF (CQY) CALL SCOPY(N,Y,1,QY,1) 197 IF (CQTY) CALL SCOPY(N,Y,1,QTY,1) 198 IF (.NOT.CQY) GO TO 70 199C 200C COMPUTE QY. 201C 202 DO 60 JJ = 1, JU 203 J = JU - JJ + 1 204 IF (QRAUX(J) .EQ. 0.0E0) GO TO 50 205 TEMP = X(J,J) 206 X(J,J) = QRAUX(J) 207 T = -SDOT(N-J+1,X(J,J),1,QY(J),1)/X(J,J) 208 CALL SAXPY(N-J+1,T,X(J,J),1,QY(J),1) 209 X(J,J) = TEMP 210 50 CONTINUE 211 60 CONTINUE 212 70 CONTINUE 213 IF (.NOT.CQTY) GO TO 100 214C 215C COMPUTE TRANS(Q)*Y. 216C 217 DO 90 J = 1, JU 218 IF (QRAUX(J) .EQ. 0.0E0) GO TO 80 219 TEMP = X(J,J) 220 X(J,J) = QRAUX(J) 221 T = -SDOT(N-J+1,X(J,J),1,QTY(J),1)/X(J,J) 222 CALL SAXPY(N-J+1,T,X(J,J),1,QTY(J),1) 223 X(J,J) = TEMP 224 80 CONTINUE 225 90 CONTINUE 226 100 CONTINUE 227C 228C SET UP TO COMPUTE B, RSD, OR XB. 229C 230 IF (CB) CALL SCOPY(K,QTY,1,B,1) 231 KP1 = K + 1 232 IF (CXB) CALL SCOPY(K,QTY,1,XB,1) 233 IF (CR .AND. K .LT. N) CALL SCOPY(N-K,QTY(KP1),1,RSD(KP1),1) 234 IF (.NOT.CXB .OR. KP1 .GT. N) GO TO 120 235 DO 110 I = KP1, N 236 XB(I) = 0.0E0 237 110 CONTINUE 238 120 CONTINUE 239 IF (.NOT.CR) GO TO 140 240 DO 130 I = 1, K 241 RSD(I) = 0.0E0 242 130 CONTINUE 243 140 CONTINUE 244 IF (.NOT.CB) GO TO 190 245C 246C COMPUTE B. 247C 248 DO 170 JJ = 1, K 249 J = K - JJ + 1 250 IF (X(J,J) .NE. 0.0E0) GO TO 150 251 INFO = J 252C ......EXIT 253 GO TO 180 254 150 CONTINUE 255 B(J) = B(J)/X(J,J) 256 IF (J .EQ. 1) GO TO 160 257 T = -B(J) 258 CALL SAXPY(J-1,T,X(1,J),1,B,1) 259 160 CONTINUE 260 170 CONTINUE 261 180 CONTINUE 262 190 CONTINUE 263 IF (.NOT.CR .AND. .NOT.CXB) GO TO 240 264C 265C COMPUTE RSD OR XB AS REQUIRED. 266C 267 DO 230 JJ = 1, JU 268 J = JU - JJ + 1 269 IF (QRAUX(J) .EQ. 0.0E0) GO TO 220 270 TEMP = X(J,J) 271 X(J,J) = QRAUX(J) 272 IF (.NOT.CR) GO TO 200 273 T = -SDOT(N-J+1,X(J,J),1,RSD(J),1)/X(J,J) 274 CALL SAXPY(N-J+1,T,X(J,J),1,RSD(J),1) 275 200 CONTINUE 276 IF (.NOT.CXB) GO TO 210 277 T = -SDOT(N-J+1,X(J,J),1,XB(J),1)/X(J,J) 278 CALL SAXPY(N-J+1,T,X(J,J),1,XB(J),1) 279 210 CONTINUE 280 X(J,J) = TEMP 281 220 CONTINUE 282 230 CONTINUE 283 240 CONTINUE 284 250 CONTINUE 285 RETURN 286 END 287