1c 2c dqrdc uses householder transformations to compute the qr 3c factorization of an n by p matrix x. column pivoting 4c based on the 2-norms of the reduced columns may be 5c performed at the users option. 6c 7c on entry 8c 9c x double precision(ldx,p), where ldx .ge. n. 10c x contains the matrix whose decomposition is to be 11c computed. 12c 13c ldx integer. 14c ldx is the leading dimension of the array x. 15c 16c n integer. 17c n is the number of rows of the matrix x. 18c 19c p integer. 20c p is the number of columns of the matrix x. 21c 22c jpvt integer(p). 23c jpvt contains integers that control the selection 24c of the pivot columns. the k-th column x(k) of x 25c is placed in one of three classes according to the 26c value of jpvt(k). 27c 28c if jpvt(k) .gt. 0, then x(k) is an initial 29c column. 30c 31c if jpvt(k) .eq. 0, then x(k) is a free column. 32c 33c if jpvt(k) .lt. 0, then x(k) is a final column. 34c 35c before the decomposition is computed, initial columns 36c are moved to the beginning of the array x and final 37c columns to the end. both initial and final columns 38c are frozen in place during the computation and only 39c free columns are moved. at the k-th stage of the 40c reduction, if x(k) is occupied by a free column 41c it is interchanged with the free column of largest 42c reduced norm. jpvt is not referenced if 43c job .eq. 0. 44c 45c work double precision(p). 46c work is a work array. work is not referenced if 47c job .eq. 0. 48c 49c job integer. 50c job is an integer that initiates column pivoting. 51c if job .eq. 0, no pivoting is done. 52c if job .ne. 0, pivoting is done. 53c 54c on return 55c 56c x x contains in its upper triangle the upper 57c triangular matrix r of the qr factorization. 58c below its diagonal x contains information from 59c which the orthogonal part of the decomposition 60c can be recovered. note that if pivoting has 61c been requested, the decomposition is not that 62c of the original matrix x but that of x 63c with its columns permuted as described by jpvt. 64c 65c qraux double precision(p). 66c qraux contains further information required to recover 67c the orthogonal part of the decomposition. 68c 69c jpvt jpvt(k) contains the index of the column of the 70c original matrix that has been interchanged into 71c the k-th column, if pivoting was requested. 72c 73c linpack. this version dated 08/14/78 . 74c g.w. stewart, university of maryland, argonne national lab. 75c 76c dqrdc uses the following functions and subprograms. 77c 78c blas daxpy,ddot,dscal,dswap,dnrm2 79c fortran dabs,dmax1,min0,dsqrt 80c 81 subroutine dqrdc(x,ldx,n,p,qraux,jpvt,work,job) 82 integer ldx,n,p,job 83 integer jpvt(p) 84 double precision x(ldx,p),qraux(p),work(p) 85c 86c internal variables 87c 88 integer j,jp,jj, l,lp1,lup,maxj,pl,pu 89 double precision maxnrm,dnrm2,tt 90 double precision ddot,nrmxl,t 91 logical negj,swapj 92c 93c 94 pl = 1 95 pu = 0 96 if (job .eq. 0) go to 60 97c 98c pivoting has been requested. rearrange the columns 99c according to jpvt. 100c 101 do 20 j = 1, p 102 swapj = jpvt(j) .gt. 0 103 negj = jpvt(j) .lt. 0 104 jpvt(j) = j 105 if (negj) jpvt(j) = -j 106 if (.not.swapj) go to 10 107 if (j .ne. pl) call dswap(n,x(1,pl),1,x(1,j),1) 108 jpvt(j) = jpvt(pl) 109 jpvt(pl) = j 110 pl = pl + 1 111 10 continue 112 20 continue 113 pu = p 114 do 50 jj = 1, p 115 j = p - jj + 1 116 if (jpvt(j) .ge. 0) go to 40 117 jpvt(j) = -jpvt(j) 118 if (j .eq. pu) go to 30 119 call dswap(n,x(1,pu),1,x(1,j),1) 120 jp = jpvt(pu) 121 jpvt(pu) = jpvt(j) 122 jpvt(j) = jp 123 30 continue 124 pu = pu - 1 125 40 continue 126 50 continue 127 60 continue 128c 129c compute the norms of the free columns. 130c 131 if (pu .lt. pl) go to 80 132 do 70 j = pl, pu 133 qraux(j) = dnrm2(n,x(1,j),1) 134 work(j) = qraux(j) 135 70 continue 136 80 continue 137c 138c perform the householder reduction of x. 139c 140 lup = min0(n,p) 141 do 200 l = 1, lup 142 if (l .lt. pl .or. l .ge. pu) go to 120 143c 144c locate the column of largest norm and bring it 145c into the pivot position. 146c 147 maxnrm = 0.0d0 148 maxj = l 149 do 100 j = l, pu 150 if (qraux(j) .le. maxnrm) go to 90 151 maxnrm = qraux(j) 152 maxj = j 153 90 continue 154 100 continue 155 if (maxj .eq. l) go to 110 156 call dswap(n,x(1,l),1,x(1,maxj),1) 157 qraux(maxj) = qraux(l) 158 work(maxj) = work(l) 159 jp = jpvt(maxj) 160 jpvt(maxj) = jpvt(l) 161 jpvt(l) = jp 162 110 continue 163 120 continue 164 qraux(l) = 0.0d0 165 if (l .eq. n) go to 190 166c 167c compute the householder transformation for column l. 168c 169 nrmxl = dnrm2(n-l+1,x(l,l),1) 170 if (nrmxl .eq. 0.0d0) go to 180 171 if (x(l,l) .ne. 0.0d0) nrmxl = dsign(nrmxl,x(l,l)) 172 call dscal(n-l+1,1.0d0/nrmxl,x(l,l),1) 173 x(l,l) = 1.0d0 + x(l,l) 174c 175c apply the transformation to the remaining columns, 176c updating the norms. 177c 178 lp1 = l + 1 179 if (p .lt. lp1) go to 170 180 do 160 j = lp1, p 181 t = -ddot(n-l+1,x(l,l),1,x(l,j),1)/x(l,l) 182 call daxpy(n-l+1,t,x(l,l),1,x(l,j),1) 183 if (j .lt. pl .or. j .gt. pu) go to 150 184 if (qraux(j) .eq. 0.0d0) go to 150 185 tt = 1.0d0 - (dabs(x(l,j))/qraux(j))**2 186 tt = dmax1(tt,0.0d0) 187 t = tt 188 tt = 1.0d0 + 0.05d0*tt*(qraux(j)/work(j))**2 189 if (tt .eq. 1.0d0) go to 130 190 qraux(j) = qraux(j)*dsqrt(t) 191 go to 140 192 130 continue 193 qraux(j) = dnrm2(n-l,x(l+1,j),1) 194 work(j) = qraux(j) 195 140 continue 196 150 continue 197 160 continue 198 170 continue 199c 200c save the transformation. 201c 202 qraux(l) = x(l,l) 203 x(l,l) = -nrmxl 204 180 continue 205 190 continue 206 200 continue 207 return 208 end 209