lopcg

Source: tests/iterative/lopcg.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/Matrix.h>
#include <myramath/dense/Vector.h>
#include <myramath/dense/DiagonalMatrix.h>
#include <myramath/dense/MatrixRange.h>
#include <myramath/dense/VectorRange.h>
#include <myramath/dense/DiagonalMatrixRange.h>
#include <myramath/sparse/SparseMatrix.h>
#include <myramath/sparse/SparseMatrixRange.h>
#include <myramath/iterative/ICholeskySolver.h>

// Algorithms.
#include <myramath/dense/gemm.h>
#include <myramath/dense/dimm.h>
#include <myramath/dense/frobenius.h>
#include <myramath/dense/euclidean.h>
#include <myramath/sparse/gemm.h>
#include <myramath/sparse/gemv.h>
#include <myramath/sparse/frobenius.h>
#include <myramath/sparse/laplacian2.h>

// Iterative solve/action stuff.
#include <myramath/iterative/Action.h>
#include <myramath/iterative/SolveAction.h>
#include <myramath/iterative/GemmAction.h>
#include <myramath/iterative/lopcg1.h>
#include <myramath/iterative/lopcgn.h>

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

using namespace myra;
using namespace myra_stlprint;

namespace {

template<class Precision> void test(int X, int Y)
  {
  myra::out() << typestring<Precision>() << std::endl;
  int N = X*Y;
  // Construct laplacian, shift it by 10^-3 to make it positive.
  auto A = laplacian2<Precision>(X,Y);
  Precision one(1);
  Precision ten(10);
  Precision sigma = std::pow(ten,-3);
  for (int n = 0; n < N; ++n)
    A(n,n) = A(n,n) - one + sigma;
  A *= one /frobenius(A);
  // Form action for A.
  auto A_action = make_GemmAction(A);
  // Form incomplete cholesky preconditioner, M.
  ICholeskySolver<Precision> M(A);
  auto M_action = make_SolveAction(M);
  // Form the (unwanted) constant eigenvector, will be deflated out.
  auto Q = Matrix<Precision>::ones(N,1);
  normalize(Q.range().vector(0));
  // Generate one eigenpair.
  auto v_mu = lopcg1(M_action,A_action,Q);
  Vector<Precision>& v = v_mu.first;
  Precision lambda = Precision(1)/v_mu.second;
  // Check residual, A*v - v*lambda
  Vector<Precision> r = A*v-v*lambda;
  Precision r_error = frobenius(r);
  myra::out() << "  |A*v - v*lambda| = " << r_error << std::endl;
  REQUIRE(r_error < 1.0e-3);
  // Generate 3 eigenpairs
  int size = 3;
  auto V_Mu = lopcgN(M_action,A_action,Q,size);
  Matrix<Precision>& V = V_Mu.first;
  DiagonalMatrix<Precision> Lambda = V_Mu.second;
  for (int ij = 0; ij < size; ++ij)
    Lambda(ij) = Precision(1) / Lambda(ij);
  // Check residual, A*V - V*Lambda
  Matrix<Precision> R = A*V - V*Lambda;
  Precision R_error = frobenius(R);
  // Check orthogonality V'V-I
  auto I = Matrix<Precision>::identity(size);
  Precision I_error = frobenius(gemm(V,'T',V,'N')-I);  
  myra::out() << "  |A*V - V*Lambda| = " << R_error << std::endl;
  myra::out() << "  |V'V - I| = " << R_error << std::endl;  
  myra::out() << "  Mu = " << V_Mu.second.range() << std::endl;
  REQUIRE(R_error < 1.0e-3);
  REQUIRE(I_error < 1.0e-3);  
  }

} // namespace

ADD_TEST("lopcg","[iterative]")
  {
  // Test on square laplacian (degenerate, then 2:1 gap)
  test<double>(32,32);
  test<float >(32,32);
  // Test on a squarish laplacian (near degenerate, then 2:1 gap)  
  test<double>(32,31);
  test<double>(32,16);
  // Test on a 2:1 laplacian (looking for 4:1 gap)
  test<float >(32,31);
  test<float >(32,16);
  }

Results: [PASS]

double
  |A*v - v*lambda| = 2.82036e-08
  |A*V - V*Lambda| = 5.71556e-08
  |V'V - I| = 5.71556e-08
  Mu = size 3 DiagonalMatrix of double: [ 13126.5 13126.5 6887.18 ]

float
  |A*v - v*lambda| = 6.89423e-05
  |A*V - V*Lambda| = 0.000183611
  |V'V - I| = 0.000183611
  Mu = size 3 DiagonalMatrix of float: [ 13124.8 9546.3 6899.35 ]

double
  |A*v - v*lambda| = 1.68718e-06
  |A*V - V*Lambda| = 6.89533e-08
  |V'V - I| = 6.89533e-08
  Mu = size 3 DiagonalMatrix of double: [ 12914.4 12191 6571.3 ]

double
  |A*v - v*lambda| = 4.84616e-08
  |A*V - V*Lambda| = 1.53363e-07
  |V'V - I| = 1.53363e-07
  Mu = size 3 DiagonalMatrix of double: [ 9160.9 2469.86 2469.86 ]

float
  |A*v - v*lambda| = 3.80715e-05
  |A*V - V*Lambda| = 0.000146797
  |V'V - I| = 0.000146797
  Mu = size 3 DiagonalMatrix of float: [ 11584.7 10065.8 7610.83 ]

float
  |A*v - v*lambda| = 4.17571e-06
  |A*V - V*Lambda| = 1.36859e-05
  |V'V - I| = 1.36859e-05
  Mu = size 3 DiagonalMatrix of float: [ 9160.15 2469.8 2469.76 ]

Go back to Summary of /test programs.