stationary.h

Go to the documentation of this file.

// ========================================================================= //
// This file is part of MyraMath, copyright (c) 2014-2019 by Ryan A Chilton  //
// and distributed by MyraCore, LLC. See LICENSE.txt for license terms.      //
// ========================================================================= //

#ifndef MYRAMATH_ITERATIVE_STATIONARY_H
#define MYRAMATH_ITERATIVE_STATIONARY_H

#include <myramath/MYRAMATH_EXPORT.h>
#include <myramath/utility/Number.h>

#include <vector>

namespace myra {

// Forward declarations.
template<class Number> class CSparseMatrixRange;
template<class Number> class CVectorRange;
template<class Number> class VectorRange;

template<class Precision> class MYRAMATH_EXPORT stationary_output
  {
  public:

    // Did the relative residual converge to requested tolerance?
    bool converged;
  
    // History of relative residual, |A*x-b|/|b|
    std::vector<Precision> history;

    // Constructs from components.
    stationary_output(bool in_converged, const std::vector<Precision>& in_history)
      {
      converged = in_converged;
      history = in_history;
      }

  };

MYRAMATH_EXPORT stationary_output<NumberS> jacobi(const CSparseMatrixRange<NumberS>& A, const CVectorRange<NumberS>& b, const VectorRange<NumberS>& x, NumberS tolerance = 1.0e-6f, int iterations = 100);
MYRAMATH_EXPORT stationary_output<NumberD> jacobi(const CSparseMatrixRange<NumberD>& A, const CVectorRange<NumberD>& b, const VectorRange<NumberD>& x, NumberD tolerance = 1.0e-6 , int iterations = 100);
MYRAMATH_EXPORT stationary_output<NumberS> jacobi(const CSparseMatrixRange<NumberC>& A, const CVectorRange<NumberC>& b, const VectorRange<NumberC>& x, NumberS tolerance = 1.0e-6f, int iterations = 100);
MYRAMATH_EXPORT stationary_output<NumberD> jacobi(const CSparseMatrixRange<NumberZ>& A, const CVectorRange<NumberZ>& b, const VectorRange<NumberZ>& x, NumberD tolerance = 1.0e-6 , int iterations = 100);

// Note that seidel() requires the transpose At, because of how it updates unknowns in a particular order.
MYRAMATH_EXPORT stationary_output<NumberS> seidel(const CSparseMatrixRange<NumberS>& At, const CVectorRange<NumberS>& b, const VectorRange<NumberS>& x, NumberS tolerance = 1.0e-6f, int iterations = 100);
MYRAMATH_EXPORT stationary_output<NumberD> seidel(const CSparseMatrixRange<NumberD>& At, const CVectorRange<NumberD>& b, const VectorRange<NumberD>& x, NumberD tolerance = 1.0e-6 , int iterations = 100);
MYRAMATH_EXPORT stationary_output<NumberS> seidel(const CSparseMatrixRange<NumberC>& At, const CVectorRange<NumberC>& b, const VectorRange<NumberC>& x, NumberS tolerance = 1.0e-6f, int iterations = 100);
MYRAMATH_EXPORT stationary_output<NumberD> seidel(const CSparseMatrixRange<NumberZ>& At, const CVectorRange<NumberZ>& b, const VectorRange<NumberZ>& x, NumberD tolerance = 1.0e-6 , int iterations = 100);

// Solves A*x=b using symmetric Seidel iteration. Vector x contains an initial guess on input, overwritten with the solution on output.
// Note that seidel() requires the transpose At, because of how it updates unknowns in a particular order.
MYRAMATH_EXPORT stationary_output<NumberS> sseidel(const CSparseMatrixRange<NumberS>& At, const CVectorRange<NumberS>& b, const VectorRange<NumberS>& x, NumberS tolerance = 1.0e-6f, int iterations = 100);
MYRAMATH_EXPORT stationary_output<NumberD> sseidel(const CSparseMatrixRange<NumberD>& At, const CVectorRange<NumberD>& b, const VectorRange<NumberD>& x, NumberD tolerance = 1.0e-6 , int iterations = 100);
MYRAMATH_EXPORT stationary_output<NumberS> sseidel(const CSparseMatrixRange<NumberC>& At, const CVectorRange<NumberC>& b, const VectorRange<NumberC>& x, NumberS tolerance = 1.0e-6f, int iterations = 100);
MYRAMATH_EXPORT stationary_output<NumberD> sseidel(const CSparseMatrixRange<NumberZ>& At, const CVectorRange<NumberZ>& b, const VectorRange<NumberZ>& x, NumberD tolerance = 1.0e-6 , int iterations = 100);

// Note that sor() requires the transpose At, because of how it updates unknowns in a particular order.
// For best results, the relaxation parameter omega should probably be between 1 and 2. Note that 1 corresponds to seidel().
MYRAMATH_EXPORT stationary_output<NumberS> sor(const CSparseMatrixRange<NumberS>& At, const CVectorRange<NumberS>& b, const VectorRange<NumberS>& x, NumberS omega = 1.5f, NumberS tolerance = 1.0e-6f, int iterations = 100);
MYRAMATH_EXPORT stationary_output<NumberD> sor(const CSparseMatrixRange<NumberD>& At, const CVectorRange<NumberD>& b, const VectorRange<NumberD>& x, NumberD omega = 1.5 , NumberD tolerance = 1.0e-6 , int iterations = 100);
MYRAMATH_EXPORT stationary_output<NumberS> sor(const CSparseMatrixRange<NumberC>& At, const CVectorRange<NumberC>& b, const VectorRange<NumberC>& x, NumberS omega = 1.5 , NumberS tolerance = 1.0e-6f, int iterations = 100);
MYRAMATH_EXPORT stationary_output<NumberD> sor(const CSparseMatrixRange<NumberZ>& At, const CVectorRange<NumberZ>& b, const VectorRange<NumberZ>& x, NumberD omega = 1.5f, NumberD tolerance = 1.0e-6 , int iterations = 100);

// Note that ssor() requires the transpose At, because of how it updates unknowns in a particular order.
// For best results, the relaxation parameter omega should probably be between 1 and 2. Note that 1 corresponsds to sseidel().
MYRAMATH_EXPORT stationary_output<NumberS> ssor(const CSparseMatrixRange<NumberS>& At, const CVectorRange<NumberS>& b, const VectorRange<NumberS>& x, NumberS omega = 1.5f, NumberS tolerance = 1.0e-6f, int iterations = 100);
MYRAMATH_EXPORT stationary_output<NumberD> ssor(const CSparseMatrixRange<NumberD>& At, const CVectorRange<NumberD>& b, const VectorRange<NumberD>& x, NumberD omega = 1.5 , NumberD tolerance = 1.0e-6 , int iterations = 100);
MYRAMATH_EXPORT stationary_output<NumberS> ssor(const CSparseMatrixRange<NumberC>& At, const CVectorRange<NumberC>& b, const VectorRange<NumberC>& x, NumberS omega = 1.5 , NumberS tolerance = 1.0e-6f, int iterations = 100);
MYRAMATH_EXPORT stationary_output<NumberD> ssor(const CSparseMatrixRange<NumberZ>& At, const CVectorRange<NumberZ>& b, const VectorRange<NumberZ>& x, NumberD omega = 1.5f, NumberD tolerance = 1.0e-6 , int iterations = 100);

} // namespace

#endif