LSsystem.h

//===========================================================================
// SINTEF LSMG library - version 1.1
//
// Copyright (C) 2000-2005 SINTEF ICT, Applied Mathematics, Norway.
//
// This program is free software; you can redistribute it and/or          
// modify it under the terms of the GNU General Public License            
// as published by the Free Software Foundation version 2 of the License. 
//
// This program is distributed in the hope that it will be useful,        
// but WITHOUT ANY WARRANTY; without even the implied warranty of         
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          
// GNU General Public License for more details.                           
//
// You should have received a copy of the GNU General Public License      
// along with this program; if not, write to the Free Software            
// Foundation, Inc.,                                                      
// 59 Temple Place - Suite 330,                                           
// Boston, MA  02111-1307, USA.                                           
//
// Contact information: e-mail: tor.dokken@sintef.no                      
// SINTEF ICT, Department of Applied Mathematics,                         
// P.O. Box 124 Blindern,                                                 
// 0314 Oslo, Norway.                                                     
//
// Other licenses are also available for this software, notably licenses
// for:
// - Building commercial software.                                        
// - Building software whose source code you wish to keep private.        
//===========================================================================
#ifndef _LSSYSTEM_H_
#define _LSSYSTEM_H_

#include <MatSparse.h>
#include <GenMatrix.h>
#include <UCBtypedef.h>

#include <map>

#include <boost/shared_ptr.hpp>

class LSsystem {
  friend class MGsystem;
  

  // Equation system A_ x_ = b_
  MatSparse A_;                                  // The system matrix, with smoothing term included
  boost::shared_ptr<GenMatrix<UCBspl_real> > x_; // The unknown B-spline coefficient matrix (corresponds to PHI_ in MBA)
  GenMatrix<UCBspl_real> b_;                     // Right hand side (also as a matrix) built from z-values

  int n1_, n2_;      // size of spline coefficient matrix: x_.noX(), x_.noY() in A_ x_= b_
  int K_,   L_;      // order of spline surface
  
  double lambda_;    // calculated from the l_2 norms of the matrices: 
  double lambdaFac_; // a factor multiplied with lambda (default = 1)
  
  // The scattered data
  boost::shared_ptr<std::vector<double> > U_;
  boost::shared_ptr<std::vector<double> > V_;
  boost::shared_ptr<std::vector<double> > Zvals_;
  boost::shared_ptr<std::map<int, double, std::less<int> > > weights_;


  std::vector<double> domain_;   // defaults is umin, vmin, umax, vmax
    
  // Internal (private) functions
  void   calculateDomain();
  double umin() const {return domain_[0];}
  double vmin() const {return domain_[1];}
  double umax() const {return domain_[2];}
  double vmax() const {return domain_[3];}
  void   buildSparseStructure();
  void   addSmoothingTerm(double lambdaVal); // (= regularization term) (to A_)
  void   setLimits();
  void   calcLambda(); // for smoothing term

  void buildRHS(); // Build right hand side from data (this is also optionally done in buildEqSystem)

  void residual(GenMatrix<UCBspl_real>& r) const;

public:
  
  LSsystem(boost::shared_ptr<std::vector<double> > U,
           boost::shared_ptr<std::vector<double> > V,
           boost::shared_ptr<std::vector<double> > Zvals,
           int noX, int noY,       // grid size and number of unknowns is noX*noY
           boost::shared_ptr<GenMatrix<UCBspl_real> > PHI, // solution vector
           double smoothingFac = 1.0);

  ~LSsystem(){}
  
  void setWeights(boost::shared_ptr<std::map<int, double, std::less<int> > > weights) {weights_ = weights;}

  void buildEqSystem(bool buildRHS = true); // A_ and b_ (and smoothing term)
  
  void relaxGaussSeidel(int noIterations);
  
  void relaxCG(int noIterations);

  void relaxCGPrecond(int noIterations);

  void relaxJacobi(int noIterations, double omega=2./3.);

  void setDomain(double umin, double vmin, double umax, double vmax);

  void getDomain(double& umin, double& vmin, double& umax, double& vmax);

  int  numberOfUnknowns() const {return x_->noX() * x_->noY();} // Which is the same as n1_*n2_

  void setSmoothingFactor(double smoothingFac);

  bool addPoint(double u, double v, double z, double weight=1.0, bool addToRHS=true);
    
  double currentSolutionNorm() const;

  double residual_l2(bool scaled=false) const;
  double residual_linf() const;

  static double rowMatVecMult(int row_no, const MatSparse& A, const GenMatrix<UCBspl_real>& x);
};

#endif

---