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