Date: Tue, 5 Dec 1995 12:46:55 -0500 From: Netlib Reply Daemon To: S.Morgan@m.cc.utah.edu To: S.Morgan@m.cc.utah.edu Subject: Re: send dqrdc from linpack 1 #!/bin/sh # Caveat receptor. Jack Dongarra - dongarra@cs.utk.edu, # Eric Grosse - ehg@research.att.com # April 8, 1995 # # # Netlib is also accessible via the World Wide Web. # Our URL is http://www.netlib.org/index.html # # If you have trouble accessing netlib files please # send mail to netlib_maintainers@netlib.org # # In addition, anonymous ftp can be used to retrieve netlib files. # Ftp to ftp.netlib.org and log in as anonymous. Use your # complete e-mail address as the password. # # For a request returned in multiple mail messages, remove the # mail header and "sh" resulting files in any order. # to unpack, sh this message in an empty directory PATH=/bin:/usr/bin test ! -d blas && mkdir blas cat > blas/daxpy.f <<'bigmail CUT HERE............' 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 modified 12/3/93, array(1) declarations changed to array(*) c double precision dx(*),dy(*),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 bigmail CUT HERE............ test ! -d blas && mkdir blas cat > blas/ddot.f <<'bigmail CUT HERE............' 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 modified 12/3/93, array(1) declarations changed to array(*) c double precision dx(*),dy(*),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 bigmail CUT HERE............ test ! -d blas && mkdir blas cat > blas/dnrm2.f <<'bigmail CUT HERE............' DOUBLE PRECISION FUNCTION DNRM2 ( N, X, INCX ) * .. Scalar Arguments .. INTEGER INCX, N * .. Array Arguments .. DOUBLE PRECISION X( * ) * .. * * DNRM2 returns the euclidean norm of a vector via the function * name, so that * * DNRM2 := sqrt( x'*x ) * * * * -- This version written on 25-October-1982. * Modified on 14-October-1993 to inline the call to DLASSQ. * Sven Hammarling, Nag Ltd. * * * .. Parameters .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. Local Scalars .. INTEGER IX DOUBLE PRECISION ABSXI, NORM, SCALE, SSQ * .. Intrinsic Functions .. INTRINSIC ABS, SQRT * .. * .. Executable Statements .. IF( N.LT.1 .OR. INCX.LT.1 )THEN NORM = ZERO ELSE IF( N.EQ.1 )THEN NORM = ABS( X( 1 ) ) ELSE SCALE = ZERO SSQ = ONE * The following loop is equivalent to this call to the LAPACK * auxiliary routine: * CALL DLASSQ( N, X, INCX, SCALE, SSQ ) * DO 10, IX = 1, 1 + ( N - 1 )*INCX, INCX IF( X( IX ).NE.ZERO )THEN ABSXI = ABS( X( IX ) ) IF( SCALE.LT.ABSXI )THEN SSQ = ONE + SSQ*( SCALE/ABSXI )**2 SCALE = ABSXI ELSE SSQ = SSQ + ( ABSXI/SCALE )**2 END IF END IF 10 CONTINUE NORM = SCALE * SQRT( SSQ ) END IF * DNRM2 = NORM RETURN * * End of DNRM2. * END bigmail CUT HERE............ test ! -d blas && mkdir blas cat > blas/dscal.f <<'bigmail CUT HERE............' 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 modified 12/3/93, array(1) declarations changed to array(*) c double precision da,dx(*) 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 bigmail CUT HERE............ test ! -d blas && mkdir blas cat > blas/dswap.f <<'bigmail CUT HERE............' 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 modified 12/3/93, array(1) declarations changed to array(*) c double precision dx(*),dy(*),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 bigmail CUT HERE............ test ! -d linpack && mkdir linpack cat > linpack/dqrdc.f <<'bigmail CUT HERE............' 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 bigmail CUT HERE............ #define END