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