stationary

Source: tests/iterative/stationary.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>

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

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

using namespace myra;

namespace {

template<class Precision> void test(int I, int J, Precision tol)
  {
  // Setup A and b for two-dimensional laplacian operator, div(grad(phi)) = 0 ..
  int N = I*J;
  auto A = laplacian2<Precision>(I,J);
  auto b = Vector<Precision>::zeros(I*J);
  // .. subjecting A and b to phi=1 boundary conditions
  Precision one(1);
  Natural2D ordering(I,J);
  for (int n = 0; n < I*J; ++n)
    A(n,n) -= one;
  // .. along j=0 and j=1 walls
  for (int i = 0; i < I; ++i)
    {
    int n0 = ordering(i,0);
    b(n0) += one;
    A(n0,n0) += one;
    int n1 = ordering(i,J-1);
    b(n1) += one;
    A(n1,n1) += one;
    }
  // .. along i=0 and i=1 walls    
  for (int j = 0; j < J; ++j)
    {
    int n0 = ordering(0,j);
    b(n0) += one;
    A(n0,n0) += one;
    int n1 = ordering(I-1,j);
    b(n1) += one;
    A(n1,n1) += one;
    }

  // Verify A is symmetric.
  myra::out() << "|A-A'| = " << frobenius(A-transpose(A)) << std::endl;

  // Verify exact solution, A*ones = b
  auto y = Vector<Precision>::ones(N);
  myra::out() << "|A*ones-b| = " << euclidean(A*y-b) << std::endl;

  // Solve A*x=b using jacobi()  
  auto x1 = Vector<Precision>::zeros(N);
  auto output1 = jacobi(A,b,x1,tol,1000);
  Precision error1 = euclidean(A*x1-b) / euclidean(b);
  myra::out() << "jacobi |A*x-b|/|b| = " << error1 << ", iterations = " << output1.history.size() << std::endl;
  REQUIRE(error1 < tol);

  // Solve A*x=b using seidel()
  auto x2 = Vector<Precision>::zeros(N);
  auto output2 = seidel(A,b,x2,tol,1000);
  Precision error2 = euclidean(A*x2-b) / euclidean(b);
  myra::out() << "seidel |A*x-b|/|b| = " << error2 << ", iterations = " << output2.history.size() << std::endl;
  REQUIRE(error2 < tol);

  // Solve A*x=b using seidel()
  auto x3 = Vector<Precision>::zeros(N);
  auto output3 = sseidel(A,b,x3,tol,1000);
  Precision error3 = euclidean(A*x3-b) / euclidean(b);
  myra::out() << "sseidel |A*x-b|/|b| = " << error3 << ", iterations = " << output3.history.size() << std::endl;
  REQUIRE(error3 < tol);

  // Solve A*x = b using sor()
  Precision omega(1.5);
  auto x4 = Vector<Precision>::zeros(N);
  auto output4 = sor(A,b,x4,omega,tol,1000);
  Precision error4 = euclidean(A*x4-b) / euclidean(b);
  myra::out() << "sor |A*x-b|/|b| = " << error4 << ", iterations = " << output4.history.size() << std::endl;
  REQUIRE(error4 < tol);

  // Solve A*x = b using ssor()
  auto x5 = Vector<Precision>::zeros(N);
  auto output5 = ssor(A,b,x5,omega,tol,1000);
  Precision error5 = euclidean(A*x5-b) / euclidean(b);
  myra::out() << "ssor |A*x-b|/|b| = " << error5 << ", iterations = " << output5.history.size() << std::endl;
  REQUIRE(error5 < tol);

  }

} // namespace

ADD_TEST("stationary","[iterative]")
  {
  test<float>(20,20,1.0e-4f);
  test<double>(20,20,1.0e-4);
  }

Results: [PASS]

|A-A'| = 0
|A*ones-b| = 0
jacobi |A*x-b|/|b| = 9.99893e-05, iterations = 597
seidel |A*x-b|/|b| = 9.81902e-05, iterations = 300
sseidel |A*x-b|/|b| = 6.78522e-05, iterations = 160
sor |A*x-b|/|b| = 9.50828e-05, iterations = 99
ssor |A*x-b|/|b| = 8.15185e-05, iterations = 56
|A-A'| = 0
|A*ones-b| = 0
jacobi |A*x-b|/|b| = 9.99856e-05, iterations = 597
seidel |A*x-b|/|b| = 9.819e-05, iterations = 300
sseidel |A*x-b|/|b| = 6.78542e-05, iterations = 160
sor |A*x-b|/|b| = 9.50743e-05, iterations = 99
ssor |A*x-b|/|b| = 8.15192e-05, iterations = 56

Go back to Summary of /test programs.