conjugate_gradient.cc File Reference

#include <random>
#include "catch2/catch_all.hpp"
#include "sopt/conjugate_gradient.h"

Include dependency graph for conjugate_gradient.cc:

Go to the source code of this file.

Functions
	TEST_CASE ("Conjugate gradient", "[cg]")

Function Documentation

TEST_CASE()

TEST_CASE	(	"Conjugate gradient"	,
		""	[cg]
	)

Definition at line 6 of file conjugate_gradient.cc.

                                        {
  using namespace sopt;
 
  ConjugateGradient const cg(std::numeric_limits<t_uint>::max(), 1e-12);
  SECTION("Real valued") {
    auto const A = Image<>::Random(10, 10).eval();
    auto const AtA = (A.transpose().matrix() * A.matrix()).eval();
    auto const expected = Array<>::Random(A.rows()).eval();
 
    auto const actual = cg(AtA, (A.transpose().matrix() * expected.matrix()).eval());
 
    CHECK(actual.niters > 0);
    CHECK(std::abs(actual.residual) < 1e-6);
    CAPTURE(actual.residual);
    CAPTURE((A.matrix() * actual.result).transpose());
    CAPTURE(expected.transpose());
    CHECK((A.matrix() * actual.result).isApprox(expected.matrix(), 1e-6));
  }
 
  SECTION("Complex valued") {
    auto const A = Image<t_complex>::Random(10, 10).eval();
    auto const AhA = (A.conjugate().transpose().matrix() * A.matrix()).eval();
    auto const expected = Array<t_complex>::Random(A.rows()).eval();
 
    auto const actual = cg(AhA, (A.conjugate().transpose().matrix() * expected.matrix()).eval());
 
    CHECK(actual.niters > 0);
    CHECK(std::abs(actual.residual) < 1e-6);
    CAPTURE(actual.residual);
    CAPTURE((A.matrix() * actual.result).transpose());
    CAPTURE(expected.transpose());
    CHECK((A.matrix() * actual.result).isApprox(expected.matrix(), 1e-6));
  }
}