MatrixSolver.h

Go to the documentation of this file.

//# MatrixSolver.h: the base class for solvers of AX=B
//# Copyright (C) 1994,1995,1999
//# Associated Universities, Inc. Washington DC, USA.
//# 
//# This library is free software; you can redistribute it and/or modify it
//# under the terms of the GNU Library General Public License as published by
//# the Free Software Foundation; either version 2 of the License, or (at your
//# option) any later version.
//# 
//# This library 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 Library General Public
//# License for more details.
//# 
//# You should have received a copy of the GNU Library General Public License
//# along with this library; if not, write to the Free Software Foundation,
//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
//# 
//# Correspondence concerning AIPS++ should be addressed as follows:
//#        Internet email: casa-feedback@nrao.edu.
//#        Postal address: AIPS++ Project Office
//#                        National Radio Astronomy Observatory
//#                        520 Edgemont Road
//#                        Charlottesville, VA 22903-2475 USA
 
#ifndef SCIMATH_MATRIXSOLVER_H
#define SCIMATH_MATRIXSOLVER_H
 
 
#include <casacore/casa/aips.h>
#include <casacore/casa/Arrays/Array.h>
#include <casacore/casa/Arrays/ArrayMath.h>
#include <casacore/casa/Arrays/Matrix.h>
#include <casacore/casa/Arrays/Vector.h>
 
#include <casacore/casa/Logging/LogSink.h>
#include <casacore/casa/Logging/LogMessage.h>
 
namespace casacore { //# NAMESPACE CASACORE - BEGIN
 
typedef Float FType;  // floating type (Float, Double)
 
//<summary>
// MatrixSolver.h: the base class for solvers of linear equations AX=B
//</summary>
 
// <use visibility=local>
 
// <reviewed reviewer="" date="",tests="" demos="">
// </reviewed>
 
// <prerequisite>
//   <li> Matrix, Vector
// </prerequisite>
//
// <etymology>
// The MatrixSolver class name reflects its use as the base class for solving
// Linear Equations of the form AX=B. This class is purely virtual
// and provides the essential implementation for derived solvers
// classes.  
// </etymology>
//
// <synopsis> 
// The MatrixSolver class is a purely virtual base class.  The programmer needs
// to define the following functions in a class derived from MatrixSolver:
// <ol>
//  <li> the derived destructor.
//  <li> <src>void setImageAndPsf(const Array<FType> & image, const 
//  Array<FType> & psf);</src> Set the image and the Point Spread Function 
//  (beam).  Setting this should reset the internal state, e.g. 
//  CurrentIter()==0.
//  <li> <src>Bool solve();</src>  Perform solution of AX=B.
//       Returns True if algorithm has converged or stop criterium reached.
// </ol>
// </synopsis> 
//
// <todo asof="">
// </todo>
 
class MatrixSolver {
public:
 
  // Default Constructor
  MatrixSolver();
  
  // Copy Constructor
  MatrixSolver(const MatrixSolver & other);
  
  // Create a MatrixSolver from a matrix A and a Vector B
  // <note role=warning> A and B are accessed by reference, so do not 
  // modify them during the lifetime of the MatrixSolver </note>
  MatrixSolver(const Matrix<FType> & A, const Vector<FType> & B);
  
  // Virtual destructor: calls all derived class destructors
  virtual ~MatrixSolver();
  
  // Assignment operator: uses reference semantics, i.e., it 
  // references the internal arrays of other
  MatrixSolver & operator=(const MatrixSolver & other);
 
  // Set A matrix and B vector
  void setAB(const Matrix<FType> & A, const Vector<FType> & B);
 
  // Set initial value of X
  void setX(const Vector<FType> & X);
  
  // Solve for the X vector.
  virtual Bool solve();
 
  // Is the current solution good enough?
  Bool accurateSolution();
 
  // Return residual vector B-AX
  const Vector<FType> & getResidual();
  
  // Return solution vector
  const Vector<FType> & getSolution();
 
  // Set the tolerance for solution
  void setTolerance(FType tol);
 
  // Return the tolerance for solution
  FType Tolerance();
  
  // Set the maximum number of iterations.
  void setMaxIters(uInt maxiters);
  
  // Return the maximum number of iterations.
  uInt MaxIters();
 
  // Set the gain for solution
  void setGain(FType g);
 
  // Return the gain for solution
  FType Gain();
 
  // Set status of solution
  void setSolved(Bool s);
 
  // Return status of solution
  Bool Solved();
 
  // Return norm of solution i.e. ||B-AX||
  FType getNorm();
  
protected:
 
  LogSink logSink_p;
  virtual LogSink& logSink() {return logSink_p;}
 
  // the A matrix data member
  Matrix<FType> AMatrix;
 
  // the constraint vector data member
  Vector<FType> BVector;
 
  // The residual vector data member
  Vector<FType> RVector;
 
  // The solution vector data member
  Vector<FType> XVector;
 
  // The solution norm i.e. ||B-AX||
  FType RNorm;
 
  // The data norm i.e. ||B||
  FType BNorm;
 
private:
 
  // Tolerance for solution i.e. ||B-AX||/||B|| must be less than this
  FType SolTolerance; 
 
  // Maximum number of iterations
  uInt MaxIterations;
 
  // Has a solution been found?
  Bool solved;
 
  // Gain
  FType gain;
 
};
 
inline void MatrixSolver::setTolerance(FType tol) 
{SolTolerance=tol;}
 
inline FType MatrixSolver::Tolerance() 
{return SolTolerance;}
 
inline void MatrixSolver::setMaxIters(uInt maxiters) 
{MaxIterations = maxiters;}
 
inline uInt MatrixSolver::MaxIters() 
{return MaxIterations;}
 
inline void MatrixSolver::setGain(FType g) 
{gain=g;}
 
inline FType MatrixSolver::Gain() 
{return gain;}
 
inline void MatrixSolver::setSolved(Bool s) 
{solved=s;}
 
inline Bool MatrixSolver::Solved() 
{return solved;}
 
inline FType MatrixSolver::getNorm()
{return RNorm;}
 
 
} //# NAMESPACE CASACORE - END
 
#endif