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