CUDA Lab Course Reference Manual 2011

/Users/sck/cuda/Praktikum_2011/include/lac/solver_cg_cuda.h

Go to the documentation of this file.

#ifndef SOLVER_CG_H
#define SOLVER_CG_H


#include <lac/solver_cuda.h>
#include <lac/tridiagonal_matrix.h>

template <class VECTOR = ::Vector<double> >
class SolverCG : public Solver<VECTOR>
{
  public:
    struct AdditionalData
    {
        bool log_coefficients;

        bool compute_condition_number;

        bool compute_all_condition_numbers;

        bool compute_eigenvalues;

        AdditionalData (const bool log_coefficients = false,
                        const bool compute_condition_number = false,
                        const bool compute_all_condition_numbers = false,
                        const bool compute_eigenvalues = false);
    };

    SolverCG (::SolverControl        &cn,
              ::VectorMemory<VECTOR> &mem,
              const AdditionalData &data = AdditionalData());

    SolverCG (::SolverControl        &cn,
              const AdditionalData &data=AdditionalData());

    virtual ~SolverCG ();

    template <class MATRIX, class PRECONDITIONER>
    void
    solve (const MATRIX         &A,
           VECTOR               &x,
           const VECTOR         &b,
           const PRECONDITIONER &precondition);

  protected:
    virtual double criterion();

    virtual void print_vectors(const unsigned int step,
                               const VECTOR& x,
                               const VECTOR& r,
                               const VECTOR& d) const;

    VECTOR *Vr;
    VECTOR *Vp;
    VECTOR *Vz;
    VECTOR *VAp;

    double res2;

    AdditionalData additional_data;

  private:
    void cleanup();
};

/*------------------------- Implementation ----------------------------*/


template <class VECTOR>
inline
SolverCG<VECTOR>::AdditionalData::
AdditionalData (const bool log_coefficients,
                const bool compute_condition_number,
                const bool compute_all_condition_numbers,
                const bool compute_eigenvalues)
                :
                log_coefficients (log_coefficients),
                compute_condition_number(compute_condition_number),
                compute_all_condition_numbers(compute_all_condition_numbers),
                compute_eigenvalues(compute_eigenvalues)
{}



template <class VECTOR>
SolverCG<VECTOR>::SolverCG (::SolverControl        &cn,
                            ::VectorMemory<VECTOR> &mem,
                            const AdditionalData &data)
                :
                Solver<VECTOR>(cn,mem),
                additional_data(data)
{}



template <class VECTOR>
SolverCG<VECTOR>::SolverCG (::SolverControl        &cn,
                            const AdditionalData &data)
                :
                Solver<VECTOR>(cn),
                additional_data(data)
{}



template <class VECTOR>
SolverCG<VECTOR>::~SolverCG ()
{}



template <class VECTOR>
double
SolverCG<VECTOR>::criterion()
{
  return std::sqrt(res2);
}



template <class VECTOR>
void
SolverCG<VECTOR>::cleanup()
{
  this->memory.free(Vr);
  this->memory.free(Vp);
  this->memory.free(Vz);
  this->memory.free(VAp);
  ::deallog.pop();
}



template <class VECTOR>
void
SolverCG<VECTOR>::print_vectors(const unsigned int,
                                const VECTOR&,
                                const VECTOR&,
                                const VECTOR&) const
{}



template <class VECTOR>
template <class MATRIX, class PRECONDITIONER>
void
SolverCG<VECTOR>::solve (const MATRIX         &A,
                         VECTOR               &x,
                         const VECTOR         &b,
                         const PRECONDITIONER &precondition)
{
  ::SolverControl::State conv=::SolverControl::iterate;

  ::deallog.push("cg");

                                   // Memory allocation
  Vr  = this->memory.alloc();
  Vp  = this->memory.alloc();
  Vz  = this->memory.alloc();
  VAp = this->memory.alloc();
                                   // Should we build the matrix for
                                   // eigenvalue computations?
  bool do_eigenvalues = additional_data.compute_condition_number
                        | additional_data.compute_all_condition_numbers
                        | additional_data.compute_eigenvalues;
  double eigen_beta_alpha = 0;

                                   // vectors used for eigenvalue
                                   // computations
  std::vector<double> diagonal;
  std::vector<double> offdiagonal;

  try {
                                     // define some aliases for simpler access
    VECTOR& g  = *Vr;
    VECTOR& h  = *Vp;
    VECTOR& d  = *Vz;
    VECTOR& Ad = *VAp;
                                     // resize the vectors, but do not set
                                     // the values since they'd be overwritten
                                     // soon anyway.
    g.reinit(x, true);
    h.reinit(x, true);
    d.reinit(x, true);
    Ad.reinit(x, true);
                                     // Implementation taken from the DEAL
                                     // library
    int  it=0;
    double res,gh,alpha,beta;

                                     // compute residual. if vector is
                                     // zero, then short-circuit the
                                     // full computation
    if (!x.all_zero())
      {
        A.vmult(g,x);
        g.sadd(-1.,1.,b);
      }
    else
      g = b;
    res = g.l2_norm();

    conv = this->control().check(0,res);
    if (conv)
      {
        cleanup();
        return;
      }

    g *= -1.;
    precondition.vmult(h,g);

    d.equ(-1.,h);

    gh = g*h;

    while (conv == ::SolverControl::iterate)
      {
        it++;
        A.vmult(Ad,d);

        alpha = d*Ad;
        alpha = gh/alpha;

        g.add(alpha,Ad);
        x.add(alpha,d );
        res = g.l2_norm();

        print_vectors(it, x, g, d);

        conv = this->control().check(it,res);
        if (conv)
          break;

        precondition.vmult(h,g);

        beta = gh;
        gh   = g*h;
        beta = gh/beta;

        if (additional_data.log_coefficients)
          ::deallog << "alpha-beta:" << alpha << '\t' << beta << std::endl;
                                         // set up the vectors
                                         // containing the diagonal
                                         // and the off diagonal of
                                         // the projected matrix.
        if (do_eigenvalues)
          {
            diagonal.push_back(1./alpha + eigen_beta_alpha);
            eigen_beta_alpha = beta/alpha;
            offdiagonal.push_back(std::sqrt(beta)/alpha);
          }

        if (additional_data.compute_all_condition_numbers && (diagonal.size()>1))
          {
            ::TridiagonalMatrix<double> T(diagonal.size(), true);
            for (unsigned int i=0;i<diagonal.size();++i)
              {
                T(i,i) = diagonal[i];
                if (i< diagonal.size()-1)
                  T(i,i+1) = offdiagonal[i];
              }
            T.compute_eigenvalues();
            ::deallog << "Condition number estimate: " <<
              T.eigenvalue(T.n()-1)/T.eigenvalue(0) << std::endl;
          }

        d.sadd(beta,-1.,h);
      }
  }
  catch (...)
    {
      cleanup();
      throw;
    }

                                   // Write eigenvalues or condition number
  if (do_eigenvalues)
    {
      ::TridiagonalMatrix<double> T(diagonal.size(), true);
      for (unsigned int i=0;i<diagonal.size();++i)
        {
          T(i,i) = diagonal[i];
          if (i< diagonal.size()-1)
            T(i,i+1) = offdiagonal[i];
        }
      T.compute_eigenvalues();
      if (additional_data.compute_condition_number
          && ! additional_data.compute_all_condition_numbers
          && (diagonal.size() > 1))
        ::deallog << "Condition number estimate: " <<
          T.eigenvalue(T.n()-1)/T.eigenvalue(0) << std::endl;
      if (additional_data.compute_eigenvalues)
        {
          for (unsigned int i=0;i<T.n();++i)
            ::deallog << ' ' << T.eigenvalue(i);
          ::deallog << std::endl;
        }
    }

                                   // Deallocate Memory
  cleanup();
                                   // in case of failure: throw
                                   // exception
  if (this->control().last_check() != ::SolverControl::success)
    throw ::SolverControl::NoConvergence (this->control().last_step(),
                                        this->control().last_value());
                                   // otherwise exit as normal
}


#endif // SOLVER_CG_H

---