352 work0.setbounds(1, n-1);
353 work1.setbounds(1, n-1);
354 work2.setbounds(1, n-1);
355 work3.setbounds(1, n-1);
359 utemp.setbounds(ustart, uend);
360 vttemp.setbounds(vstart, vend);
361 ctemp.setbounds(cstart, cend);
369 etemp.setbounds(1, n);
370 for(
i=1;
i<=n-1;
i++)
375 for(
i=1;
i<=n-1;
i++)
394 for(
i=1;
i<=n-1;
i++)
396 rotations::generaterotation<Precision>(d(
i), e(
i), cs, sn, r);
409 rotations::applyrotationsfromtheright<Precision>(fwddir, ustart, uend, 1+ustart-1, n+ustart-1, work0, work1, u, utemp);
413 rotations::applyrotationsfromtheleft<Precision>(fwddir, 1+cstart-1, n+cstart-1, cstart, cend, work0, work1, c, ctemp);
422 tolmul = amp::maximum<Precision>(10, amp::minimum<Precision>(100, amp::pow<Precision>(eps, -
amp::ampf<Precision>(
"0.125"))));
424 if( !isfractionalaccuracyrequired )
435 smax = amp::maximum<Precision>(smax, amp::abs<Precision>(d(
i)));
437 for(
i=1;
i<=n-1;
i++)
439 smax = amp::maximum<Precision>(smax, amp::abs<Precision>(e(
i)));
448 sminoa = amp::abs<Precision>(d(1));
454 mu = amp::abs<Precision>(d(
i))*(
mu/(
mu+amp::abs<Precision>(e(
i-1))));
455 sminoa = amp::minimum<Precision>(sminoa,
mu);
462 sminoa = sminoa/amp::sqrt<Precision>(n);
463 thresh = amp::maximum<Precision>(tol*sminoa, maxitr*n*n*unfl);
471 thresh = amp::maximum<Precision>(amp::abs<Precision>(tol)*smax, maxitr*n*n*unfl);
511 if( tol<0 && amp::abs<Precision>(d(
m))<=thresh )
515 smax = amp::abs<Precision>(d(
m));
517 matrixsplitflag =
false;
518 for(lll=1; lll<=
m-1; lll++)
521 abss = amp::abs<Precision>(d(ll));
522 abse = amp::abs<Precision>(e(ll));
523 if( tol<0 && abss<=thresh )
529 matrixsplitflag =
true;
532 smin = amp::minimum<Precision>(smin, abss);
533 smax = amp::maximum<Precision>(smax, amp::maximum<Precision>(abss, abse));
535 if( !matrixsplitflag )
567 svdv2x2<Precision>(d(
m-1), e(
m-1), d(
m), sigmn, sigmx, sinr, cosr, sinl, cosl);
578 mm1 =
m-1+(vstart-1);
579 ap::vmove(vttemp.getvector(vstart, vend), vt.
getrow(mm1, vstart, vend), cosr);
580 ap::vadd(vttemp.getvector(vstart, vend), vt.
getrow(mm0, vstart, vend), sinr);
583 ap::vmove(vt.
getrow(mm1, vstart, vend), vttemp.getvector(vstart, vend));
599 ap::vmove(ctemp.getvector(cstart, cend), c.
getrow(mm1, cstart, cend), cosl);
600 ap::vadd(ctemp.getvector(cstart, cend), c.
getrow(mm0, cstart, cend), sinl);
627 if( ll!=oldll ||
m!=oldm || bchangedir )
629 if( amp::abs<Precision>(d(ll))>=amp::abs<Precision>(d(
m)) )
657 if( amp::abs<Precision>(e(
m-1))<=amp::abs<Precision>(tol)*amp::abs<Precision>(d(
m)) || tol<0 && amp::abs<Precision>(e(
m-1))<=thresh )
669 mu = amp::abs<Precision>(d(ll));
672 for(lll=ll; lll<=
m-1; lll++)
674 if( amp::abs<Precision>(e(lll))<=tol*
mu )
681 mu = amp::abs<Precision>(d(lll+1))*(
mu/(
mu+amp::abs<Precision>(e(lll))));
682 sminl = amp::minimum<Precision>(sminl,
mu);
697 if( amp::abs<Precision>(e(ll))<=amp::abs<Precision>(tol)*amp::abs<Precision>(d(ll)) || tol<0 && amp::abs<Precision>(e(ll))<=thresh )
709 mu = amp::abs<Precision>(d(
m));
712 for(lll=
m-1; lll>=ll; lll--)
714 if( amp::abs<Precision>(e(lll))<=tol*
mu )
721 mu = amp::abs<Precision>(d(lll))*(
mu/(
mu+amp::abs<Precision>(e(lll))));
722 sminl = amp::minimum<Precision>(sminl,
mu);
737 if( tol>=0 && n*tol*(sminl/smax)<=amp::maximum<Precision>(eps,
amp::ampf<Precision>(
"0.01")*tol) )
753 sll = amp::abs<Precision>(d(ll));
754 svd2x2<Precision>(d(
m-1), e(
m-1), d(
m), shift, r);
758 sll = amp::abs<Precision>(d(
m));
759 svd2x2<Precision>(d(ll), e(ll), d(ll+1), shift, r);
767 if( amp::sqr<Precision>(shift/sll)<eps )
793 for(
i=ll;
i<=
m-1;
i++)
795 rotations::generaterotation<Precision>(d(
i)*cs, e(
i), cs, sn, r);
800 rotations::generaterotation<Precision>(oldcs*r, d(
i+1)*sn, oldcs, oldsn, tmp);
804 work2(
i-ll+1) = oldcs;
805 work3(
i-ll+1) = oldsn;
816 rotations::applyrotationsfromtheleft<Precision>(fwddir, ll+vstart-1,
m+vstart-1, vstart, vend, work0, work1, vt, vttemp);
820 rotations::applyrotationsfromtheright<Precision>(fwddir, ustart, uend, ll+ustart-1,
m+ustart-1, work2, work3, u, utemp);
824 rotations::applyrotationsfromtheleft<Precision>(fwddir, ll+cstart-1,
m+cstart-1, cstart, cend, work2, work3, c, ctemp);
830 if( amp::abs<Precision>(e(
m-1))<=thresh )
844 for(
i=
m;
i>=ll+1;
i--)
846 rotations::generaterotation<Precision>(d(
i)*cs, e(
i-1), cs, sn, r);
851 rotations::generaterotation<Precision>(oldcs*r, d(
i-1)*sn, oldcs, oldsn, tmp);
856 work3(
i-ll) = -oldsn;
867 rotations::applyrotationsfromtheleft<Precision>(!fwddir, ll+vstart-1,
m+vstart-1, vstart, vend, work2, work3, vt, vttemp);
871 rotations::applyrotationsfromtheright<Precision>(!fwddir, ustart, uend, ll+ustart-1,
m+ustart-1, work0, work1, u, utemp);
875 rotations::applyrotationsfromtheleft<Precision>(!fwddir, ll+cstart-1,
m+cstart-1, cstart, cend, work0, work1, c, ctemp);
881 if( amp::abs<Precision>(e(ll))<=thresh )
900 f = (amp::abs<Precision>(d(ll))-shift)*(extsignbdsqr<Precision>(1, d(ll))+shift/d(ll));
902 for(
i=ll;
i<=
m-1;
i++)
904 rotations::generaterotation<Precision>(
f,
g, cosr, sinr, r);
909 f = cosr*d(
i)+sinr*e(
i);
910 e(
i) = cosr*e(
i)-sinr*d(
i);
912 d(
i+1) = cosr*d(
i+1);
913 rotations::generaterotation<Precision>(
f,
g, cosl, sinl, r);
915 f = cosl*e(
i)+sinl*d(
i+1);
916 d(
i+1) = cosl*d(
i+1)-sinl*e(
i);
920 e(
i+1) = cosl*e(
i+1);
922 work0(
i-ll+1) = cosr;
923 work1(
i-ll+1) = sinr;
924 work2(
i-ll+1) = cosl;
925 work3(
i-ll+1) = sinl;
934 rotations::applyrotationsfromtheleft<Precision>(fwddir, ll+vstart-1,
m+vstart-1, vstart, vend, work0, work1, vt, vttemp);
938 rotations::applyrotationsfromtheright<Precision>(fwddir, ustart, uend, ll+ustart-1,
m+ustart-1, work2, work3, u, utemp);
942 rotations::applyrotationsfromtheleft<Precision>(fwddir, ll+cstart-1,
m+cstart-1, cstart, cend, work2, work3, c, ctemp);
948 if( amp::abs<Precision>(e(
m-1))<=thresh )
960 f = (amp::abs<Precision>(d(
m))-shift)*(extsignbdsqr<Precision>(1, d(
m))+shift/d(
m));
962 for(
i=
m;
i>=ll+1;
i--)
964 rotations::generaterotation<Precision>(
f,
g, cosr, sinr, r);
969 f = cosr*d(
i)+sinr*e(
i-1);
970 e(
i-1) = cosr*e(
i-1)-sinr*d(
i);
972 d(
i-1) = cosr*d(
i-1);
973 rotations::generaterotation<Precision>(
f,
g, cosl, sinl, r);
975 f = cosl*e(
i-1)+sinl*d(
i-1);
976 d(
i-1) = cosl*d(
i-1)-sinl*e(
i-1);
980 e(
i-2) = cosl*e(
i-2);
992 if( amp::abs<Precision>(e(ll))<=thresh )
1002 rotations::applyrotationsfromtheleft<Precision>(!fwddir, ll+vstart-1,
m+vstart-1, vstart, vend, work2, work3, vt, vttemp);
1006 rotations::applyrotationsfromtheright<Precision>(!fwddir, ustart, uend, ll+ustart-1,
m+ustart-1, work0, work1, u, utemp);
1010 rotations::applyrotationsfromtheleft<Precision>(!fwddir, ll+cstart-1,
m+cstart-1, cstart, cend, work0, work1, c, ctemp);
1044 for(
i=1;
i<=n-1;
i++)
1052 for(
j=2;
j<=n+1-
i;
j++)
1071 ap::vmove(vttemp.getvector(vstart, vend), vt.
getrow(isub+vstart-1, vstart, vend));
1073 ap::vmove(vt.
getrow(
j+vstart-1, vstart, vend), vttemp.getvector(vstart, vend));
1085 ap::vmove(ctemp.getvector(cstart, cend), c.
getrow(isub+cstart-1, cstart, cend));
1087 ap::vmove(c.
getrow(
j+cstart-1, cstart, cend), ctemp.getvector(cstart, cend));
1095 template<
unsigned int Precision>
1104 result = amp::abs<Precision>(a);
1108 result = -amp::abs<Precision>(a);
1114 template<
unsigned int Precision>