Commit e088e3d2 by Francois Gygi

cleanup use of Matrix::polar in ortho_align

git-svn-id: http://qboxcode.org/svn/qb/trunk@1186 cba15fb0-1239-40c8-b417-11db7ca47a34
parent 796e919d
......@@ -765,110 +765,9 @@ void SlaterDet::ortho_align(const SlaterDet& sd)
tmap["polar"].reset();
tmap["polar"].start();
#endif
#define MATRIX_POLAR 1
#if MATRIX_POLAR
x.polar();
#else
DoubleMatrix xp(ctxt_,c_.n(),c_.n(),c_.nb(),c_.nb());
DoubleMatrix t(ctxt_,c_.n(),c_.n(),c_.nb(),c_.nb());
DoubleMatrix q(ctxt_,c_.n(),c_.n(),c_.nb(),c_.nb());
DoubleMatrix qt(ctxt_,c_.n(),c_.n(),c_.nb(),c_.nb());
#if 1
// compute the polar decomposition of X = L^-1 B
// Bjork-Bowie version
//
const int niter = 3;
for ( int iter = 0; iter < niter; iter++ )
{
// q = I - x^T * x
q.identity();
q.syrk('l','t',-1.0,x,1.0);
q.symmetrize('l');
cout << " q.nrm2 = " << q.nrm2() << endl;
// compute xp = x * ( I + 0.5*q * ( I + 0.75 * q ) )
// t = ( I + 0.75 * q )
t.identity();
t.axpy(0.75,q);
// compute q*t
qt.gemm('n','n',1.0,q,t,0.0);
// pk = x * ( I + 0.5*q * ( I + 0.75 * q ) )
// = x * ( I + 0.5 * qt )
// Use t to store (I + 0.5 * qt)
t.identity();
t.axpy(0.5,qt);
// t now contains (I + 0.5 * qt)
// xp = x * t
xp.gemm('n','n',1.0,x,t,0.0);
// update x
x = xp;
}
// at convergence, x contains the polar factor of L^-1 B
#else
// Higham version
// compute the polar decomposition of L^-1 B
double diff = 1.0;
const double epsilon = 1.e-10;
const int maxiter = 3;
int iter = 0;
#if TIMING
tmap["transpose"].reset();
tmap["inverse"].reset();
#endif
while ( iter < maxiter && diff > epsilon )
{
// t = X^T
#if TIMING
tmap["transpose"].start();
#endif
t.transpose(1.0,x,0.0);
#if TIMING
tmap["transpose"].stop();
tmap["inverse"].start();
#endif
t.inverse();
#if TIMING
tmap["inverse"].stop();
#endif
// t now contains X^-T
// xp = 0.5 * ( x + x^-T );
for ( int i = 0; i < x.size(); i++ )
xp[i] = 0.5 * ( x[i] + t[i] );
// Next lines: use t as temporary to compute || x - xp ||_F
for ( int i = 0; i < t.size(); i++ )
t[i] = x[i] - xp[i];
#if TIMING
tmap["nrm2"].start();
#endif
diff = t.nrm2();
#if TIMING
tmap["nrm2"].stop();
#endif
//cout << " SlaterDet::ortho_align: diff=" << diff << endl;
x = xp;
//cout << "SlaterDet::ortho_align: X=\n" << x << endl;
iter++;
}
#endif // Bjork-Bowie vs Higham
#endif // MATRIX_POLAR
const double tol = 1.e-6;
const int maxiter = 2;
x.polar(tol,maxiter);
#if TIMING
tmap["polar"].stop();
#endif
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment