1 subroutine sqrdc(x,ldx,n,p,qraux,jpvt,work,job) 2 integer ldx,n,p,job 3 integer jpvt(1) 4 real x(ldx,1),qraux(1),work(1) 5c 6c sqrdc 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 real(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 real(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 orthogonal 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 real(p). 70c qraux contains further information required to recover 71c the orthogonal 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 sqrdc uses the following functions and subprograms. 81c 82c blas saxpy,sdot,sscal,sswap,snrm2 83c fortran abs,amax1,min0,sqrt 84c 85c internal variables 86c 87 integer j,jp,l,lp1,lup,maxj,pl,pu 88 real maxnrm,snrm2,tt 89 real sdot,nrmxl,t 90 logical negj,swapj 91c 92c 93 pl = 1 94 pu = 0 95 if (job .eq. 0) go to 60 96c 97c pivoting has been requested. rearrange the columns 98c according to jpvt. 99c 100 do 20 j = 1, p 101 swapj = jpvt(j) .gt. 0 102 negj = jpvt(j) .lt. 0 103 jpvt(j) = j 104 if (negj) jpvt(j) = -j 105 if (.not.swapj) go to 10 106 if (j .ne. pl) call sswap(n,x(1,pl),1,x(1,j),1) 107 jpvt(j) = jpvt(pl) 108 jpvt(pl) = j 109 pl = pl + 1 110 10 continue 111 20 continue 112 pu = p 113 do 50 jj = 1, p 114 j = p - jj + 1 115 if (jpvt(j) .ge. 0) go to 40 116 jpvt(j) = -jpvt(j) 117 if (j .eq. pu) go to 30 118 call sswap(n,x(1,pu),1,x(1,j),1) 119 jp = jpvt(pu) 120 jpvt(pu) = jpvt(j) 121 jpvt(j) = jp 122 30 continue 123 pu = pu - 1 124 40 continue 125 50 continue 126 60 continue 127c 128c compute the norms of the free columns. 129c 130 if (pu .lt. pl) go to 80 131 do 70 j = pl, pu 132 qraux(j) = snrm2(n,x(1,j),1) 133 work(j) = qraux(j) 134 70 continue 135 80 continue 136c 137c perform the householder reduction of x. 138c 139 lup = min0(n,p) 140 do 200 l = 1, lup 141 if (l .lt. pl .or. l .ge. pu) go to 120 142c 143c locate the column of largest norm and bring it 144c into the pivot position. 145c 146 maxnrm = 0.0e0 147 maxj = l 148 do 100 j = l, pu 149 if (qraux(j) .le. maxnrm) go to 90 150 maxnrm = qraux(j) 151 maxj = j 152 90 continue 153 100 continue 154 if (maxj .eq. l) go to 110 155 call sswap(n,x(1,l),1,x(1,maxj),1) 156 qraux(maxj) = qraux(l) 157 work(maxj) = work(l) 158 jp = jpvt(maxj) 159 jpvt(maxj) = jpvt(l) 160 jpvt(l) = jp 161 110 continue 162 120 continue 163 qraux(l) = 0.0e0 164 if (l .eq. n) go to 190 165c 166c compute the householder transformation for column l. 167c 168 nrmxl = snrm2(n-l+1,x(l,l),1) 169 if (nrmxl .eq. 0.0e0) go to 180 170 if (x(l,l) .ne. 0.0e0) nrmxl = sign(nrmxl,x(l,l)) 171 call sscal(n-l+1,1.0e0/nrmxl,x(l,l),1) 172 x(l,l) = 1.0e0 + x(l,l) 173c 174c apply the transformation to the remaining columns, 175c updating the norms. 176c 177 lp1 = l + 1 178 if (p .lt. lp1) go to 170 179 do 160 j = lp1, p 180 t = -sdot(n-l+1,x(l,l),1,x(l,j),1)/x(l,l) 181 call saxpy(n-l+1,t,x(l,l),1,x(l,j),1) 182 if (j .lt. pl .or. j .gt. pu) go to 150 183 if (qraux(j) .eq. 0.0e0) go to 150 184 tt = 1.0e0 - (abs(x(l,j))/qraux(j))**2 185 tt = amax1(tt,0.0e0) 186 t = tt 187 tt = 1.0e0 + 0.05e0*tt*(qraux(j)/work(j))**2 188 if (tt .eq. 1.0e0) go to 130 189 qraux(j) = qraux(j)*sqrt(t) 190 go to 140 191 130 continue 192 qraux(j) = snrm2(n-l,x(l+1,j),1) 193 work(j) = qraux(j) 194 140 continue 195 150 continue 196 160 continue 197 170 continue 198c 199c save the transformation. 200c 201 qraux(l) = x(l,l) 202 x(l,l) = -nrmxl 203 180 continue 204 190 continue 205 200 continue 206 return 207 end 208