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);
619 if( idir==1 && amp::abs<Precision>(d(ll))<
amp::ampf<Precision>(
"1.0E-3")*amp::abs<Precision>(d(m)) )
623 if( idir==2 && amp::abs<Precision>(d(m))<
amp::ampf<Precision>(
"1.0E-3")*amp::abs<Precision>(d(ll)) )
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)