12345qiupeng
/
gtsam
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Add the citation of the reference.

release/4.3a0
Sungtae An 2014-12-05 15:33:22 -05:00
parent 06b82ce3e3
commit 1c9a1f7cdb
1 changed files with 14 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

24
gtsam/linear/ConjugateGradientSolver.h
View File

@ -93,6 +93,10 @@ public:
* System class should support residual(v, g), multiply(v,Av), scal(alpha,v), dot(v,v), axpy(alpha,x,y)
* leftPrecondition(v, L^{-1}v, rightPrecondition(v, L^{-T}v) where preconditioner M = L*L^T
* Note that the residual is in the preconditioned domain. Refer to Section 9.2 of Saad's book.
*
** REFERENCES:
* [1] Y. Saad, "Preconditioned Iterations," in Iterative Methods for Sparse Linear Systems,
* 2nd ed. SIAM, 2003, ch. 9, sec. 2, pp.276-281.
*/
template <class S, class V>
V preconditionedConjugateGradient(const S &system, const V &initial, const ConjugateGradientParameters &parameters) {
@ -101,8 +105,8 @@ V preconditionedConjugateGradient(const S &system, const V &initial, const Conju
estimate = residual = direction = q1 = q2 = initial;
system.residual(estimate, q1); /* q1 = b-Ax */
system.leftPrecondition(q1, residual); /* r = S^{-T} (b-Ax) */
system.rightPrecondition(residual, direction);/* d = S^{-1} r */
system.leftPrecondition(q1, residual); /* r = L^{-1} (b-Ax) */
system.rightPrecondition(residual, direction);/* p = L^{-T} r */
double currentGamma = system.dot(residual, residual), prevGamma, alpha, beta;
@ -125,20 +129,20 @@ V preconditionedConjugateGradient(const S &system, const V &initial, const Conju
if ( k % iReset == 0 ) {
system.residual(estimate, q1); /* q1 = b-Ax */
system.leftPrecondition(q1, residual); /* r = L^{-1} (b-Ax) */
system.rightPrecondition(residual, direction); /* d = L^{-T} r */
system.rightPrecondition(residual, direction); /* p = L^{-T} r */
currentGamma = system.dot(residual, residual);
}
system.multiply(direction, q1); /* q1 = A d */
alpha = currentGamma / system.dot(direction, q1); /* alpha = gamma / (d' A d) */
system.axpy(alpha, direction, estimate); /* estimate += alpha * direction */
system.multiply(direction, q1); /* q1 = A p */
alpha = currentGamma / system.dot(direction, q1); /* alpha = gamma / (p' A p) */
system.axpy(alpha, direction, estimate); /* estimate += alpha * p */
system.leftPrecondition(q1, q2); /* q2 = L^{-1} * q1 */
system.axpy(-alpha, q2, residual); /* residual -= alpha * q2 */
system.axpy(-alpha, q2, residual); /* r -= alpha * q2 */
prevGamma = currentGamma;
currentGamma = system.dot(residual, residual); /* gamma = |residual|^2 */
currentGamma = system.dot(residual, residual); /* gamma = |r|^2 */
beta = currentGamma / prevGamma;
system.rightPrecondition(residual, q1); /* q1 = L^{-T} residual */
system.rightPrecondition(residual, q1); /* q1 = L^{-T} r */
system.scal(beta, direction);
system.axpy(1.0, q1, direction); /* direction = q1 + beta * direction */
system.axpy(1.0, q1, direction); /* p = q1 + beta * p */
if (parameters.verbosity() >= ConjugateGradientParameters::ERROR )
std::cout << "[PCG] k = " << k