Line data    Source code

       1             : /* 
       2             :  * File:   CGSolver.cpp
       3             :  * Author: Pantelis Sopasakis
       4             :  * 
       5             :  * Created on November 11, 2015, 12:41 AM
       6             :  * 
       7             :  * ForBES is free software: you can redistribute it and/or modify
       8             :  * it under the terms of the GNU Lesser General Public License as published by
       9             :  * the Free Software Foundation, either version 3 of the License, or
      10             :  * (at your option) any later version.
      11             :  *  
      12             :  * ForBES is distributed in the hope that it will be useful,
      13             :  * but WITHOUT ANY WARRANTY; without even the implied warranty of
      14             :  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
      15             :  * GNU Lesser General Public License for more details.
      16             :  * 
      17             :  * You should have received a copy of the GNU Lesser General Public License
      18             :  * along with ForBES. If not, see <http://www.gnu.org/licenses/>.
      19             :  */
      20             : 
      21             : #include "CGSolver.h"
      22             : #include <lapacke.h>
      23             : #include <math.h>
      24             : #include <iostream>
      25             : 
      26           1 : CGSolver::CGSolver(LinearOperator& linop) : LinOpSolver(linop) {
      27           1 :     init();
      28           1 :     m_precond = NULL;
      29           1 : }
      30             : 
      31           0 : CGSolver::CGSolver(LinearOperator& linop, LinearOperator& preconditioner) : LinOpSolver(linop), m_precond(&preconditioner) {
      32           0 :     init();
      33           0 : }
      34             : 
      35           3 : CGSolver::CGSolver(LinearOperator& linop, LinearOperator& preconditioner, double tolerance, size_t max_iterations)
      36           3 : : LinOpSolver(linop), m_precond(&preconditioner), m_tolerance(tolerance), m_max_iterations(max_iterations) {
      37           3 :     m_err = NAN;
      38           3 :     m_iterations_count = 0;
      39           3 : }
      40             : 
      41           4 : CGSolver::~CGSolver() {
      42           4 : }
      43             : 
      44           1 : void CGSolver::init() {
      45           1 :     m_tolerance = 1e-4;
      46           1 :     m_max_iterations = 500;
      47           1 :     m_err = NAN;
      48           1 :     m_iterations_count = 0;
      49           1 : }
      50             : 
      51           4 : int CGSolver::solve(Matrix& b, Matrix& solution) {
      52           4 :     m_iterations_count = 0;
      53           4 :     m_err = NAN;
      54           4 :     Matrix r(b);
      55           8 :     Matrix z = m_precond != NULL ? m_precond->call(r) : MatrixFactory::ShallowMatrix(r);
      56           8 :     Matrix p = MatrixFactory::ShallowVector(z, 0);
      57           4 :     bool keepgoing = true;
      58           8 :     Matrix Ap;
      59         122 :     while (keepgoing) {
      60         117 :         Ap = m_linop->call(p);                      // Ap = A * p
      61         117 :         double a_denom = (p * Ap).get(0, 0);        //
      62         117 :         double a_numer = (r * z).get(0, 0);         //
      63         117 :         double alpha = a_numer / a_denom;           // alpha = (r,z)/(p, Ap);
      64         117 :         Matrix::add(solution, alpha, p, 1.0);       // x = x + alpha * p
      65         117 :         Matrix r_new = r;                           // r_new = r
      66             : 
      67             :         m_err = LAPACKE_dlange(
      68             :                 LAPACK_COL_MAJOR,
      69             :                 'I',
      70         117 :                 r_new.getNrows(),
      71             :                 1,
      72         117 :                 r_new.getData(),
      73         234 :                 r_new.getNrows());                  // compute ||r||_{inf}
      74             : 
      75         117 :         if (m_err < m_tolerance) {
      76           3 :             break;                                  // stop if tolerance reached
      77             :         }
      78             : 
      79         114 :         Matrix::add(r_new, -alpha, Ap, 1.0);        // r_new = r - alpha A p;
      80         114 :         Matrix z_new = m_precond != NULL 
      81          48 :              ? m_precond->call(r_new)
      82         276 :              : MatrixFactory::ShallowMatrix(r_new); // z_new = P(r_new)
      83         114 :         double zr = (r_new * z_new).get(0, 0);
      84         114 :         double beta = zr / a_numer;                 // beta = (r_new, z_nwq)/(r, z)
      85         114 :         Matrix::add(p, 1.0, z_new, beta);           // p = beta p + z_new
      86         114 :         z = z_new;                                  // z = z_new
      87         114 :         r = r_new;                                  // r = r_new
      88         114 :         m_iterations_count++;                       // k = k + 1
      89         114 :         if (m_iterations_count == m_max_iterations) {
      90           1 :             keepgoing = false;
      91             :         }
      92         114 :     }
      93           4 :     if (m_iterations_count>=m_max_iterations){
      94           1 :         return ForBESUtils::STATUS_MAX_ITERATIONS_REACHED;
      95             :     }
      96           7 :     return ForBESUtils::STATUS_OK;
      97             : }
      98             : 
      99           2 : double CGSolver::last_error() const {
     100           2 :     return m_err;
     101             : }
     102             : 
     103           3 : size_t CGSolver::last_num_iter() const {
     104           3 :     return m_iterations_count;
     105          30 : }
     106             : 
     107             : 
     108             :