Loading...
Searching...
No Matches
bdsvd Namespace Reference

## Functions

template<unsigned int Precision>
bool | rmatrixbdsvd (ap::template_1d_array< amp::ampf< Precision > > &d, ap::template_1d_array< amp::ampf< Precision > > e, int n, bool isupper, bool isfractionalaccuracyrequired, ap::template_2d_array< amp::ampf< Precision > > &u, int nru, ap::template_2d_array< amp::ampf< Precision > > &c, int ncc, ap::template_2d_array< amp::ampf< Precision > > &vt, int ncvt)
template<unsigned int Precision>
bool | bidiagonalsvddecomposition (ap::template_1d_array< amp::ampf< Precision > > &d, ap::template_1d_array< amp::ampf< Precision > > e, int n, bool isupper, bool isfractionalaccuracyrequired, ap::template_2d_array< amp::ampf< Precision > > &u, int nru, ap::template_2d_array< amp::ampf< Precision > > &c, int ncc, ap::template_2d_array< amp::ampf< Precision > > &vt, int ncvt)
template<unsigned int Precision>

bool

bidiagonalsvddecompositioninternal (ap::template_1d_array< amp::ampf< Precision > > &d, ap::template_1d_array< amp::ampf< Precision > > e, int n, bool isupper, bool isfractionalaccuracyrequired, ap::template_2d_array< amp::ampf< Precision > > &u, int ustart, int nru, ap::template_2d_array< amp::ampf< Precision > > &c, int cstart, int ncc, ap::template_2d_array< amp::ampf< Precision > > &vt, int vstart, int ncvt)

template<unsigned int Precision>
amp::ampf< Precision > | extsignbdsqr (amp::ampf< Precision > a, amp::ampf< Precision > b)
template<unsigned int Precision>
void | svd2x2 (amp::ampf< Precision > f, amp::ampf< Precision > g, amp::ampf< Precision > h, amp::ampf< Precision > &ssmin, amp::ampf< Precision > &ssmax)
template<unsigned int Precision>
void | svdv2x2 (amp::ampf< Precision > f, amp::ampf< Precision > g, amp::ampf< Precision > h, amp::ampf< Precision > &ssmin, amp::ampf< Precision > &ssmax, amp::ampf< Precision > &snr, amp::ampf< Precision > &csr, amp::ampf< Precision > &snl, amp::ampf< Precision > &csl)

## Function Documentation

## bidiagonalsvddecomposition()

template<unsigned int Precision>

bool bdsvd::bidiagonalsvddecomposition | ( | ap::template_1d_array< amp::ampf< Precision > > & | d,
|  | ap::template_1d_array< amp::ampf< Precision > > | e,
|  | int | n,
|  | bool | isupper,
|  | bool | isfractionalaccuracyrequired,
|  | ap::template_2d_array< amp::ampf< Precision > > & | u,
|  | int | nru,
|  | ap::template_2d_array< amp::ampf< Precision > > & | c,
|  | int | ncc,
|  | ap::template_2d_array< amp::ampf< Precision > > & | vt,
|  | int | ncvt )

Definition at line 229 of file bdsvd.h.

{
bool result;
result = bidiagonalsvddecompositioninternal<Precision>(d, e, n, isupper, isfractionalaccuracyrequired, u, 1, nru, c, 1, ncc, vt, 1, ncvt);
return result;
}
result
return result
Definition facAbsBiFact.cc:76
bdsvd::bidiagonalsvddecompositioninternal
bool bidiagonalsvddecompositioninternal(ap::template_1d_array< amp::ampf< Precision > > &d, ap::template_1d_array< amp::ampf< Precision > > e, int n, bool isupper, bool isfractionalaccuracyrequired, ap::template_2d_array< amp::ampf< Precision > > &u, int ustart, int nru, ap::template_2d_array< amp::ampf< Precision > > &c, int cstart, int ncc, ap::template_2d_array< amp::ampf< Precision > > &vt, int vstart, int ncvt)
Definition bdsvd.h:253

## bidiagonalsvddecompositioninternal()

template<unsigned int Precision>

bool bdsvd::bidiagonalsvddecompositioninternal | ( | ap::template_1d_array< amp::ampf< Precision > > & | d,
|  | ap::template_1d_array< amp::ampf< Precision > > | e,
|  | int | n,
|  | bool | isupper,
|  | bool | isfractionalaccuracyrequired,
|  | ap::template_2d_array< amp::ampf< Precision > > & | u,
|  | int | ustart,
|  | int | nru,
|  | ap::template_2d_array< amp::ampf< Precision > > & | c,
|  | int | cstart,
|  | int | ncc,
|  | ap::template_2d_array< amp::ampf< Precision > > & | vt,
|  | int | vstart,
|  | int | ncvt )

Definition at line 253 of file bdsvd.h.

{
bool result;
int i;
int idir;
int isub;
int iter;
int j;
int ll;
int lll;
int m;
int maxit;
int oldll;
int oldm;
amp::ampf<Precision> abse;
amp::ampf<Precision> abss;
amp::ampf<Precision> cosl;
amp::ampf<Precision> cosr;
amp::ampf<Precision> cs;
amp::ampf<Precision> eps;
amp::ampf<Precision> f;
amp::ampf<Precision> g;
amp::ampf<Precision> h;
amp::ampf<Precision> mu;
amp::ampf<Precision> oldcs;
amp::ampf<Precision> oldsn;
amp::ampf<Precision> r;
amp::ampf<Precision> shift;
amp::ampf<Precision> sigmn;
amp::ampf<Precision> sigmx;
amp::ampf<Precision> sinl;
amp::ampf<Precision> sinr;
amp::ampf<Precision> sll;
amp::ampf<Precision> smax;
amp::ampf<Precision> smin;
amp::ampf<Precision> sminl;
amp::ampf<Precision> sminlo;
amp::ampf<Precision> sminoa;
amp::ampf<Precision> sn;
amp::ampf<Precision> thresh;
amp::ampf<Precision> tol;
amp::ampf<Precision> tolmul;
amp::ampf<Precision> unfl;
ap::template_1d_array< amp::ampf<Precision> > work0;
ap::template_1d_array< amp::ampf<Precision> > work1;
ap::template_1d_array< amp::ampf<Precision> > work2;
ap::template_1d_array< amp::ampf<Precision> > work3;
int maxitr;
bool matrixsplitflag;
bool iterflag;
ap::template_1d_array< amp::ampf<Precision> > utemp;
ap::template_1d_array< amp::ampf<Precision> > vttemp;
ap::template_1d_array< amp::ampf<Precision> > ctemp;
ap::template_1d_array< amp::ampf<Precision> > etemp;
bool rightside;
bool fwddir;
amp::ampf<Precision> tmp;
int mm1;
int mm0;
bool bchangedir;
int uend;
int cend;
int vend;
result = true;
if( n==0 )
{
return result;
}
if( n==1 )
{
if( d(1)<0 )
{
d(1) = -d(1);
if( ncvt>0 )
{
ap::vmul(vt.getrow(vstart, vstart, vstart+ncvt-1), -1);
}
}
return result;
}
//
// init
//
work0.setbounds(1, n-1);
work1.setbounds(1, n-1);
work2.setbounds(1, n-1);
work3.setbounds(1, n-1);
uend = ustart+ap::maxint(nru-1, 0);
vend = vstart+ap::maxint(ncvt-1, 0);
cend = cstart+ap::maxint(ncc-1, 0);
utemp.setbounds(ustart, uend);
vttemp.setbounds(vstart, vend);
ctemp.setbounds(cstart, cend);
maxitr = 12;
rightside = true;
fwddir = true;
//
// resize E from N-1 to N
//
etemp.setbounds(1, n);
for(i=1; i<=n-1; i++)
{
etemp(i) = e(i);
}
e.setbounds(1, n);
for(i=1; i<=n-1; i++)
{
e(i) = etemp(i);
}
e(n) = 0;
idir = 0;
//
// Get machine constants
//
eps = amp::ampf<Precision>::getAlgoPascalEpsilon();
unfl = amp::ampf<Precision>::getAlgoPascalMinNumber();
//
// If matrix lower bidiagonal, rotate to be upper bidiagonal
// by applying Givens rotations on the left
//
if( !isupper )
{
for(i=1; i<=n-1; i++)
{
rotations::generaterotation<Precision>(d(i), e(i), cs, sn, r);
d(i) = r;
e(i) = sn*d(i+1);
d(i+1) = cs*d(i+1);
work0(i) = cs;
work1(i) = sn;
}
//
// Update singular vectors if desired
//
if( nru>0 )
{
rotations::applyrotationsfromtheright<Precision>(fwddir, ustart, uend, 1+ustart-1, n+ustart-1, work0, work1, u, utemp);
}
if( ncc>0 )
{
rotations::applyrotationsfromtheleft<Precision>(fwddir, 1+cstart-1, n+cstart-1, cstart, cend, work0, work1, c, ctemp);
}
}
//
// Compute singular values to relative accuracy TOL
// (By setting TOL to be negative, algorithm will compute
// singular values to absolute accuracy ABS(TOL)*norm(input matrix))
//
tolmul = amp::maximum<Precision>(10, amp::minimum<Precision>(100, amp::pow<Precision>(eps, -amp::ampf<Precision>("0.125"))));
tol = tolmul*eps;
if( !isfractionalaccuracyrequired )
{
tol = -tol;
}
//
// Compute approximate maximum, minimum singular values
//
smax = 0;
for(i=1; i<=n; i++)
{
smax = amp::maximum<Precision>(smax, amp::abs<Precision>(d(i)));
}
for(i=1; i<=n-1; i++)
{
smax = amp::maximum<Precision>(smax, amp::abs<Precision>(e(i)));
}
sminl = 0;
if( tol>=0 )
{
//
// Relative accuracy desired
//
sminoa = amp::abs<Precision>(d(1));
if( sminoa!=0 )
{
mu = sminoa;
for(i=2; i<=n; i++)
{
mu = amp::abs<Precision>(d(i))*(mu/(mu+amp::abs<Precision>(e(i-1))));
sminoa = amp::minimum<Precision>(sminoa, mu);
if( sminoa==0 )
{
break;
}
}
}
sminoa = sminoa/amp::sqrt<Precision>(n);
thresh = amp::maximum<Precision>(tol*sminoa, maxitr*n*n*unfl);
}
else
{
//
// Absolute accuracy desired
//
thresh = amp::maximum<Precision>(amp::abs<Precision>(tol)*smax, maxitr*n*n*unfl);
}
//
// Prepare for main iteration loop for the singular values
// (MAXIT is the maximum number of passes through the inner
// loop permitted before nonconvergence signalled.)
//
maxit = maxitr*n*n;
iter = 0;
oldll = -1;
oldm = -1;
//
// M points to last element of unconverged part of matrix
//
m = n;
//
// Begin main iteration loop
//
while( true )
{
//
// Check for convergence or exceeding iteration count
//
if( m<=1 )
{
break;
}
if( iter>maxit )
{
result = false;
return result;
}
//
// Find diagonal block of matrix to work on
//
if( tol<0 && amp::abs<Precision>(d(m))<=thresh )
{
d(m) = 0;
}
smax = amp::abs<Precision>(d(m));
smin = smax;
matrixsplitflag = false;
for(lll=1; lll<=m-1; lll++)
{
ll = m-lll;
abss = amp::abs<Precision>(d(ll));
abse = amp::abs<Precision>(e(ll));
if( tol<0 && abss<=thresh )
{
d(ll) = 0;
}
if( abse<=thresh )
{
matrixsplitflag = true;
break;
}
smin = amp::minimum<Precision>(smin, abss);
smax = amp::maximum<Precision>(smax, amp::maximum<Precision>(abss, abse));
}
if( !matrixsplitflag )
{
ll = 0;
}
else
{
//
// Matrix splits since E(LL) = 0
//
e(ll) = 0;
if( ll==m-1 )
{
//
// Convergence of bottom singular value, return to top of loop
//
m = m-1;
continue;
}
}
ll = ll+1;
//
// E(LL) through E(M-1) are nonzero, E(LL-1) is zero
//
if( ll==m-1 )
{
//
// 2 by 2 block, handle separately
//
svdv2x2<Precision>(d(m-1), e(m-1), d(m), sigmn, sigmx, sinr, cosr, sinl, cosl);
d(m-1) = sigmx;
e(m-1) = 0;
d(m) = sigmn;
//
// Compute singular vectors, if desired
//
if( ncvt>0 )
{
mm0 = m+(vstart-1);
mm1 = m-1+(vstart-1);
ap::vmove(vttemp.getvector(vstart, vend), vt.getrow(mm1, vstart, vend), cosr);
ap::vadd(vttemp.getvector(vstart, vend), vt.getrow(mm0, vstart, vend), sinr);
ap::vmul(vt.getrow(mm0, vstart, vend), cosr);
ap::vsub(vt.getrow(mm0, vstart, vend), vt.getrow(mm1, vstart, vend), sinr);
ap::vmove(vt.getrow(mm1, vstart, vend), vttemp.getvector(vstart, vend));
}
if( nru>0 )
{
mm0 = m+ustart-1;
mm1 = m-1+ustart-1;
ap::vmove(utemp.getvector(ustart, uend), u.getcolumn(mm1, ustart, uend), cosl);
ap::vadd(utemp.getvector(ustart, uend), u.getcolumn(mm0, ustart, uend), sinl);
ap::vmul(u.getcolumn(mm0, ustart, uend), cosl);
ap::vsub(u.getcolumn(mm0, ustart, uend), u.getcolumn(mm1, ustart, uend), sinl);
ap::vmove(u.getcolumn(mm1, ustart, uend), utemp.getvector(ustart, uend));
}
if( ncc>0 )
{
mm0 = m+cstart-1;
mm1 = m-1+cstart-1;
ap::vmove(ctemp.getvector(cstart, cend), c.getrow(mm1, cstart, cend), cosl);
ap::vadd(ctemp.getvector(cstart, cend), c.getrow(mm0, cstart, cend), sinl);
ap::vmul(c.getrow(mm0, cstart, cend), cosl);
ap::vsub(c.getrow(mm0, cstart, cend), c.getrow(mm1, cstart, cend), sinl);
ap::vmove(c.getrow(mm1, cstart, cend), ctemp.getvector(cstart, cend));
}
m = m-2;
continue;
}
//
// If working on new submatrix, choose shift direction
// (from larger end diagonal element towards smaller)
//
// Previously was
// "if (LL>OLDM) or (M<OLDLL) then"
// fixed thanks to Michael Rolle < m@rolle.name >
// Very strange that LAPACK still contains it.
//
bchangedir = false;
if( idir==1 && amp::abs<Precision>(d(ll))<amp::ampf<Precision>("1.0E-3")*amp::abs<Precision>(d(m)) )
{
bchangedir = true;
}
if( idir==2 && amp::abs<Precision>(d(m))<amp::ampf<Precision>("1.0E-3")*amp::abs<Precision>(d(ll)) )
{
bchangedir = true;
}
if( ll!=oldll || m!=oldm || bchangedir )
{
if( amp::abs<Precision>(d(ll))>=amp::abs<Precision>(d(m)) )
{
//
// Chase bulge from top (big end) to bottom (small end)
//
idir = 1;
}
else
{
//
// Chase bulge from bottom (big end) to top (small end)
//
idir = 2;
}
}
//
// Apply convergence tests
//
if( idir==1 )
{
//
// Run convergence test in forward direction
// First apply standard test to bottom of matrix
//
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 )
{
e(m-1) = 0;
continue;
}
if( tol>=0 )
{
//
// If relative accuracy desired,
// apply convergence criterion forward
//
mu = amp::abs<Precision>(d(ll));
sminl = mu;
iterflag = false;
for(lll=ll; lll<=m-1; lll++)
{
if( amp::abs<Precision>(e(lll))<=tol*mu )
{
e(lll) = 0;
iterflag = true;
break;
}
sminlo = sminl;
mu = amp::abs<Precision>(d(lll+1))*(mu/(mu+amp::abs<Precision>(e(lll))));
sminl = amp::minimum<Precision>(sminl, mu);
}
if( iterflag )
{
continue;
}
}
}
else
{
//
// Run convergence test in backward direction
// First apply standard test to top of matrix
//
if( amp::abs<Precision>(e(ll))<=amp::abs<Precision>(tol)*amp::abs<Precision>(d(ll)) || tol<0 && amp::abs<Precision>(e(ll))<=thresh )
{
e(ll) = 0;
continue;
}
if( tol>=0 )
{
//
// If relative accuracy desired,
// apply convergence criterion backward
//
mu = amp::abs<Precision>(d(m));
sminl = mu;
iterflag = false;
for(lll=m-1; lll>=ll; lll--)
{
if( amp::abs<Precision>(e(lll))<=tol*mu )
{
e(lll) = 0;
iterflag = true;
break;
}
sminlo = sminl;
mu = amp::abs<Precision>(d(lll))*(mu/(mu+amp::abs<Precision>(e(lll))));
sminl = amp::minimum<Precision>(sminl, mu);
}
if( iterflag )
{
continue;
}
}
}
oldll = ll;
oldm = m;
//
// Compute shift. First, test if shifting would ruin relative
// accuracy, and if so set the shift to zero.
//
if( tol>=0 && n*tol*(sminl/smax)<=amp::maximum<Precision>(eps, amp::ampf<Precision>("0.01")*tol) )
{
//
// Use a zero shift to avoid loss of relative accuracy
//
shift = 0;
}
else
{
//
// Compute the shift from 2-by-2 block at end of matrix
//
if( idir==1 )
{
sll = amp::abs<Precision>(d(ll));
svd2x2<Precision>(d(m-1), e(m-1), d(m), shift, r);
}
else
{
sll = amp::abs<Precision>(d(m));
svd2x2<Precision>(d(ll), e(ll), d(ll+1), shift, r);
}
//
// Test if shift negligible, and if so set to zero
//
if( sll>0 )
{
if( amp::sqr<Precision>(shift/sll)<eps )
{
shift = 0;
}
}
}
//
// Increment iteration count
//
iter = iter+m-ll;
//
// If SHIFT = 0, do simplified QR iteration
//
if( shift==0 )
{
if( idir==1 )
{
//
// Chase bulge from top to bottom
// Save cosines and sines for later singular vector updates
//
cs = 1;
oldcs = 1;
for(i=ll; i<=m-1; i++)
{
rotations::generaterotation<Precision>(d(i)*cs, e(i), cs, sn, r);
if( i>ll )
{
e(i-1) = oldsn*r;
}
rotations::generaterotation<Precision>(oldcs*r, d(i+1)*sn, oldcs, oldsn, tmp);
d(i) = tmp;
work0(i-ll+1) = cs;
work1(i-ll+1) = sn;
work2(i-ll+1) = oldcs;
work3(i-ll+1) = oldsn;
}
h = d(m)*cs;
d(m) = h*oldcs;
e(m-1) = h*oldsn;
//
// Update singular vectors
//
if( ncvt>0 )
{
rotations::applyrotationsfromtheleft<Precision>(fwddir, ll+vstart-1, m+vstart-1, vstart, vend, work0, work1, vt, vttemp);
}
if( nru>0 )
{
rotations::applyrotationsfromtheright<Precision>(fwddir, ustart, uend, ll+ustart-1, m+ustart-1, work2, work3, u, utemp);
}
if( ncc>0 )
{
rotations::applyrotationsfromtheleft<Precision>(fwddir, ll+cstart-1, m+cstart-1, cstart, cend, work2, work3, c, ctemp);
}
//
// Test convergence
//
if( amp::abs<Precision>(e(m-1))<=thresh )
{
e(m-1) = 0;
}
}
else
{
//
// Chase bulge from bottom to top
// Save cosines and sines for later singular vector updates
//
cs = 1;
oldcs = 1;
for(i=m; i>=ll+1; i--)
{
rotations::generaterotation<Precision>(d(i)*cs, e(i-1), cs, sn, r);
if( i<m )
{
e(i) = oldsn*r;
}
rotations::generaterotation<Precision>(oldcs*r, d(i-1)*sn, oldcs, oldsn, tmp);
d(i) = tmp;
work0(i-ll) = cs;
work1(i-ll) = -sn;
work2(i-ll) = oldcs;
work3(i-ll) = -oldsn;
}
h = d(ll)*cs;
d(ll) = h*oldcs;
e(ll) = h*oldsn;
//
// Update singular vectors
//
if( ncvt>0 )
{
rotations::applyrotationsfromtheleft<Precision>(!fwddir, ll+vstart-1, m+vstart-1, vstart, vend, work2, work3, vt, vttemp);
}
if( nru>0 )
{
rotations::applyrotationsfromtheright<Precision>(!fwddir, ustart, uend, ll+ustart-1, m+ustart-1, work0, work1, u, utemp);
}
if( ncc>0 )
{
rotations::applyrotationsfromtheleft<Precision>(!fwddir, ll+cstart-1, m+cstart-1, cstart, cend, work0, work1, c, ctemp);
}
//
// Test convergence
//
if( amp::abs<Precision>(e(ll))<=thresh )
{
e(ll) = 0;
}
}
}
else
{
//
// Use nonzero shift
//
if( idir==1 )
{
//
// Chase bulge from top to bottom
// Save cosines and sines for later singular vector updates
//
f = (amp::abs<Precision>(d(ll))-shift)*(extsignbdsqr<Precision>(1, d(ll))+shift/d(ll));
g = e(ll);
for(i=ll; i<=m-1; i++)
{
rotations::generaterotation<Precision>(f, g, cosr, sinr, r);
if( i>ll )
{
e(i-1) = r;
}
f = cosr*d(i)+sinr*e(i);
e(i) = cosr*e(i)-sinr*d(i);
g = sinr*d(i+1);
d(i+1) = cosr*d(i+1);
rotations::generaterotation<Precision>(f, g, cosl, sinl, r);
d(i) = r;
f = cosl*e(i)+sinl*d(i+1);
d(i+1) = cosl*d(i+1)-sinl*e(i);
if( i<m-1 )
{
g = sinl*e(i+1);
e(i+1) = cosl*e(i+1);
}
work0(i-ll+1) = cosr;
work1(i-ll+1) = sinr;
work2(i-ll+1) = cosl;
work3(i-ll+1) = sinl;
}
e(m-1) = f;
//
// Update singular vectors
//
if( ncvt>0 )
{
rotations::applyrotationsfromtheleft<Precision>(fwddir, ll+vstart-1, m+vstart-1, vstart, vend, work0, work1, vt, vttemp);
}
if( nru>0 )
{
rotations::applyrotationsfromtheright<Precision>(fwddir, ustart, uend, ll+ustart-1, m+ustart-1, work2, work3, u, utemp);
}
if( ncc>0 )
{
rotations::applyrotationsfromtheleft<Precision>(fwddir, ll+cstart-1, m+cstart-1, cstart, cend, work2, work3, c, ctemp);
}
//
// Test convergence
//
if( amp::abs<Precision>(e(m-1))<=thresh )
{
e(m-1) = 0;
}
}
else
{
//
// Chase bulge from bottom to top
// Save cosines and sines for later singular vector updates
//
f = (amp::abs<Precision>(d(m))-shift)*(extsignbdsqr<Precision>(1, d(m))+shift/d(m));
g = e(m-1);
for(i=m; i>=ll+1; i--)
{
rotations::generaterotation<Precision>(f, g, cosr, sinr, r);
if( i<m )
{
e(i) = r;
}
f = cosr*d(i)+sinr*e(i-1);
e(i-1) = cosr*e(i-1)-sinr*d(i);
g = sinr*d(i-1);
d(i-1) = cosr*d(i-1);
rotations::generaterotation<Precision>(f, g, cosl, sinl, r);
d(i) = r;
f = cosl*e(i-1)+sinl*d(i-1);
d(i-1) = cosl*d(i-1)-sinl*e(i-1);
if( i>ll+1 )
{
g = sinl*e(i-2);
e(i-2) = cosl*e(i-2);
}
work0(i-ll) = cosr;
work1(i-ll) = -sinr;
work2(i-ll) = cosl;
work3(i-ll) = -sinl;
}
e(ll) = f;
//
// Test convergence
//
if( amp::abs<Precision>(e(ll))<=thresh )
{
e(ll) = 0;
}
//
// Update singular vectors if desired
//
if( ncvt>0 )
{
rotations::applyrotationsfromtheleft<Precision>(!fwddir, ll+vstart-1, m+vstart-1, vstart, vend, work2, work3, vt, vttemp);
}
if( nru>0 )
{
rotations::applyrotationsfromtheright<Precision>(!fwddir, ustart, uend, ll+ustart-1, m+ustart-1, work0, work1, u, utemp);
}
if( ncc>0 )
{
rotations::applyrotationsfromtheleft<Precision>(!fwddir, ll+cstart-1, m+cstart-1, cstart, cend, work0, work1, c, ctemp);
}
}
}
//
// QR iteration finished, go back and check convergence
//
continue;
}
//
// All singular values converged, so make them positive
//
for(i=1; i<=n; i++)
{
if( d(i)<0 )
{
d(i) = -d(i);
//
// Change sign of singular vectors, if desired
//
if( ncvt>0 )
{
ap::vmul(vt.getrow(i+vstart-1, vstart, vend), -1);
}
}
}
//
// Sort the singular values into decreasing order (insertion sort on
// singular values, but only one transposition per singular vector)
//
for(i=1; i<=n-1; i++)
{
//
// Scan for smallest D(I)
//
isub = 1;
smin = d(1);
for(j=2; j<=n+1-i; j++)
{
if( d(j)<=smin )
{
isub = j;
smin = d(j);
}
}
if( isub!=n+1-i )
{
//
// Swap singular values and vectors
//
d(isub) = d(n+1-i);
d(n+1-i) = smin;
if( ncvt>0 )
{
j = n+1-i;
ap::vmove(vttemp.getvector(vstart, vend), vt.getrow(isub+vstart-1, vstart, vend));
ap::vmove(vt.getrow(isub+vstart-1, vstart, vend), vt.getrow(j+vstart-1, vstart, vend));
ap::vmove(vt.getrow(j+vstart-1, vstart, vend), vttemp.getvector(vstart, vend));
}
if( nru>0 )
{
j = n+1-i;
ap::vmove(utemp.getvector(ustart, uend), u.getcolumn(isub+ustart-1, ustart, uend));
ap::vmove(u.getcolumn(isub+ustart-1, ustart, uend), u.getcolumn(j+ustart-1, ustart, uend));
ap::vmove(u.getcolumn(j+ustart-1, ustart, uend), utemp.getvector(ustart, uend));
}
if( ncc>0 )
{
j = n+1-i;
ap::vmove(ctemp.getvector(cstart, cend), c.getrow(isub+cstart-1, cstart, cend));
ap::vmove(c.getrow(isub+cstart-1, cstart, cend), c.getrow(j+cstart-1, cstart, cend));
ap::vmove(c.getrow(j+cstart-1, cstart, cend), ctemp.getvector(cstart, cend));
}
}
}
return result;
}
m
int m
Definition cfEzgcd.cc:128
i
int i
Definition cfEzgcd.cc:132
g
g
Definition cfModGcd.cc:4098
f
FILE * f
Definition checklibs.c:9
amp::ampf
Definition amp.h:82
amp::ampf::getAlgoPascalEpsilon
static const ampf getAlgoPascalEpsilon()
Definition amp.h:565
amp::ampf::getAlgoPascalMinNumber
static const ampf getAlgoPascalMinNumber()
Definition amp.h:582
ap::template_1d_array
Definition ap.h:657
ap::template_1d_array::getvector
raw_vector< T > getvector(int iStart, int iEnd)
Definition ap.h:776
ap::template_1d_array::setbounds
void setbounds(int iLow, int iHigh)
Definition ap.h:735
ap::template_2d_array::getrow
raw_vector< T > getrow(int iRow, int iColumnStart, int iColumnEnd)
Definition ap.h:939
ap::template_2d_array::getcolumn
raw_vector< T > getcolumn(int iColumn, int iRowStart, int iRowEnd)
Definition ap.h:931
iter
CFFListIterator iter
Definition facAbsBiFact.cc:54
j
int j
Definition facHensel.cc:110
h
STATIC_VAR Poly * h
Definition janet.cc:971
mu
static matrix mu(matrix A, const ring R)
Definition matpol.cc:2028
amp::abs
const ampf< Precision > abs(const ampf< Precision > &x)
Definition amp.h:713
amp::maximum
const ampf< Precision > maximum(const ampf< Precision > &x, const ampf< Precision > &y)
Definition amp.h:722
amp::sqrt
const ampf< Precision > sqrt(const ampf< Precision > &x)
Definition amp.h:740
amp::pow
const ampf< Precision > pow(const ampf< Precision > &x, const ampf< Precision > &y)
Definition amp.h:1062
amp::sqr
const ampf< Precision > sqr(const ampf< Precision > &x)
Definition amp.h:693
amp::minimum
const ampf< Precision > minimum(const ampf< Precision > &x, const ampf< Precision > &y)
Definition amp.h:731
ap::maxint
int maxint(int m1, int m2)
Definition ap.cpp:162
ap::vmove
void vmove(raw_vector< T > vdst, const_raw_vector< T > vsrc)
Definition ap.h:237
ap::vadd
void vadd(raw_vector< T > vdst, const_raw_vector< T > vsrc)
Definition ap.h:413
ap::vmul
void vmul(raw_vector< T > vdst, T2 alpha)
Definition ap.h:603
ap::vsub
void vsub(raw_vector< T > vdst, const_raw_vector< T > vsrc)
Definition ap.h:533
bdsvd::extsignbdsqr
amp::ampf< Precision > extsignbdsqr(amp::ampf< Precision > a, amp::ampf< Precision > b)
Definition bdsvd.h:1096
bdsvd::svdv2x2
void svdv2x2(amp::ampf< Precision > f, amp::ampf< Precision > g, amp::ampf< Precision > h, amp::ampf< Precision > &ssmin, amp::ampf< Precision > &ssmax, amp::ampf< Precision > &snr, amp::ampf< Precision > &csr, amp::ampf< Precision > &snl, amp::ampf< Precision > &csl)
Definition bdsvd.h:1189
bdsvd::svd2x2
void svd2x2(amp::ampf< Precision > f, amp::ampf< Precision > g, amp::ampf< Precision > h, amp::ampf< Precision > &ssmin, amp::ampf< Precision > &ssmax)
Definition bdsvd.h:1115
rotations::applyrotationsfromtheright
void applyrotationsfromtheright(bool isforward, int m1, int m2, int n1, int n2, const ap::template_1d_array< amp::ampf< Precision > > &c, const ap::template_1d_array< amp::ampf< Precision > > &s, ap::template_2d_array< amp::ampf< Precision > > &a, ap::template_1d_array< amp::ampf< Precision > > &work)
Definition rotations.h:238
rotations::generaterotation
void generaterotation(amp::ampf< Precision > f, amp::ampf< Precision > g, amp::ampf< Precision > &cs, amp::ampf< Precision > &sn, amp::ampf< Precision > &r)
Definition rotations.h:355
rotations::applyrotationsfromtheleft
void applyrotationsfromtheleft(bool isforward, int m1, int m2, int n1, int n2, const ap::template_1d_array< amp::ampf< Precision > > &c, const ap::template_1d_array< amp::ampf< Precision > > &s, ap::template_2d_array< amp::ampf< Precision > > &a, ap::template_1d_array< amp::ampf< Precision > > &work)
Definition rotations.h:100

## extsignbdsqr()

template<unsigned int Precision>

amp::ampf< Precision > bdsvd::extsignbdsqr | ( | amp::ampf< Precision > | a,
|  | amp::ampf< Precision > | b )

Definition at line 1096 of file bdsvd.h.

{
amp::ampf<Precision> result;
if( b>=0 )
{
result = amp::abs<Precision>(a);
}
else
{
result = -amp::abs<Precision>(a);
}
return result;
}
b
CanonicalForm b
Definition cfModGcd.cc:4111

## rmatrixbdsvd()

template<unsigned int Precision>

bool bdsvd::rmatrixbdsvd | ( | ap::template_1d_array< amp::ampf< Precision > > & | d,
|  | ap::template_1d_array< amp::ampf< Precision > > | e,
|  | int | n,
|  | bool | isupper,
|  | bool | isfractionalaccuracyrequired,
|  | ap::template_2d_array< amp::ampf< Precision > > & | u,
|  | int | nru,
|  | ap::template_2d_array< amp::ampf< Precision > > & | c,
|  | int | ncc,
|  | ap::template_2d_array< amp::ampf< Precision > > & | vt,
|  | int | ncvt )

Definition at line 186 of file bdsvd.h.

{
bool result;
ap::template_1d_array< amp::ampf<Precision> > d1;
ap::template_1d_array< amp::ampf<Precision> > e1;
d1.setbounds(1, n);
ap::vmove(d1.getvector(1, n), d.getvector(0, n-1));
if( n>1 )
{
e1.setbounds(1, n-1);
ap::vmove(e1.getvector(1, n-1), e.getvector(0, n-2));
}
result = bidiagonalsvddecompositioninternal<Precision>(d1, e1, n, isupper, isfractionalaccuracyrequired, u, 0, nru, c, 0, ncc, vt, 0, ncvt);
ap::vmove(d.getvector(0, n-1), d1.getvector(1, n));
return result;
}

## svd2x2()

template<unsigned int Precision>

void bdsvd::svd2x2 | ( | amp::ampf< Precision > | f,
|  | amp::ampf< Precision > | g,
|  | amp::ampf< Precision > | h,
|  | amp::ampf< Precision > & | ssmin,
|  | amp::ampf< Precision > & | ssmax )

Definition at line 1115 of file bdsvd.h.

{
amp::ampf<Precision> aas;
amp::ampf<Precision> at;
amp::ampf<Precision> au;
amp::ampf<Precision> c;
amp::ampf<Precision> fa;
amp::ampf<Precision> fhmn;
amp::ampf<Precision> fhmx;
amp::ampf<Precision> ga;
amp::ampf<Precision> ha;
fa = amp::abs<Precision>(f);
ga = amp::abs<Precision>(g);
ha = amp::abs<Precision>(h);
fhmn = amp::minimum<Precision>(fa, ha);
fhmx = amp::maximum<Precision>(fa, ha);
if( fhmn==0 )
{
ssmin = 0;
if( fhmx==0 )
{
ssmax = ga;
}
else
{
ssmax = amp::maximum<Precision>(fhmx, ga)*amp::sqrt<Precision>(1+amp::sqr<Precision>(amp::minimum<Precision>(fhmx, ga)/amp::maximum<Precision>(fhmx, ga)));
}
}
else
{
if( ga<fhmx )
{
aas = 1+fhmn/fhmx;
at = (fhmx-fhmn)/fhmx;
au = amp::sqr<Precision>(ga/fhmx);
c = 2/(amp::sqrt<Precision>(aas*aas+au)+amp::sqrt<Precision>(at*at+au));
ssmin = fhmn*c;
ssmax = fhmx/c;
}
else
{
au = fhmx/ga;
if( au==0 )
{
//
// Avoid possible harmful underflow if exponent range
// asymmetric (true SSMIN may not underflow even if
// AU underflows)
//
ssmin = fhmn*fhmx/ga;
ssmax = ga;
}
else
{
aas = 1+fhmn/fhmx;
at = (fhmx-fhmn)/fhmx;
c = 1/(amp::sqrt<Precision>(1+amp::sqr<Precision>(aas*au))+amp::sqrt<Precision>(1+amp::sqr<Precision>(at*au)));
ssmin = fhmn*c*au;
ssmin = ssmin+ssmin;
ssmax = ga/(c+c);
}
}
}
}
fa
static BOOLEAN fa(leftv res, leftv args)
Definition cohomo.cc:3814

## svdv2x2()

template<unsigned int Precision>

void bdsvd::svdv2x2 | ( | amp::ampf< Precision > | f,
|  | amp::ampf< Precision > | g,
|  | amp::ampf< Precision > | h,
|  | amp::ampf< Precision > & | ssmin,
|  | amp::ampf< Precision > & | ssmax,
|  | amp::ampf< Precision > & | snr,
|  | amp::ampf< Precision > & | csr,
|  | amp::ampf< Precision > & | snl,
|  | amp::ampf< Precision > & | csl )

Definition at line 1189 of file bdsvd.h.

{
bool gasmal;
bool swp;
int pmax;
amp::ampf<Precision> a;
amp::ampf<Precision> clt;
amp::ampf<Precision> crt;
amp::ampf<Precision> d;
amp::ampf<Precision> fa;
amp::ampf<Precision> ft;
amp::ampf<Precision> ga;
amp::ampf<Precision> gt;
amp::ampf<Precision> ha;
amp::ampf<Precision> ht;
amp::ampf<Precision> l;
amp::ampf<Precision> m;
amp::ampf<Precision> mm;
amp::ampf<Precision> r;
amp::ampf<Precision> s;
amp::ampf<Precision> slt;
amp::ampf<Precision> srt;
amp::ampf<Precision> t;
amp::ampf<Precision> temp;
amp::ampf<Precision> tsign;
amp::ampf<Precision> tt;
amp::ampf<Precision> v;
ft = f;
fa = amp::abs<Precision>(ft);
ht = h;
ha = amp::abs<Precision>(h);
//
// PMAX points to the maximum absolute element of matrix
// PMAX = 1 if F largest in absolute values
// PMAX = 2 if G largest in absolute values
// PMAX = 3 if H largest in absolute values
//
pmax = 1;
swp = ha>fa;
if( swp )
{
//
// Now FA .ge. HA
//
pmax = 3;
temp = ft;
ft = ht;
ht = temp;
temp = fa;
fa = ha;
ha = temp;
}
gt = g;
ga = amp::abs<Precision>(gt);
if( ga==0 )
{
//
// Diagonal matrix
//
ssmin = ha;
ssmax = fa;
clt = 1;
crt = 1;
slt = 0;
srt = 0;
}
else
{
gasmal = true;
if( ga>fa )
{
pmax = 2;
if( fa/ga<amp::ampf<Precision>::getAlgoPascalEpsilon() )
{
//
// Case of very large GA
//
gasmal = false;
ssmax = ga;
if( ha>1 )
{
v = ga/ha;
ssmin = fa/v;
}
else
{
v = fa/ga;
ssmin = v*ha;
}
clt = 1;
slt = ht/gt;
srt = 1;
crt = ft/gt;
}
}
if( gasmal )
{
//
// Normal case
//
d = fa-ha;
if( d==fa )
{
l = 1;
}
else
{
l = d/fa;
}
m = gt/ft;
t = 2-l;
mm = m*m;
tt = t*t;
s = amp::sqrt<Precision>(tt+mm);
if( l==0 )
{
r = amp::abs<Precision>(m);
}
else
{
r = amp::sqrt<Precision>(l*l+mm);
}
a = amp::ampf<Precision>("0.5")*(s+r);
ssmin = ha/a;
ssmax = fa*a;
if( mm==0 )
{
//
// Note that M is very tiny
//
if( l==0 )
{
t = extsignbdsqr<Precision>(2, ft)*extsignbdsqr<Precision>(1, gt);
}
else
{
t = gt/extsignbdsqr<Precision>(d, ft)+m/t;
}
}
else
{
t = (m/(s+t)+m/(r+l))*(1+a);
}
l = amp::sqrt<Precision>(t*t+4);
crt = 2/l;
srt = t/l;
clt = (crt+srt*m)/a;
v = ht/ft;
slt = v*srt/a;
}
}
if( swp )
{
csl = srt;
snl = crt;
csr = slt;
snr = clt;
}
else
{
csl = clt;
snl = slt;
csr = crt;
snr = srt;
}
//
// Correct signs of SSMAX and SSMIN
//
if( pmax==1 )
{
tsign = extsignbdsqr<Precision>(1, csr)*extsignbdsqr<Precision>(1, csl)*extsignbdsqr<Precision>(1, f);
}
if( pmax==2 )
{
tsign = extsignbdsqr<Precision>(1, snr)*extsignbdsqr<Precision>(1, csl)*extsignbdsqr<Precision>(1, g);
}
if( pmax==3 )
{
tsign = extsignbdsqr<Precision>(1, snr)*extsignbdsqr<Precision>(1, snl)*extsignbdsqr<Precision>(1, h);
}
ssmax = extsignbdsqr<Precision>(ssmax, tsign);
ssmin = extsignbdsqr<Precision>(ssmin, tsign*extsignbdsqr<Precision>(1, f)*extsignbdsqr<Precision>(1, h));
}
l
int l
Definition cfEzgcd.cc:100
s
const CanonicalForm int s
Definition facAbsFact.cc:51
v
const Variable & v
< [in] a sqrfree bivariate poly
Definition facBivar.h:39