opt_CG.h

Go to the documentation of this file.

#ifndef OPT_CG_H
#define OPT_CG_H
 
#include <math.h>
 
#include <iostream>
 
#include "global_function.h"
#include "source_base/parallel_reduce.h"
 
namespace ModuleBase
{
class Opt_CG
{
  public:
    Opt_CG();
    ~Opt_CG();
 
    void init_b(double* pinp_b // b in the linear equation Ax = b
    );
    void allocate(int nx // length of the solution array x
    );
    void set_para(double dV);
    void refresh(int nx_new = 0,          // length of new x, default 0 means the length doesn't change
                 double* pinp_b = nullptr // new b in Ax = b, default nullptr means we are dealing with general case
    );
 
    void next_direct(double* pgradient, // Ad for linear equaiont Ax=b, and gradient for general case
                     int label,         // 0 for solve Ax=b, 1 for PR form, 2 for HZ form
                     double* rdirect    // next direct
    );
    double step_length(double* pAd,     // Ad for Ax=b
                       double* pdirect, // direct
                       int& ifPD        // if postive definit
    );
 
    double get_residual()
    {
        return sqrt(this->gg_);
    };
    int get_iter()
    {
        return this->iter_;
    }
 
    // void ZEROS(double *x, int n)
    // {
    //     for (int i = 0; i < n; ++i) x[i] =0;
    // }
 
  private:
    double dV_ = 1.;
    int nx_ = 0;                      // length of the sulotion array x
    int iter_ = 0;                    // number of iteration
    double gg_ = 1000;                // gradient dot gradient
    double beta_ = 0.;                // d = -g + beta * d
    double eta_ = 0.01;               // a constand used in HZ form
    double* pdirect_old_ = nullptr;   // direction of last step
    double* pgradient_old_ = nullptr; // gradient, for meth=0, gradient is minus residual r.
 
    // only for standard CG
    double alpha_ = 0.;    // step length in standard CG
    double* pb_ = nullptr; // b in Ax=b, only for standard CG
 
    void stantard_CGdirect(double* pAd,    // Ad for Ax=b
                           double* rdirect // next direct
    );
    void PR_beta(double* pgradient // df(x)/dx
    );
    void HZ_beta(double* pgradient // df(x)/dx
    );
    double inner_product(double* pa, double* pb, int length)
    {
        double innerproduct = BlasConnector::dot(length, pa, 1, pb, 1);
        innerproduct *= this->dV_;
        return innerproduct;
    }
};
} // namespace ModuleBase
#endif