sseidel_symmetry

Source: tests/iterative/stationary_symmetry.cpp

// ========================================================================= //
// 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.      //
// ========================================================================= //

// Containers.
#include <myramath/utility/Number.h>
#include <myramath/dense/Vector.h>
#include <myramath/dense/VectorRange.h>
#include <myramath/dense/Matrix.h>
#include <myramath/dense/MatrixRange.h>
#include <myramath/sparse/SparseMatrix.h>
#include <myramath/sparse/SparseMatrixRange.h>
#include <myramath/iterative/Action.h>
#include <myramath/iterative/UserAction.h>

// Algorithms.
#include <myramath/dense/frobenius.h>
#include <myramath/dense/transpose.h>
#include <myramath/sparse/laplacian2.h>
#include <myramath/iterative/stationary.h>

// Reporting.
#include <tests/myratest.h>

using namespace myra;

namespace {

// Action that applies seidel()
template<class Precision> class seidel_Action
  {
  public:

    // Required typedef.
    typedef Precision Number;

    seidel_Action(const CSparseMatrixRange<Precision>& in_At, Precision in_tol, int in_iter)
      : At(in_At), tol(in_tol), iter(in_iter) { }

    std::pair<int,int> size() const
      { return At.size(); }

    void multiply (const CMatrixRange<Precision>& B, const MatrixRange<Precision>& X) const
      { 
      X.zero();
      for (int j = 0; j < B.J; ++j)
        seidel(At,B.vector(j),X.vector(j),tol,iter); 
      }

  private:

    // Captured environment for seidel()
    CSparseMatrixRange<Precision> At;
    Precision tol;
    int iter;
    
  };

// Action that applies sseidel()
template<class Precision> class sseidel_Action
  {
  public:

    // Required typedef.
    typedef Precision Number;

    sseidel_Action(const CSparseMatrixRange<Precision>& in_At, Precision in_tol, int in_iter)
      : At(in_At), tol(in_tol), iter(in_iter) { }

    std::pair<int,int> size() const
      { return At.size(); }

    void multiply (const CMatrixRange<Precision>& B, const MatrixRange<Precision>& X) const
      { 
      X.zero();
      for (int j = 0; j < B.J; ++j)
        sseidel(At,B.vector(j),X.vector(j),tol,iter); 
      }

  private:

    // Captured environment for sseidel()
    CSparseMatrixRange<Precision> At;
    Precision tol;
    int iter;
    
  };

// Tests seidel() and sseidel()
template<class Precision> void test1(int I, int J, Precision tol)
  {
  auto A = laplacian2<Precision>(I,J);
  // Look at the assymetry of seidel()
  Action<Precision> S = make_UserAction( seidel_Action<Precision>(A,0.0,10) );
  Matrix<Precision> S0 = S.make_Matrix();
  Precision seidel_assymetry = frobenius(S0-transpose(S0));
  myra::out() << "| seidel -  seidel'| = " << seidel_assymetry << std::endl;
  // Verify that sseidel() is a symmetric Action.
  Action<Precision> SS = make_UserAction( sseidel_Action<Precision>(A,0.0,10) );
  Matrix<Precision> SS0 = SS.make_Matrix();
  Precision sseidel_assymetry = frobenius(SS0-transpose(SS0));
  myra::out() << "|sseidel - sseidel'| = " << sseidel_assymetry << std::endl;
  REQUIRE( sseidel_assymetry < tol );
  }

// Action that applies sor()
template<class Precision> class sor_Action
  {
  public:

    // Required typedef.
    typedef Precision Number;

    sor_Action(const CSparseMatrixRange<Precision>& in_At, Precision in_omega, Precision in_tol, int in_iter)
      : At(in_At), omega(in_omega), tol(in_tol), iter(in_iter) { }

    std::pair<int,int> size() const
      { return At.size(); }

    void multiply (const CMatrixRange<Precision>& B, const MatrixRange<Precision>& X) const
      { 
      X.zero();
      for (int j = 0; j < B.J; ++j)
        sor(At,B.vector(j),X.vector(j),omega,tol,iter); 
      }

  private:

    // Captured environment for sor()
    CSparseMatrixRange<Precision> At;
    Precision omega;
    Precision tol;
    int iter;
    
  };

// Action that applies ssor()
template<class Precision> class ssor_Action
  {
  public:

    // Required typedef.
    typedef Precision Number;

    ssor_Action(const CSparseMatrixRange<Precision>& in_At, Precision in_omega, Precision in_tol, int in_iter)
      : At(in_At), omega(in_omega), tol(in_tol), iter(in_iter) { }

    std::pair<int,int> size() const
      { return At.size(); }

    void multiply (const CMatrixRange<Precision>& B, const MatrixRange<Precision>& X) const
      { 
      X.zero();
      for (int j = 0; j < B.J; ++j)
        ssor(At,B.vector(j),X.vector(j),omega,tol,iter); 
      }

  private:

    // Captured environment for ssor()
    CSparseMatrixRange<Precision> At;
    Precision omega;
    Precision tol;
    int iter;
    
  };

// Tests sor() and ssor()
template<class Precision> void test2(int I, int J, Precision tol)
  {
  auto A = laplacian2<Precision>(I,J);
  Precision omega(1.5);
  // Look at the assymetry of sor()
  Action<Precision> S = make_UserAction( sor_Action<Precision>(A,omega,0.0,10) );
  Matrix<Precision> S0 = S.make_Matrix();
  Precision sor_assymetry = frobenius(S0-transpose(S0));
  myra::out() << "| sor -  sor'| = " << sor_assymetry << std::endl;
  // Verify that ssor() is a symmetric Action.
  Action<Precision> SS = make_UserAction( ssor_Action<Precision>(A,omega,0.0,10) );
  Matrix<Precision> SS0 = SS.make_Matrix();
  Precision ssor_assymetry = frobenius(SS0-transpose(SS0));
  myra::out() << "|ssor - ssor'| = " << ssor_assymetry << std::endl;
  REQUIRE( ssor_assymetry < tol );
  }


} // namespace

ADD_TEST("sseidel_symmetry","[iterative]")
  {
  test1<float>(20,20,1.0e-5f);
  test1<double>(20,20,1.0e-12);
  }

ADD_TEST("ssor_symmetry","[iterative]")
  {
  test2<float>(20,20,1.0e-5f);
  test2<double>(20,20,1.0e-12);
  }

Results: [PASS]

| seidel -  seidel'| = 0.0376709
|sseidel - sseidel'| = 3.23786e-07
| seidel -  seidel'| = 0.0376727
|sseidel - sseidel'| = 6.37134e-16

Go back to Summary of /test programs.