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