Received: (from mfrancl@localhost) by kelvin.brynmawr.edu (8.6.12/8.6.12) id LAA02523 for srusso@cc; Fri, 19 May 1995 11:06:00 -0400 Date: Fri, 19 May 1995 11:06:00 -0400
From: Francl Michelle M 
Message-Id: <199505191506.LAA02523@kelvin.brynmawr.edu> To: srusso
X-UIDL: 800896205.001

cc
cc	a calling routine to do QR with column pivoting on the chelp matrix
cc
implicit real*8 (a-h,o-z)
dimension x(100,100),qraux(100),work(100),jpvt(100) c
c	for job.ne.0 then pivoting is done on all columns of x where
c	jpvt(i) is 0
c	in this case, don't pivot through last col, it's the constraint!
c	so make it the first column
c
read(5,*) ldx
job=1
c
read(5,*) nrows,ncol
do 100 k=1,nrows
read(5,*) (x(k,mu),mu=1,ncol)
100 continue
c
c	set up the pivoting constraints, make all col's free except
c	for the one corresponding to the constraint, which should be
c	kept!
c
do 150 k=1,nrows-1
jpvt(k) = 0
150 continue
jpvt(nrows) = 1
c
call dqrdc(x,ldx,nrows,ncol,qraux,jpvt,work,job) c
write(6,*) 'r matrix'
do 200 k=1,nrows
write(6,1000) (x(k,mu),mu=1,nrows)
1000 format(1x,10f10.4)
200 continue
c
write(6,*) 'pivots list'
write(6,*) (jpvt(k),k=1,ncol)
c
write(6,*) 'QR aux'
write(6,1000) (qraux(mu),mu=1,ncol)
c
end
subroutine dqrdc(x,ldx,n,p,qraux,jpvt,work,job) integer ldx,n,p,job
integer jpvt(1)
double precision x(ldx,1),qraux(1),work(1) c
c	dqrdc uses householder transformations to compute the qr
c	factorization of an n by p matrix x. column pivoting
c	based on the 2-norms of the reduced columns may be
c	performed at the users option.
c
c	on entry
c
c	x	double precision(ldx,p), where ldx .ge. n.
c	x contains the matrix whose decomposition is to be
c	computed.
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 x.
c
c	p	integer.
c	p is the number of columns of the matrix x.
c
c	jpvt integer(p).
c	jpvt contains integers that control the selection
c	of the pivot columns. the k-th column x(k) of x
c	is placed in one of three classes according to the
c	value of jpvt(k).
c
c	if jpvt(k) .gt. 0, then x(k) is an initial
c	column.
c
c	if jpvt(k) .eq. 0, then x(k) is a free column.
c
c	if jpvt(k) .lt. 0, then x(k) is a final column.
c
c	before the decomposition is computed, initial columns
c	are moved to the beginning of the array x and final
c	columns to the end. both initial and final columns
c	are frozen in place during the computation and only
c	free columns are moved. at the k-th stage of the
c	reduction, if x(k) is occupied by a free column
c	it is interchanged with the free column of largest
c	reduced norm. jpvt is not referenced if
c	job .eq. 0.
c
c	work double precision(p).
c	work is a work array. work is not referenced if
c	job .eq. 0.
c
c	job	integer.
c	job is an integer that initiates column pivoting.
c	if job .eq. 0, no pivoting is done.
c	if job .ne. 0, pivoting is done.
c
c	on return
c
c	x	x contains in its upper triangle the upper
c	triangular matrix r of the qr factorization.
c	below its diagonal x contains information from
c	which the orthogonal part of the decomposition
c	can be recovered. note that if pivoting has
c	been requested, the decomposition is not that
c	of the original matrix x but that of x
c	with its columns permuted as described by jpvt.
c
c	qraux double precision(p).
c	qraux contains further information required to recover
c	the orthogonal part of the decomposition.
c
c	jpvt jpvt(k) contains the index of the column of the
c	original matrix that has been interchanged into
c	the k-th column, if pivoting was requested.
c
c	linpack. this version dated 08/14/78 .
c	g.w. stewart, university of maryland, argonne national lab.
c
c	dqrdc uses the following functions and subprograms.
c
c	blas daxpy,ddot,dscal,dswap,dnrm2
c	fortran dabs,dmax1,min0,dsqrt
c
c	internal variables
c
integer j,jp,l,lp1,lup,maxj,pl,pu
double precision maxnrm,dnrm2,tt
double precision ddot,nrmxl,t
logical negj,swapj
c
c
pl = 1
pu = 0
if (job .eq. 0) go to 60
c
c	pivoting has been requested. rearrange the columns
c	according to jpvt.
c
do 20 j = 1, p
swapj = jpvt(j) .gt. 0
negj = jpvt(j) .lt. 0
jpvt(j) = j
if (negj) jpvt(j) = -j
if (.not.swapj) go to 10
if (j .ne. pl) call dswap(n,x(1,pl),1,x(1,j),1) jpvt(j) = jpvt(pl)
jpvt(pl) = j
pl = pl + 1
10	continue
20 continue
pu = p
do 50 jj = 1, p
j = p - jj + 1
if (jpvt(j) .ge. 0) go to 40
jpvt(j) = -jpvt(j)
if (j .eq. pu) go to 30
call dswap(n,x(1,pu),1,x(1,j),1)
jp = jpvt(pu)
jpvt(pu) = jpvt(j)
jpvt(j) = jp
30	continue
pu = pu - 1
40	continue
50 continue
60 continue
c
c	compute the norms of the free columns.
c
if (pu .lt. pl) go to 80
do 70 j = pl, pu
qraux(j) = dnrm2(n,x(1,j),1)
work(j) = qraux(j)
70 continue
80 continue
c
c	perform the householder reduction of x.
c
lup = min0(n,p)
do 200 l = 1, lup
if (l .lt. pl .or. l .ge. pu) go to 120
c
c	locate the column of largest norm and bring it
c	into the pivot position.
c
maxnrm = 0.0d0
maxj = l
do 100 j = l, pu
if (qraux(j) .le. maxnrm) go to 90
maxnrm = qraux(j)
maxj = j
90	continue
100	continue
if (maxj .eq. l) go to 110
call dswap(n,x(1,l),1,x(1,maxj),1)
qraux(maxj) = qraux(l)
work(maxj) = work(l)
jp = jpvt(maxj)
jpvt(maxj) = jpvt(l)
jpvt(l) = jp
110	continue
120 continue
qraux(l) = 0.0d0
if (l .eq. n) go to 190
c
c	compute the householder transformation for column l.
c
nrmxl = dnrm2(n-l+1,x(l,l),1)
if (nrmxl .eq. 0.0d0) go to 180
if (x(l,l) .ne. 0.0d0) nrmxl = dsign(nrmxl,x(l,l)) call dscal(n-l+1,1.0d0/nrmxl,x(l,l),1)
x(l,l) = 1.0d0 + x(l,l)
c
c	apply the transformation to the remaining columns,
c	updating the norms.
c
lp1 = l + 1
if (p .lt. lp1) go to 170
do 160 j = lp1, p
t = -ddot(n-l+1,x(l,l),1,x(l,j),1)/x(l,l) call daxpy(n-l+1,t,x(l,l),1,x(l,j),1)
if (j .lt. pl .or. j .gt. pu) go to 150
if (qraux(j) .eq. 0.0d0) go to 150
tt = 1.0d0 - (dabs(x(l,j))/qraux(j))**2
tt = dmax1(tt,0.0d0)
t = tt
tt = 1.0d0 + 0.05d0*tt*(qraux(j)/work(j))**2
if (tt .eq. 1.0d0) go to 130
qraux(j) = qraux(j)*dsqrt(t)
go to 140
130	continue
qraux(j) = dnrm2(n-l,x(l+1,j),1)
work(j) = qraux(j)
140	continue
150	continue
160	continue
170	continue
c
c	save the transformation.
c
qraux(l) = x(l,l)
x(l,l) = -nrmxl
180	continue
190 continue
200 continue
return
end
subroutine daxpy(n,da,dx,incx,dy,incy)
c
c	constant times a vector plus a vector.
c	uses unrolled loops for increments equal to one.
c	jack dongarra, linpack, 3/11/78.
c
double precision dx(1),dy(1),da
integer i,incx,incy,ix,iy,m,mp1,n
c
if(n.le.0)return
if (da .eq. 0.0d0) return
if(incx.eq.1.and.incy.eq.1)go to 20
c
c	code for unequal increments or equal increments
c	not equal to 1
c
ix = 1
iy = 1
if(incx.lt.0)ix = (-n+1)*incx + 1
if(incy.lt.0)iy = (-n+1)*incy + 1
do 10 i = 1,n
dy(iy) = dy(iy) + da*dx(ix)
ix = ix + incx
iy = iy + incy
10 continue
return
c
c	code for both increments equal to 1
c
c
c	clean-up loop
c
20 m = mod(n,4)
if( m .eq. 0 ) go to 40
do 30 i = 1,m
dy(i) = dy(i) + da*dx(i)
30 continue
if( n .lt. 4 ) return
40 mp1 = m + 1
do 50 i = mp1,n,4
dy(i) = dy(i) + da*dx(i)
dy(i + 1) = dy(i + 1) + da*dx(i + 1)
dy(i + 2) = dy(i + 2) + da*dx(i + 2)
dy(i + 3) = dy(i + 3) + da*dx(i + 3)
50 continue
return
end
double precision function ddot(n,dx,incx,dy,incy) c
c	forms the dot product of two vectors.
c	uses unrolled loops for increments equal to one.
c	jack dongarra, linpack, 3/11/78.
c
double precision dx(1),dy(1),dtemp
integer i,incx,incy,ix,iy,m,mp1,n
c
ddot = 0.0d0
dtemp = 0.0d0
if(n.le.0)return
if(incx.eq.1.and.incy.eq.1)go to 20
c
c	code for unequal increments or equal increments
c	not equal to 1
c
ix = 1
iy = 1
if(incx.lt.0)ix = (-n+1)*incx + 1
if(incy.lt.0)iy = (-n+1)*incy + 1
do 10 i = 1,n
dtemp = dtemp + dx(ix)*dy(iy)
ix = ix + incx
iy = iy + incy
10 continue
ddot = dtemp
return
c
c	code for both increments equal to 1
c
c
c	clean-up loop
c
20 m = mod(n,5)
if( m .eq. 0 ) go to 40
do 30 i = 1,m
dtemp = dtemp + dx(i)*dy(i)
30 continue
if( n .lt. 5 ) go to 60
40 mp1 = m + 1
do 50 i = mp1,n,5
dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) + * dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4) 50 continue
60 ddot = dtemp
return
end
subroutine dscal(n,da,dx,incx)
c
c	scales a vector by a constant.
c	uses unrolled loops for increment equal to one.
c	jack dongarra, linpack, 3/11/78.
c	modified 3/93 to return if incx .le. 0.
c
double precision da,dx(1)
integer i,incx,m,mp1,n,nincx
c
if( n.le.0 .or. incx.le.0 )return
if(incx.eq.1)go to 20
c
c	code for increment not equal to 1
c
nincx = n*incx
do 10 i = 1,nincx,incx
dx(i) = da*dx(i)
10 continue
return
c
c	code for increment equal to 1
c
c
c	clean-up loop
c
20 m = mod(n,5)
if( m .eq. 0 ) go to 40
do 30 i = 1,m
dx(i) = da*dx(i)
30 continue
if( n .lt. 5 ) return
40 mp1 = m + 1
do 50 i = mp1,n,5
dx(i) = da*dx(i)
dx(i + 1) = da*dx(i + 1)
dx(i + 2) = da*dx(i + 2)
dx(i + 3) = da*dx(i + 3)
dx(i + 4) = da*dx(i + 4)
50 continue
return
end
double precision function dnrm2 ( n, dx, incx) integer i, incx, ix, j, n, next
double precision dx(1), cutlo, cuthi, hitest, sum, xmax,zero,one data zero, one /0.0d0, 1.0d0/
c
c	euclidean norm of the n-vector stored in dx() with storage
c	increment incx .
c	if n .le. 0 return with result = 0.
c	if n .ge. 1 then incx must be .ge. 1
c
c	c.l.lawson, 1978 jan 08
c	modified to correct failure to update ix, 1/25/92.
c	modified 3/93 to return if incx .le. 0.
c
c	four phase method	using two built-in constants that are
c	hopefully applicable to all machines.
c	cutlo = maximum of dsqrt(u/eps) over all known machines.
c	cuthi = minimum of dsqrt(v)	over all known machines.
c	where
c	eps = smallest no. such that eps + 1. .gt. 1.
c	u = smallest positive no. (underflow limit)
c	v = largest no.	(overflow limit)
c
c	brief outline of algorithm..
c
c	phase 1 scans zero components.
c	move to phase 2 when a component is nonzero and .le. cutlo
c	move to phase 3 when a component is .gt. cutlo
c	move to phase 4 when a component is .ge. cuthi/m
c	where m = n for x() real and m = 2*n for complex.
c
c	values for cutlo and cuthi..
c	from the environmental parameters listed in the imsl converter
c	document the limiting values are as follows..
c	cutlo, s.p. u/eps = 2**(-102) for honeywell. close seconds are
c	univac and dec at 2**(-103)
c	thus cutlo = 2**(-51) = 4.44089e-16
c	cuthi, s.p. v = 2**127 for univac, honeywell, and dec.
c	thus cuthi = 2**(63.5) = 1.30438e19
c	cutlo, d.p. u/eps = 2**(-67) for honeywell and dec.
c	thus cutlo = 2**(-33.5) = 8.23181d-11
c	cuthi, d.p. same as s.p. cuthi = 1.30438d19
c	data cutlo, cuthi / 8.232d-11, 1.304d19 /
c	data cutlo, cuthi / 4.441e-16, 1.304e19 /
data cutlo, cuthi / 8.232d-11, 1.304d19 / c
if(n .gt. 0 .and. incx.gt.0) go to 10
dnrm2 = zero
go to 300
c
10 assign 30 to next
sum = zero
i = 1
ix = 1
c	begin main loop
20 go to next,(30, 50, 70, 110)
30 if( dabs(dx(i)) .gt. cutlo) go to 85
assign 50 to next
xmax = zero
c
c	phase 1. sum is zero
c
50 if( dx(i) .eq. zero) go to 200
if( dabs(dx(i)) .gt. cutlo) go to 85
c
c	prepare for phase 2.
assign 70 to next
go to 105
c
c	prepare for phase 4.
c
100 continue
ix = j
assign 110 to next
sum = (sum / dx(i)) / dx(i)
105 xmax = dabs(dx(i))
go to 115
c
c	phase 2. sum is small.
c	scale to avoid destructive underflow.
c
70 if( dabs(dx(i)) .gt. cutlo ) go to 75 c
c	common code for phases 2 and 4.
c	in phase 4 sum is large. scale to avoid overflow.
c
110 if( dabs(dx(i)) .le. xmax ) go to 115 
sum = one + sum * (xmax / dx(i))**2
xmax = dabs(dx(i))
go to 200
c
115 sum = sum + (dx(i)/xmax)**2
go to 200
c
c
c	prepare for phase 3.
c
75 sum = (sum * xmax) * xmax
c
c
c	for real or d.p. set hitest = cuthi/n
c	for complex	set hitest = cuthi/(2*n)
c
85 hitest = cuthi/float( n )
c
c	phase 3. sum is mid-range. no scaling.
c
do 95 j = ix,n
if(dabs(dx(i)) .ge. hitest) go to 100
sum = sum + dx(i)**2
i = i + incx
95 continue
dnrm2 = dsqrt( sum )
go to 300
c
200 continue
ix = ix + 1
i = i + incx
if( ix .le. n ) go to 20
c
c	end of main loop.
c
c	compute square root and adjust for scaling.
c
dnrm2 = xmax * dsqrt(sum)
300 continue
return
end
subroutine dswap (n,dx,incx,dy,incy)
c
c	interchanges two vectors.
c	uses unrolled loops for increments equal one.
c	jack dongarra, linpack, 3/11/78.
c
double precision dx(1),dy(1),dtemp
integer i,incx,incy,ix,iy,m,mp1,n
c
if(n.le.0)return
if(incx.eq.1.and.incy.eq.1)go to 20
c
c	code for unequal increments or equal increments not equal
c	to 1
c
ix = 1
iy = 1
if(incx.lt.0)ix = (-n+1)*incx + 1
if(incy.lt.0)iy = (-n+1)*incy + 1
do 10 i = 1,n
dtemp = dx(ix)
dx(ix) = dy(iy)
dy(iy) = dtemp
ix = ix + incx
iy = iy + incy
10 continue
return
c
c	code for both increments equal to 1
c
c
c	clean-up loop
c
20 m = mod(n,3)
if( m .eq. 0 ) go to 40
do 30 i = 1,m
dtemp = dx(i)
dx(i) = dy(i)
dy(i) = dtemp
30 continue
if( n .lt. 3 ) return
40 mp1 = m + 1
do 50 i = mp1,n,3
dtemp = dx(i)
dx(i) = dy(i)
dy(i) = dtemp
dtemp = dx(i + 1)
dx(i + 1) = dy(i + 1)
dy(i + 1) = dtemp
dtemp = dx(i + 2)
dx(i + 2) = dy(i + 2)
dy(i + 2) = dtemp
50 continue
return
end