Math-LOESS
view release on metacpan or search on metacpan
loess/linpack_lite.f view on Meta::CPAN
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)
loess/linpack_lite.f view on Meta::CPAN
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
( run in 0.736 second using v1.01-cache-2.11-cpan-71847e10f99 )