Quellcodebibliothek Statistik Leitseite products/sources/formale Sprachen/C/MySQL/unsupported/test/   (MySQL Server Version 8.1-8.4©)  Datei vom 12.11.2025 mit Größe 7 kB image not shown  

Quelle  matrix_function.cpp   Sprache: C

 
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

#include "main.h"
#include <unsupported/Eigen/MatrixFunctions>

// Variant of VERIFY_IS_APPROX which uses absolute error instead of
// relative error.
#define VERIFY_IS_APPROX_ABS(a, b) VERIFY(test_isApprox_abs(a, b))

template<typename Type1, typename Type2>
inline bool test_isApprox_abs(const Type1& a, const Type2& b)
{
  return ((a-b).array().abs() < test_precision<typename Type1::RealScalar>()).all();
}


// Returns a matrix with eigenvalues clustered around 0, 1 and 2.
template<typename MatrixType>
MatrixType randomMatrixWithRealEivals(const Index size)
{
  typedef typename MatrixType::Scalar Scalar;
  typedef typename MatrixType::RealScalar RealScalar;
  MatrixType diag = MatrixType::Zero(size, size);
  for (Index i = 0; i < size; ++i) {
    diag(i, i) = Scalar(RealScalar(internal::random<int>(0,2)))
      + internal::random<Scalar>() * Scalar(RealScalar(0.01));
  }
  MatrixType A = MatrixType::Random(size, size);
  HouseholderQR<MatrixType> QRofA(A);
  return QRofA.householderQ().inverse() * diag * QRofA.householderQ();
}

template <typename MatrixType, int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex>
struct randomMatrixWithImagEivals
{
  // Returns a matrix with eigenvalues clustered around 0 and +/- i.
  static MatrixType run(const Index size);
};

// Partial specialization for real matrices
template<typename MatrixType>
struct randomMatrixWithImagEivals<MatrixType, 0>
{
  static MatrixType run(const Index size)
  {
    typedef typename MatrixType::Scalar Scalar;
    MatrixType diag = MatrixType::Zero(size, size);
    Index i = 0;
    while (i < size) {
      Index randomInt = internal::random<Index>(-1, 1);
      if (randomInt == 0 || i == size-1) {
        diag(i, i) = internal::random<Scalar>() * Scalar(0.01);
        ++i;
      } else {
        Scalar alpha = Scalar(randomInt) + internal::random<Scalar>() * Scalar(0.01);
        diag(i, i+1) = alpha;
        diag(i+1, i) = -alpha;
        i += 2;
      }
    }
    MatrixType A = MatrixType::Random(size, size);
    HouseholderQR<MatrixType> QRofA(A);
    return QRofA.householderQ().inverse() * diag * QRofA.householderQ();
  }
};

// Partial specialization for complex matrices
template<typename MatrixType>
struct randomMatrixWithImagEivals<MatrixType, 1>
{
  static MatrixType run(const Index size)
  {
    typedef typename MatrixType::Scalar Scalar;
    typedef typename MatrixType::RealScalar RealScalar;
    const Scalar imagUnit(0, 1);
    MatrixType diag = MatrixType::Zero(size, size);
    for (Index i = 0; i < size; ++i) {
      diag(i, i) = Scalar(RealScalar(internal::random<Index>(-1, 1))) * imagUnit
        + internal::random<Scalar>() * Scalar(RealScalar(0.01));
    }
    MatrixType A = MatrixType::Random(size, size);
    HouseholderQR<MatrixType> QRofA(A);
    return QRofA.householderQ().inverse() * diag * QRofA.householderQ();
  }
};


template<typename MatrixType>
void testMatrixExponential(const MatrixType& A)
{
  typedef typename internal::traits<MatrixType>::Scalar Scalar;
  typedef typename NumTraits<Scalar>::Real RealScalar;
  typedef std::complex<RealScalar> ComplexScalar;

  VERIFY_IS_APPROX(A.exp(), A.matrixFunction(internal::stem_function_exp<ComplexScalar>));
}

template<typename MatrixType>
void testMatrixLogarithm(const MatrixType& A)
{
  typedef typename internal::traits<MatrixType>::Scalar Scalar;
  typedef typename NumTraits<Scalar>::Real RealScalar;

  MatrixType scaledA;
  RealScalar maxImagPartOfSpectrum = A.eigenvalues().imag().cwiseAbs().maxCoeff();
  if (maxImagPartOfSpectrum >= RealScalar(0.9L * EIGEN_PI))
    scaledA = A * RealScalar(0.9L * EIGEN_PI) / maxImagPartOfSpectrum;
  else
    scaledA = A;

  // identity X.exp().log() = X only holds if Im(lambda) < pi for all eigenvalues of X
  MatrixType expA = scaledA.exp();
  MatrixType logExpA = expA.log();
  VERIFY_IS_APPROX(logExpA, scaledA);
}

template<typename MatrixType>
void testHyperbolicFunctions(const MatrixType& A)
{
  // Need to use absolute error because of possible cancellation when
  // adding/subtracting expA and expmA.
  VERIFY_IS_APPROX_ABS(A.sinh(), (A.exp() - (-A).exp()) / 2);
  VERIFY_IS_APPROX_ABS(A.cosh(), (A.exp() + (-A).exp()) / 2);
}

template<typename MatrixType>
void testGonioFunctions(const MatrixType& A)
{
  typedef typename MatrixType::Scalar Scalar;
  typedef typename NumTraits<Scalar>::Real RealScalar;
  typedef std::complex<RealScalar> ComplexScalar;
  typedef Matrix<ComplexScalar, MatrixType::RowsAtCompileTime, 
                 MatrixType::ColsAtCompileTime, MatrixType::Options> ComplexMatrix;

  ComplexScalar imagUnit(0,1);
  ComplexScalar two(2,0);

  ComplexMatrix Ac = A.template cast<ComplexScalar>();
  
  ComplexMatrix exp_iA = (imagUnit * Ac).exp();
  ComplexMatrix exp_miA = (-imagUnit * Ac).exp();
  
  ComplexMatrix sinAc = A.sin().template cast<ComplexScalar>();
  VERIFY_IS_APPROX_ABS(sinAc, (exp_iA - exp_miA) / (two*imagUnit));
  
  ComplexMatrix cosAc = A.cos().template cast<ComplexScalar>();
  VERIFY_IS_APPROX_ABS(cosAc, (exp_iA + exp_miA) / 2);
}

template<typename MatrixType>
void testMatrix(const MatrixType& A)
{
  testMatrixExponential(A);
  testMatrixLogarithm(A);
  testHyperbolicFunctions(A);
  testGonioFunctions(A);
}

template<typename MatrixType>
void testMatrixType(const MatrixType& m)
{
  // Matrices with clustered eigenvalue lead to different code paths
  // in MatrixFunction.h and are thus useful for testing.

  const Index size = m.rows();
  for (int i = 0; i < g_repeat; i++) {
    testMatrix(MatrixType::Random(size, size).eval());
    testMatrix(randomMatrixWithRealEivals<MatrixType>(size));
    testMatrix(randomMatrixWithImagEivals<MatrixType>::run(size));
  }
}

template<typename MatrixType>
void testMapRef(const MatrixType& A)
{
  // Test if passing Ref and Map objects is possible
  // (Regression test for Bug #1796)
  Index size = A.rows();
  MatrixType X; X.setRandom(size, size);
  MatrixType Y(size,size);
  Ref<      MatrixType> R(Y);
  Ref<const MatrixType> Rc(X);
  Map<      MatrixType> M(Y.data(), size, size);
  Map<const MatrixType> Mc(X.data(), size, size);

  X = X*X; // make sure sqrt is possible
  Y = X.sqrt();
  R = Rc.sqrt();
  M = Mc.sqrt();
  Y = X.exp();
  R = Rc.exp();
  M = Mc.exp();
  X = Y; // make sure log is possible
  Y = X.log();
  R = Rc.log();
  M = Mc.log();

  Y = X.cos() + Rc.cos() + Mc.cos();
  Y = X.sin() + Rc.sin() + Mc.sin();

  Y = X.cosh() + Rc.cosh() + Mc.cosh();
  Y = X.sinh() + Rc.sinh() + Mc.sinh();
}


EIGEN_DECLARE_TEST(matrix_function)
{
  CALL_SUBTEST_1(testMatrixType(Matrix<float,1,1>()));
  CALL_SUBTEST_2(testMatrixType(Matrix3cf()));
  CALL_SUBTEST_3(testMatrixType(MatrixXf(8,8)));
  CALL_SUBTEST_4(testMatrixType(Matrix2d()));
  CALL_SUBTEST_5(testMatrixType(Matrix<double,5,5,RowMajor>()));
  CALL_SUBTEST_6(testMatrixType(Matrix4cd()));
  CALL_SUBTEST_7(testMatrixType(MatrixXd(13,13)));

  CALL_SUBTEST_1(testMapRef(Matrix<float,1,1>()));
  CALL_SUBTEST_2(testMapRef(Matrix3cf()));
  CALL_SUBTEST_3(testMapRef(MatrixXf(8,8)));
  CALL_SUBTEST_7(testMapRef(MatrixXd(13,13)));
}

93%


¤ Dauer der Verarbeitung: 0.1 Sekunden  (vorverarbeitet)  ¤

*© Formatika GbR, Deutschland






Wurzel

Suchen

Beweissystem der NASA

Beweissystem Isabelle

NIST Cobol Testsuite

Cephes Mathematical Library

Wiener Entwicklungsmethode

Haftungshinweis

Die Informationen auf dieser Webseite wurden nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit, noch Qualität der bereit gestellten Informationen zugesichert.

Bemerkung:

Die farbliche Syntaxdarstellung ist noch experimentell.