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
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
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
subroutine dpodi(a,lda,n,det,job)
integer lda,n,job
double precision a(lda,1)
double precision det(2)
c
c dpodi computes the determinant and inverse of a certain
c double precision symmetric positive definite matrix (see below)
c using the factors computed by dpoco, dpofa or dqrdc.
c
c on entry
c
c a double precision(lda, n)
c the output a from dpoco or dpofa
c or the output x from dqrdc.
c
c lda integer
c the leading dimension of the array a .
c
c n integer
c the order of the matrix a .
c
c job integer
c = 11 both determinant and inverse.
c = 01 inverse only.
c = 10 determinant only.
c
c on return
c
c a if dpoco or dpofa was used to factor a then
c dpodi produces the upper half of inverse(a) .
c if dqrdc was used to decompose x then
c dpodi produces the upper half of inverse(trans(x)*x)
c where trans(x) is the transpose.
c elements of a below the diagonal are unchanged.
c if the units digit of job is zero, a is unchanged.
c
c det double precision(2)
c determinant of a or of trans(x)*x if requested.
c otherwise not referenced.
c determinant = det(1) * 10.0**det(2)
c with 1.0 .le. det(1) .lt. 10.0
c or det(1) .eq. 0.0 .
c
c error condition
c
c a division by zero will occur if the input factor contains
c a zero on the diagonal and the inverse is requested.
c it will not occur if the subroutines are called correctly
c and if dpoco or dpofa has set info .eq. 0 .
c
c linpack. this version dated 08/14/78 .
c cleve moler, university of new mexico, argonne national lab.
c
c subroutines and functions
c
c blas daxpy,dscal
c fortran mod
c
c internal variables
c
double precision t
double precision s
integer i,j,jm1,k,kp1
c
c compute determinant
c
if (job/10 .eq. 0) go to 70
det(1) = 1.0d0
det(2) = 0.0d0
s = 10.0d0
do 50 i = 1, n
det(1) = a(i,i)**2*det(1)
c ...exit
if (det(1) .eq. 0.0d0) go to 60
10 if (det(1) .ge. 1.0d0) go to 20
det(1) = s*det(1)
det(2) = det(2) - 1.0d0
go to 10
20 continue
30 if (det(1) .lt. s) go to 40
det(1) = det(1)/s
det(2) = det(2) + 1.0d0
go to 30
40 continue
50 continue
60 continue
70 continue
c
c compute inverse(r)
c
if (mod(job,10) .eq. 0) go to 140
do 100 k = 1, n
a(k,k) = 1.0d0/a(k,k)
t = -a(k,k)
call dscal(k-1,t,a(1,k),1)
kp1 = k + 1
if (n .lt. kp1) go to 90
do 80 j = kp1, n
t = a(k,j)
a(k,j) = 0.0d0
call daxpy(k,t,a(1,k),1,a(1,j),1)
80 continue
90 continue
100 continue
c
c form inverse(r) * trans(inverse(r))
c
do 130 j = 1, n
jm1 = j - 1
if (jm1 .lt. 1) go to 120
do 110 k = 1, jm1
t = a(k,j)
call daxpy(k,t,a(1,j),1,a(1,k),1)
110 continue
120 continue
t = a(j,j)
call dscal(j,t,a(1,j),1)
130 continue
140 continue
return
end
subroutine dpofa(a,lda,n,info)
integer lda,n,info
double precision a(lda,1)
c
c dpofa factors a double precision symmetric positive definite
c matrix.
c
c dpofa is usually called by dpoco, but it can be called
c directly with a saving in time if rcond is not needed.
c (time for dpoco) = (1 + 18/n)*(time for dpofa) .
c
c on entry
c
c a double precision(lda, n)
c the symmetric matrix to be factored. only the
c diagonal and upper triangle are used.
c
c lda integer
c the leading dimension of the array a .
c
c n integer
c the order of the matrix a .
c
c on return
c
c a an upper triangular matrix r so that a = trans(r)*r
c where trans(r) is the transpose.
c the strict lower triangle is unaltered.
c if info .ne. 0 , the factorization is not complete.
c
c info integer
c = 0 for normal return.
c = k signals an error condition. the leading minor
c of order k is not positive definite.
c
c linpack. this version dated 08/14/78 .
c cleve moler, university of new mexico, argonne national lab.
c
c subroutines and functions
c
c blas ddot
c fortran dsqrt
c
c internal variables
c
double precision ddot,t
double precision s
integer j,jm1,k
c begin block with ...exits to 40
c
c
do 30 j = 1, n
info = j
s = 0.0d0
jm1 = j - 1
if (jm1 .lt. 1) go to 20
do 10 k = 1, jm1
t = a(k,j) - ddot(k-1,a(1,k),1,a(1,j),1)
t = t/a(k,k)
a(k,j) = t
s = s + t*t
10 continue
20 continue
s = a(j,j) - s
c ......exit
if (s .le. 0.0d0) go to 40
a(j,j) = dsqrt(s)
30 continue
info = 0
40 continue
return
end
subroutine drot (n,dx,incx,dy,incy,c,s)
c
c applies a plane rotation.
c jack dongarra, linpack, 3/11/78.
c modified 12/3/93, array(1) declarations changed to array(*)
c
double precision dx(*),dy(*),dtemp,c,s
integer i,incx,incy,ix,iy,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 = c*dx(ix) + s*dy(iy)
dy(iy) = c*dy(iy) - s*dx(ix)
dx(ix) = dtemp
ix = ix + incx
iy = iy + incy
10 continue
return
c
c code for both increments equal to 1
c
20 do 30 i = 1,n
dtemp = c*dx(i) + s*dy(i)
dy(i) = c*dy(i) - s*dx(i)
dx(i) = dtemp
30 continue
return
end
subroutine drotg(da,db,c,s)
c
c construct givens plane rotation.
c jack dongarra, linpack, 3/11/78.
c
double precision da,db,c,s,roe,scale,r,z
c
roe = db
if( dabs(da) .gt. dabs(db) ) roe = da
scale = dabs(da) + dabs(db)
if( scale .ne. 0.0d0 ) go to 10
c = 1.0d0
s = 0.0d0
r = 0.0d0
z = 0.0d0
go to 20
10 r = scale*dsqrt((da/scale)**2 + (db/scale)**2)
r = dsign(1.0d0,roe)*r
c = da/r
s = db/r
z = 1.0d0
if( dabs(da) .gt. dabs(db) ) z = s
if( dabs(db) .ge. dabs(da) .and. c .ne. 0.0d0 ) z = 1.0d0/c
20 da = r
db = z
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 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
subroutine dsvdc(x,ldx,n,p,s,e,u,ldu,v,ldv,work,job,info)
integer ldx,n,p,ldu,ldv,job,info
double precision x(ldx,1),s(1),e(1),u(ldu,1),v(ldv,1),work(1)
c
c
c dsvdc is a subroutine to reduce a double precision nxp matrix x
c by orthogonal transformations u and v to diagonal form. the
c diagonal elements s(i) are the singular values of x. the
c columns of u are the corresponding left singular vectors,
c and the columns of v the right singular vectors.
c
c on entry
c
c x double precision(ldx,p), where ldx.ge.n.
c x contains the matrix whose singular value
c decomposition is to be computed. x is
c destroyed by dsvdc.
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 ldu integer.
c ldu is the leading dimension of the array u.
c (see below).
c
c ldv integer.
c ldv is the leading dimension of the array v.
c (see below).
c
c work double precision(n).
c work is a scratch array.
c
c job integer.
c job controls the computation of the singular
c vectors. it has the decimal expansion ab
c with the following meaning
c
c a.eq.0 do not compute the left singular
c vectors.
c a.eq.1 return the n left singular vectors
c in u.
c a.ge.2 return the first min(n,p) singular
c vectors in u.
c b.eq.0 do not compute the right singular
c vectors.
c b.eq.1 return the right singular vectors
c in v.
c
c on return
c
c s double precision(mm), where mm=min(n+1,p).
c the first min(n,p) entries of s contain the
c singular values of x arranged in descending
c order of magnitude.
c
c e double precision(p),
c e ordinarily contains zeros. however see the
c discussion of info for exceptions.
c
c u double precision(ldu,k), where ldu.ge.n. if
c joba.eq.1 then k.eq.n, if joba.ge.2
c then k.eq.min(n,p).
c u contains the matrix of left singular vectors.
c u is not referenced if joba.eq.0. if n.le.p
c or if joba.eq.2, then u may be identified with x
c in the subroutine call.
c
c v double precision(ldv,p), where ldv.ge.p.
c v contains the matrix of right singular vectors.
c v is not referenced if job.eq.0. if p.le.n,
c then v may be identified with x in the
c subroutine call.
c
c info integer.
c the singular values (and their corresponding
c singular vectors) s(info+1),s(info+2),...,s(m)
c are correct (here m=min(n,p)). thus if
c info.eq.0, all the singular values and their
c vectors are correct. in any event, the matrix
c b = trans(u)*x*v is the bidiagonal matrix
c with the elements of s on its diagonal and the
c elements of e on its super-diagonal (trans(u)
c is the transpose of u). thus the singular
c values of x and b are the same.
c
c linpack. this version dated 08/14/78 .
c correction made to shift 2/84.
c g.w. stewart, university of maryland, argonne national lab.
c
c dsvdc uses the following functions and subprograms.
c
c external drot
c blas daxpy,ddot,dscal,dswap,dnrm2,drotg
c fortran dabs,dmax1,max0,min0,mod,dsqrt
c
c internal variables
c
integer i,iter,j,jobu,k,kase,kk,l,ll,lls,lm1,lp1,ls,lu,m,maxit,
* mm,mm1,mp1,nct,nctp1,ncu,nrt,nrtp1
double precision ddot,t,r
double precision b,c,cs,el,emm1,f,g,dnrm2,scale,shift,sl,sm,sn,
* smm1,t1,test,ztest
logical wantu,wantv
c
c
c set the maximum number of iterations.
c
maxit = 30
c
c determine what is to be computed.
c
wantu = .false.
wantv = .false.
jobu = mod(job,100)/10
ncu = n
if (jobu .gt. 1) ncu = min0(n,p)
if (jobu .ne. 0) wantu = .true.
if (mod(job,10) .ne. 0) wantv = .true.
c
c reduce x to bidiagonal form, storing the diagonal elements
c in s and the super-diagonal elements in e.
c
info = 0
nct = min0(n-1,p)
nrt = max0(0,min0(p-2,n))
lu = max0(nct,nrt)
if (lu .lt. 1) go to 170
do 160 l = 1, lu
lp1 = l + 1
if (l .gt. nct) go to 20
c
c compute the transformation for the l-th column and
c place the l-th diagonal in s(l).
c
s(l) = dnrm2(n-l+1,x(l,l),1)
if (s(l) .eq. 0.0d0) go to 10
if (x(l,l) .ne. 0.0d0) s(l) = dsign(s(l),x(l,l))
call dscal(n-l+1,1.0d0/s(l),x(l,l),1)
x(l,l) = 1.0d0 + x(l,l)
10 continue
s(l) = -s(l)
20 continue
if (p .lt. lp1) go to 50
do 40 j = lp1, p
if (l .gt. nct) go to 30
if (s(l) .eq. 0.0d0) go to 30
c
c apply the transformation.
c
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)
30 continue
c
c place the l-th row of x into e for the
c subsequent calculation of the row transformation.
c
e(j) = x(l,j)
40 continue
50 continue
if (.not.wantu .or. l .gt. nct) go to 70
c
c place the transformation in u for subsequent back
c multiplication.
c
do 60 i = l, n
u(i,l) = x(i,l)
60 continue
70 continue
if (l .gt. nrt) go to 150
c
c compute the l-th row transformation and place the
c l-th super-diagonal in e(l).
c
e(l) = dnrm2(p-l,e(lp1),1)
if (e(l) .eq. 0.0d0) go to 80
if (e(lp1) .ne. 0.0d0) e(l) = dsign(e(l),e(lp1))
call dscal(p-l,1.0d0/e(l),e(lp1),1)
e(lp1) = 1.0d0 + e(lp1)
80 continue
e(l) = -e(l)
if (lp1 .gt. n .or. e(l) .eq. 0.0d0) go to 120
c
c apply the transformation.
c
do 90 i = lp1, n
work(i) = 0.0d0
90 continue
do 100 j = lp1, p
call daxpy(n-l,e(j),x(lp1,j),1,work(lp1),1)
100 continue
do 110 j = lp1, p
call daxpy(n-l,-e(j)/e(lp1),work(lp1),1,x(lp1,j),1)
110 continue
120 continue
if (.not.wantv) go to 140
c
c place the transformation in v for subsequent
c back multiplication.
c
do 130 i = lp1, p
v(i,l) = e(i)
130 continue
140 continue
150 continue
160 continue
170 continue
c
c set up the final bidiagonal matrix or order m.
c
m = min0(p,n+1)
nctp1 = nct + 1
nrtp1 = nrt + 1
if (nct .lt. p) s(nctp1) = x(nctp1,nctp1)
if (n .lt. m) s(m) = 0.0d0
if (nrtp1 .lt. m) e(nrtp1) = x(nrtp1,m)
e(m) = 0.0d0
c
c if required, generate u.
c
if (.not.wantu) go to 300
if (ncu .lt. nctp1) go to 200
do 190 j = nctp1, ncu
do 180 i = 1, n
u(i,j) = 0.0d0
180 continue
u(j,j) = 1.0d0
190 continue
200 continue
if (nct .lt. 1) go to 290
do 280 ll = 1, nct
l = nct - ll + 1
if (s(l) .eq. 0.0d0) go to 250
lp1 = l + 1
if (ncu .lt. lp1) go to 220
do 210 j = lp1, ncu
t = -ddot(n-l+1,u(l,l),1,u(l,j),1)/u(l,l)
call daxpy(n-l+1,t,u(l,l),1,u(l,j),1)
210 continue
220 continue
call dscal(n-l+1,-1.0d0,u(l,l),1)
u(l,l) = 1.0d0 + u(l,l)
lm1 = l - 1
if (lm1 .lt. 1) go to 240
do 230 i = 1, lm1
u(i,l) = 0.0d0
230 continue
240 continue
go to 270
250 continue
do 260 i = 1, n
u(i,l) = 0.0d0
260 continue
u(l,l) = 1.0d0
270 continue
280 continue
290 continue
300 continue
c
c if it is required, generate v.
c
if (.not.wantv) go to 350
do 340 ll = 1, p
l = p - ll + 1
lp1 = l + 1
if (l .gt. nrt) go to 320
if (e(l) .eq. 0.0d0) go to 320
do 310 j = lp1, p
t = -ddot(p-l,v(lp1,l),1,v(lp1,j),1)/v(lp1,l)
call daxpy(p-l,t,v(lp1,l),1,v(lp1,j),1)
310 continue
320 continue
do 330 i = 1, p
v(i,l) = 0.0d0
330 continue
v(l,l) = 1.0d0
340 continue
350 continue
c
c main iteration loop for the singular values.
c
mm = m
iter = 0
360 continue
c
c quit if all the singular values have been found.
c
c ...exit
if (m .eq. 0) go to 620
c
c if too many iterations have been performed, set
c flag and return.
c
if (iter .lt. maxit) go to 370
info = m
c ......exit
go to 620
370 continue
c
c this section of the program inspects for
c negligible elements in the s and e arrays. on
c completion the variables kase and l are set as follows.
c
c kase = 1 if s(m) and e(l-1) are negligible and l.lt.m
c kase = 2 if s(l) is negligible and l.lt.m
c kase = 3 if e(l-1) is negligible, l.lt.m, and
c s(l), ..., s(m) are not negligible (qr step).
c kase = 4 if e(m-1) is negligible (convergence).
c
do 390 ll = 1, m
l = m - ll
c ...exit
if (l .eq. 0) go to 400
test = dabs(s(l)) + dabs(s(l+1))
ztest = test + dabs(e(l))
if (ztest .ne. test) go to 380
e(l) = 0.0d0
c ......exit
go to 400
380 continue
390 continue
400 continue
if (l .ne. m - 1) go to 410
kase = 4
go to 480
410 continue
lp1 = l + 1
mp1 = m + 1
do 430 lls = lp1, mp1
ls = m - lls + lp1
c ...exit
if (ls .eq. l) go to 440
test = 0.0d0
if (ls .ne. m) test = test + dabs(e(ls))
if (ls .ne. l + 1) test = test + dabs(e(ls-1))
ztest = test + dabs(s(ls))
if (ztest .ne. test) go to 420
s(ls) = 0.0d0
c ......exit
go to 440
420 continue
430 continue
440 continue
if (ls .ne. l) go to 450
kase = 3
go to 470
450 continue
if (ls .ne. m) go to 460
kase = 1
go to 470
460 continue
kase = 2
l = ls
470 continue
480 continue
l = l + 1
c
c perform the task indicated by kase.
c
go to (490,520,540,570), kase
c
c deflate negligible s(m).
c
490 continue
mm1 = m - 1
f = e(m-1)
e(m-1) = 0.0d0
do 510 kk = l, mm1
k = mm1 - kk + l
t1 = s(k)
call drotg(t1,f,cs,sn)
s(k) = t1
if (k .eq. l) go to 500
f = -sn*e(k-1)
e(k-1) = cs*e(k-1)
500 continue
if (wantv) call drot(p,v(1,k),1,v(1,m),1,cs,sn)
510 continue
go to 610
c
c split at negligible s(l).
c
520 continue
f = e(l-1)
e(l-1) = 0.0d0
do 530 k = l, m
t1 = s(k)
call drotg(t1,f,cs,sn)
s(k) = t1
f = -sn*e(k)
e(k) = cs*e(k)
if (wantu) call drot(n,u(1,k),1,u(1,l-1),1,cs,sn)
530 continue
go to 610
c
c perform one qr step.
c
540 continue
c
c calculate the shift.
c
scale = dmax1(dabs(s(m)),dabs(s(m-1)),dabs(e(m-1)),
* dabs(s(l)),dabs(e(l)))
sm = s(m)/scale
smm1 = s(m-1)/scale
emm1 = e(m-1)/scale
sl = s(l)/scale
el = e(l)/scale
b = ((smm1 + sm)*(smm1 - sm) + emm1**2)/2.0d0
c = (sm*emm1)**2
shift = 0.0d0
if (b .eq. 0.0d0 .and. c .eq. 0.0d0) go to 550
shift = dsqrt(b**2+c)
if (b .lt. 0.0d0) shift = -shift
shift = c/(b + shift)
550 continue
f = (sl + sm)*(sl - sm) + shift
g = sl*el
c
c chase zeros.
c
mm1 = m - 1
do 560 k = l, mm1
call drotg(f,g,cs,sn)
if (k .ne. l) e(k-1) = f
f = cs*s(k) + sn*e(k)
e(k) = cs*e(k) - sn*s(k)
g = sn*s(k+1)
s(k+1) = cs*s(k+1)
if (wantv) call drot(p,v(1,k),1,v(1,k+1),1,cs,sn)
call drotg(f,g,cs,sn)
s(k) = f
f = cs*e(k) + sn*s(k+1)
s(k+1) = -sn*e(k) + cs*s(k+1)
g = sn*e(k+1)
e(k+1) = cs*e(k+1)
if (wantu .and. k .lt. n)
* call drot(n,u(1,k),1,u(1,k+1),1,cs,sn)
560 continue
e(m-1) = f
iter = iter + 1
go to 610
c
c convergence.
c
570 continue
c
c make the singular value positive.
c
if (s(l) .ge. 0.0d0) go to 580
s(l) = -s(l)
if (wantv) call dscal(p,-1.0d0,v(1,l),1)
580 continue
c
c order the singular value.
c
590 if (l .eq. mm) go to 600
c ...exit
if (s(l) .ge. s(l+1)) go to 600
t = s(l)
s(l) = s(l+1)
s(l+1) = t
if (wantv .and. l .lt. p)
* call dswap(p,v(1,l),1,v(1,l+1),1)
if (wantu .and. l .lt. n)
* call dswap(n,u(1,l),1,u(1,l+1),1)
l = l + 1
go to 590
600 continue
iter = 0
m = m - 1
610 continue
go to 360
620 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 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