| c |
| c dqrsl applies the output of dqrdc to compute coordinate |
| c transformations, projections, and least squares solutions. |
| c for k .le. min(n,p), let xk be the matrix |
| c |
| c xk = (x(jpvt(1)),x(jpvt(2)), ... ,x(jpvt(k))) |
| c |
| c formed from columns jpvt(1), ... ,jpvt(k) of the original |
| c n x p matrix x that was input to dqrdc (if no pivoting was |
| c done, xk consists of the first k columns of x in their |
| c original order). dqrdc produces a factored orthogonal matrix q |
| c and an upper triangular matrix r such that |
| c |
| c xk = q * (r) |
| c (0) |
| c |
| c this information is contained in coded form in the arrays |
| c x and qraux. |
| c |
| c on entry |
| c |
| c x double precision(ldx,p). |
| c x contains the output of dqrdc. |
| c |
| c ldx integer. |
| c ldx is the leading dimension of the array x. |
| c |
| c n integer. |
| c n is the number of rows of the matrix xk. it must |
| c have the same value as n in dqrdc. |
| c |
| c k integer. |
| c k is the number of columns of the matrix xk. k |
| c must not be greater than min(n,p), where p is the |
| c same as in the calling sequence to dqrdc. |
| c |
| c qraux double precision(p). |
| c qraux contains the auxiliary output from dqrdc. |
| c |
| c y double precision(n) |
| c y contains an n-vector that is to be manipulated |
| c by dqrsl. |
| c |
| c job integer. |
| c job specifies what is to be computed. job has |
| c the decimal expansion abcde, with the following |
| c meaning. |
| c |
| c if a.ne.0, compute qy. |
| c if b,c,d, or e .ne. 0, compute qty. |
| c if c.ne.0, compute b. |
| c if d.ne.0, compute rsd. |
| c if e.ne.0, compute xb. |
| c |
| c note that a request to compute b, rsd, or xb |
| c automatically triggers the computation of qty, for |
| c which an array must be provided in the calling |
| c sequence. |
| c |
| c on return |
| c |
| c qy double precision(n). |
| c qy contains q*y, if its computation has been |
| c requested. |
| c |
| c qty double precision(n). |
| c qty contains trans(q)*y, if its computation has |
| c been requested. here trans(q) is the |
| c transpose of the matrix q. |
| c |
| c b double precision(k) |
| c b contains the solution of the least squares problem |
| c |
| c minimize norm2(y - xk*b), |
| c |
| c if its computation has been requested. (note that |
| c if pivoting was requested in dqrdc, the j-th |
| c component of b will be associated with column jpvt(j) |
| c of the original matrix x that was input into dqrdc.) |
| c |
| c rsd double precision(n). |
| c rsd contains the least squares residual y - xk*b, |
| c if its computation has been requested. rsd is |
| c also the orthogonal projection of y onto the |
| c orthogonal complement of the column space of xk. |
| c |
| c xb double precision(n). |
| c xb contains the least squares approximation xk*b, |
| c if its computation has been requested. xb is also |
| c the orthogonal projection of y onto the column space |
| c of x. |
| c |
| c info integer. |
| c info is zero unless the computation of b has |
| c been requested and r is exactly singular. in |
| c this case, info is the index of the first zero |
| c diagonal element of r and b is left unaltered. |
| c |
| c the parameters qy, qty, b, rsd, and xb are not referenced |
| c if their computation is not requested and in this case |
| c can be replaced by dummy variables in the calling program. |
| c to save storage, the user may in some cases use the same |
| c array for different parameters in the calling sequence. a |
| c frequently occuring example is when one wishes to compute |
| c any of b, rsd, or xb and does not need y or qty. in this |
| c case one may identify y, qty, and one of b, rsd, or xb, while |
| c providing separate arrays for anything else that is to be |
| c computed. thus the calling sequence |
| c |
| c call dqrsl(x,ldx,n,k,qraux,y,dum,y,b,y,dum,110,info) |
| c |
| c will result in the computation of b and rsd, with rsd |
| c overwriting y. more generally, each item in the following |
| c list contains groups of permissible identifications for |
| c a single calling sequence. |
| c |
| c 1. (y,qty,b) (rsd) (xb) (qy) |
| c |
| c 2. (y,qty,rsd) (b) (xb) (qy) |
| c |
| c 3. (y,qty,xb) (b) (rsd) (qy) |
| c |
| c 4. (y,qy) (qty,b) (rsd) (xb) |
| c |
| c 5. (y,qy) (qty,rsd) (b) (xb) |
| c |
| c 6. (y,qy) (qty,xb) (b) (rsd) |
| c |
| c in any group the value returned in the array allocated to |
| c the group corresponds to the last member of the group. |
| c |
| c linpack. this version dated 08/14/78 . |
| c g.w. stewart, university of maryland, argonne national lab. |
| c |
| c dqrsl uses the following functions and subprograms. |
| c |
| c BLAS daxpy,dcopy,ddot |
| c Fortran dabs,min0,mod |
| c |
| subroutine dqrsl(x,ldx,n,k,qraux,y,qy,qty,b,rsd,xb,job,info) |
| integer ldx,n,k,job,info |
| double precision x(ldx,k),qraux(k),y(n),qy(n),qty(n),b(k),rsd(n), |
| * xb(n) |
| c |
| c internal variables |
| c |
| integer i,j,jj,ju,kp1 |
| double precision ddot,t,temp |
| logical cb,cqy,cqty,cr,cxb |
| c |
| c |
| c set info flag. |
| c |
| info = 0 |
| c |
| c determine what is to be computed. |
| c |
| cqy = job/10000 .ne. 0 |
| cqty = mod(job,10000) .ne. 0 |
| cb = mod(job,1000)/100 .ne. 0 |
| cr = mod(job,100)/10 .ne. 0 |
| cxb = mod(job,10) .ne. 0 |
| ju = min0(k,n-1) |
| c |
| c special action when n=1. |
| c |
| if (ju .ne. 0) go to 40 |
| if (cqy) qy(1) = y(1) |
| if (cqty) qty(1) = y(1) |
| if (cxb) xb(1) = y(1) |
| if (.not.cb) go to 30 |
| if (x(1,1) .ne. 0.0d0) go to 10 |
| info = 1 |
| go to 20 |
| 10 continue |
| b(1) = y(1)/x(1,1) |
| 20 continue |
| 30 continue |
| if (cr) rsd(1) = 0.0d0 |
| go to 250 |
| 40 continue |
| c |
| c set up to compute qy or qty. |
| c |
| if (cqy) call dcopy(n,y,1,qy,1) |
| if (cqty) call dcopy(n,y,1,qty,1) |
| if (.not.cqy) go to 70 |
| c |
| c compute qy. |
| c |
| do 60 jj = 1, ju |
| j = ju - jj + 1 |
| if (qraux(j) .eq. 0.0d0) go to 50 |
| temp = x(j,j) |
| x(j,j) = qraux(j) |
| t = -ddot(n-j+1,x(j,j),1,qy(j),1)/x(j,j) |
| call daxpy(n-j+1,t,x(j,j),1,qy(j),1) |
| x(j,j) = temp |
| 50 continue |
| 60 continue |
| 70 continue |
| if (.not.cqty) go to 100 |
| c |
| c compute trans(q)*y. |
| c |
| do 90 j = 1, ju |
| if (qraux(j) .eq. 0.0d0) go to 80 |
| temp = x(j,j) |
| x(j,j) = qraux(j) |
| t = -ddot(n-j+1,x(j,j),1,qty(j),1)/x(j,j) |
| call daxpy(n-j+1,t,x(j,j),1,qty(j),1) |
| x(j,j) = temp |
| 80 continue |
| 90 continue |
| 100 continue |
| c |
| c set up to compute b, rsd, or xb. |
| c |
| if (cb) call dcopy(k,qty,1,b,1) |
| kp1 = k + 1 |
| if (cxb) call dcopy(k,qty,1,xb,1) |
| if (cr .and. k .lt. n) call dcopy(n-k,qty(kp1),1,rsd(kp1),1) |
| if (.not.cxb .or. kp1 .gt. n) go to 120 |
| do 110 i = kp1, n |
| xb(i) = 0.0d0 |
| 110 continue |
| 120 continue |
| if (.not.cr) go to 140 |
| do 130 i = 1, k |
| rsd(i) = 0.0d0 |
| 130 continue |
| 140 continue |
| if (.not.cb) go to 190 |
| c |
| c compute b. |
| c |
| do 170 jj = 1, k |
| j = k - jj + 1 |
| if (x(j,j) .ne. 0.0d0) go to 150 |
| info = j |
| c ......exit |
| go to 180 |
| 150 continue |
| b(j) = b(j)/x(j,j) |
| if (j .eq. 1) go to 160 |
| t = -b(j) |
| call daxpy(j-1,t,x(1,j),1,b,1) |
| 160 continue |
| 170 continue |
| 180 continue |
| 190 continue |
| if (.not.cr .and. .not.cxb) go to 240 |
| c |
| c compute rsd or xb as required. |
| c |
| do 230 jj = 1, ju |
| j = ju - jj + 1 |
| if (qraux(j) .eq. 0.0d0) go to 220 |
| temp = x(j,j) |
| x(j,j) = qraux(j) |
| if (.not.cr) go to 200 |
| t = -ddot(n-j+1,x(j,j),1,rsd(j),1)/x(j,j) |
| call daxpy(n-j+1,t,x(j,j),1,rsd(j),1) |
| 200 continue |
| if (.not.cxb) go to 210 |
| t = -ddot(n-j+1,x(j,j),1,xb(j),1)/x(j,j) |
| call daxpy(n-j+1,t,x(j,j),1,xb(j),1) |
| 210 continue |
| x(j,j) = temp |
| 220 continue |
| 230 continue |
| 240 continue |
| 250 continue |
| return |
| end |
