314 bool matrixsplitflag;
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);
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++)
static const ampf getAlgoPascalMinNumber()
void mu(int **points, int sizePoints)
void vsub(raw_vector< T > vdst, const_raw_vector< T > vsrc)
static const ampf getAlgoPascalEpsilon()
void vmul(raw_vector< T > vdst, T2 alpha)
void vadd(raw_vector< T > vdst, const_raw_vector< T > vsrc)
raw_vector< T > getrow(int iRow, int iColumnStart, int iColumnEnd)
raw_vector< T > getcolumn(int iColumn, int iRowStart, int iRowEnd)
int maxint(int m1, int m2)
void setbounds(int iLow, int iHigh)
raw_vector< T > getvector(int iStart, int iEnd)
void vmove(raw_vector< T > vdst, const_raw_vector< T > vsrc)