diff options
Diffstat (limited to 'math/llsq/qrbd.f')
| -rw-r--r-- | math/llsq/qrbd.f | 208 |
1 files changed, 208 insertions, 0 deletions
diff --git a/math/llsq/qrbd.f b/math/llsq/qrbd.f new file mode 100644 index 00000000..d9a12a24 --- /dev/null +++ b/math/llsq/qrbd.f @@ -0,0 +1,208 @@ +c subroutine qrbd (ipass,q,e,nn,v,mdv,nrv,c,mdc,ncc) +c c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 12 +c to appear in 'solving least squares problems', prentice-hall, 1974 +c qr algorithm for singular values of a bidiagonal matrix. +c +c the bidiagonal matrix +c +c (q1,e2,0... ) +c ( q2,e3,0... ) +c d= ( . ) +c ( . 0) +c ( .en) +c ( 0,qn) +c +c is pre and post multiplied by +c elementary rotation matrices +c ri and pi so that +c +c rk...r1*d*p1**(t)...pk**(t) = diag(s1,...,sn) +c +c to within working accuracy. +c +c 1. ei and qi occupy e(i) and q(i) as input. +c +c 2. rm...r1*c replaces 'c' in storage as output. +c +c 3. v*p1**(t)...pm**(t) replaces 'v' in storage as output. +c +c 4. si occupies q(i) as output. +c +c 5. the si's are nonincreasing and nonnegative. +c +c this code is based on the paper and 'algol' code.. +c ref.. +c 1. reinsch,c.h. and golub,g.h. 'singular value decomposition +c and least squares solutions' (numer. math.), vol. 14,(1970). +c + subroutine qrbd (ipass,q,e,nn,v,mdv,nrv,c,mdc,ncc) + logical wntv ,havers,fail + dimension q(nn),e(nn),v(mdv,nn),c(mdc,ncc) + zero=0. + one=1. + two=2. +c + n=nn + ipass=1 + if (n.le.0) return + n10=10*n + wntv=nrv.gt.0 + havers=ncc.gt.0 + fail=.false. + nqrs=0 + e(1)=zero + dnorm=zero + do 10 j=1,n + 10 dnorm=amax1(abs(q(j))+abs(e(j)),dnorm) + do 200 kk=1,n + k=n+1-kk +c +c test for splitting or rank deficiencies.. +c first make test for last diagonal term, q(k), being small. + 20 if(k.eq.1) go to 50 + if(diff(dnorm+q(k),dnorm)) 50,25,50 +c +c since q(k) is small we will make a special pass to +c transform e(k) to zero. +c + 25 cs=zero + sn=-one + do 40 ii=2,k + i=k+1-ii + f=-sn*e(i+1) + e(i+1)=cs*e(i+1) + call g1 (q(i),f,cs,sn,q(i)) +c transformation constructed to zero position (i,k). +c + if (.not.wntv) go to 40 + do 30 j=1,nrv + 30 call g2 (cs,sn,v(j,i),v(j,k)) +c accumulate rt. transformations in v. +c + 40 continue +c +c the matrix is now bidiagonal, and of lower order +c since e(k) .eq. zero.. +c + 50 do 60 ll=1,k + l=k+1-ll + if(diff(dnorm+e(l),dnorm)) 55,100,55 + 55 if(diff(dnorm+q(l-1),dnorm)) 60,70,60 + 60 continue +c this loop can't complete since e(1) = zero. +c + go to 100 +c +c cancellation of e(l), l.gt.1. + 70 cs=zero + sn=-one + do 90 i=l,k + f=-sn*e(i) + e(i)=cs*e(i) + if(diff(dnorm+f,dnorm)) 75,100,75 + 75 call g1 (q(i),f,cs,sn,q(i)) + if (.not.havers) go to 90 + do 80 j=1,ncc + 80 call g2 (cs,sn,c(i,j),c(l-1,j)) + 90 continue +c +c test for convergence.. + 100 z=q(k) + if (l.eq.k) go to 170 +c +c shift from bottom 2 by 2 minor of b**(t)*b. + x=q(l) + y=q(k-1) + g=e(k-1) + h=e(k) + f=((y-z)*(y+z)+(g-h)*(g+h))/(two*h*y) + g=sqrt(one+f**2) + if (f.lt.zero) go to 110 + t=f+g + go to 120 + 110 t=f-g + 120 f=((x-z)*(x+z)+h*(y/t-h))/x +c +c next qr sweep.. + cs=one + sn=one + lp1=l+1 + do 160 i=lp1,k + g=e(i) + y=q(i) + h=sn*g + g=cs*g + call g1 (f,h,cs,sn,e(i-1)) + f=x*cs+g*sn + g=-x*sn+g*cs + h=y*sn + y=y*cs + if (.not.wntv) go to 140 +c +c accumulate rotations (from the right) in 'v' + do 130 j=1,nrv + 130 call g2 (cs,sn,v(j,i-1),v(j,i)) + 140 call g1 (f,h,cs,sn,q(i-1)) + f=cs*g+sn*y + x=-sn*g+cs*y + if (.not.havers) go to 160 + do 150 j=1,ncc + 150 call g2 (cs,sn,c(i-1,j),c(i,j)) +c apply rotations from the left to +c right hand sides in 'c'.. +c + 160 continue + e(l)=zero + e(k)=f + q(k)=x + nqrs=nqrs+1 + if (nqrs.le.n10) go to 20 +c return to 'test for splitting'. +c + fail=.true. +c .. +c cutoff for convergence failure. 'nqrs' will be 2*n usually. + 170 if (z.ge.zero) go to 190 + q(k)=-z + if (.not.wntv) go to 190 + do 180 j=1,nrv + 180 v(j,k)=-v(j,k) + 190 continue +c convergence. q(k) is made nonnegative.. +c + 200 continue + if (n.eq.1) return + do 210 i=2,n + if (q(i).gt.q(i-1)) go to 220 + 210 continue + if (fail) ipass=2 + return +c .. +c every singular value is in order.. + 220 do 270 i=2,n + t=q(i-1) + k=i-1 + do 230 j=i,n + if (t.ge.q(j)) go to 230 + t=q(j) + k=j + 230 continue + if (k.eq.i-1) go to 270 + q(k)=q(i-1) + q(i-1)=t + if (.not.havers) go to 250 + do 240 j=1,ncc + t=c(i-1,j) + c(i-1,j)=c(k,j) + 240 c(k,j)=t + 250 if (.not.wntv) go to 270 + do 260 j=1,nrv + t=v(j,i-1) + v(j,i-1)=v(j,k) + 260 v(j,k)=t + 270 continue +c end of ordering algorithm. +c + if (fail) ipass=2 + return + end |