compose_action2

Source: tests/iterative/compose2.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 for A,P and and inv(A)
#include <myramath/utility/Number.h>
#include <myramath/dense/Matrix.h>
#include <myramath/dense/MatrixRange.h>

// Container for x/b
#include <myramath/dense/Vector.h>
#include <myramath/dense/VectorRange.h>

// For examining eigenvalues of inv(A)*(A+P*P')
#include <myramath/dense/geevd.h>
#include <myramath/dense/DiagonalMatrix.h>
#include <myramath/dense/DiagonalMatrixRange.h>
#include <myramath/dense/euclidean.h>
#include <myramath/dense/expr.h>
#include <myramath/expression/functions.h>
#include <myramath/expression/print.h>

// An underlying direct solver for A, wrapped by make_SolveAction.
#include <myramath/dense/LUSolver.h>

// Iterative solve/action stuff.
#include <myramath/iterative/gmres.h>
#include <myramath/iterative/Action.h>
#include <myramath/iterative/SumAction.h>
#include <myramath/iterative/ProductAction.h>
#include <myramath/iterative/GemmAction.h>
#include <myramath/iterative/SolveAction.h>

// Reporting.
#include <myramath/utility/stlprint.h>
#include <tests/myratest.h>

using namespace myra;
using namespace myra_stlprint;

namespace {

// Tests Action composition using gmres()
template<class Number> void test(int N, int K, typename ReflectPrecision<Number>::type tolerance)
  {
  typedef typename ReflectPrecision<Number>::type Precision;
  myra::out() << typestring<Number>() << std::endl;
  // Form random A and a direct solver (lu) for it.
  auto A = Matrix<Number>::random(N,N);
  LUSolver<Number> solver(A.add_const());
  // Form forward action (A+PP', a rank-K perturbation of A).
  auto P = Matrix<Number>::random(N,K);
  auto forward = make_GemmAction(A) + make_GemmAction(P)*make_GemmAction(P,'T');
  auto preconditioner = make_SolveAction(solver);
  // Examine eigenvalues of preconditioner*forward, should be mostly 1's with K non-1's
  /*
  auto DRL = geevd( (preconditioner*forward).make_Matrix() );
  myra::out() << "|lambda| = " << abs(expr1(DRL.D)) << std::endl;
  */
  // Solve forward*x = b, using gmres(preconditioner,forward)
  auto b = Vector<Number>::random(N);
  auto x = Vector<Number>::zeros(N);
  auto result = gmres(preconditioner,forward,b,x,tolerance);
  // Check residual is small.
  Precision residual = euclidean(forward*x-b) / euclidean(b);
  myra::out() << "residual = " << residual << std::endl;
  // Since (preconditioner*action) is a rank-K perturbation of identity, gmres() should converge in O(K) iterations.
  myra::out() << "history = " << result.history << std::endl;
  // Perform assertions.  
  REQUIRE(residual < tolerance*100.0);
  REQUIRE(result.history.size() <= K+5);
  }

} // namespace

ADD_TEST("compose_action2","[iterative]")
  {
  // Test gmres() on a nonsymmetric composition.
  test<NumberD>(100,5,1.0e-8);
  test<NumberZ>(100,5,1.0e-8);
  }
  

Results: [PASS]

double
residual = 8.17273e-14
history = [ 47.2166 7.17701 2.92443 0.855566 0.196902 0.196818 7.43497e-14 ] (7)
std::complex<double>
residual = 1.96979e-13
history = [ 35.5167 1.26435 1.11192 1.02587 0.637327 0.545325 1.1242e-13 ] (7)

Go back to Summary of /test programs.