algo_test.h

一个通用的数学库

#ifndef MTL_ALGO_TEST_H
#define MTL_ALGO_TEST_H

#include "mtl/mtl.h"
#include "mtl/dense1D.h"
#include "fill_matrix.h"

using mtl::dense1D;

template <class Matrix>
bool
mat_algo_test(Matrix& A, std::string test_name)
{
  typedef typename mtl::matrix_traits<Matrix>::value_type T;
  typedef typename mtl::matrix_traits<Matrix>::size_type Int;
  typename Matrix::iterator i;
  typename Matrix::OneD::iterator j;
  T s;
  T sum, norm;
  bool success = true;

  mtl::insert_zero_matrix(A);

  mtl::matrix<T>::type C(A.nrows(), A.ncols());
  mtl::zero_matrix(C);

  // set
  mtl::set_value(A, T(5));

  for (i = A.begin(); i != A.end(); ++i)
    for (j = (*i).begin(); j != (*i).end(); ++j)
      if (*j != T(5)) {
        cerr << "**** FAILED: (matrix set) " << test_name.c_str() << " ****" << endl;
        success = false;
        break;
      }

  // scale
  mtl::scale(A, T(5));

  for (i = A.begin(); i != A.end(); ++i)
    for (j = (*i).begin(); j != (*i).end(); ++j)
      if (*j != T(25)) {
        cerr << "**** FAILED: (matrix set) " << test_name.c_str() << " ****" << endl;
        success = false;      
        break;
      }

  matrix_fill(A);

  mtl::copy(A, C);

  // one_norm
  {
    Int i, j;
    s = mtl::one_norm(A);

    if (C.ncols() > 0) {
      i = 0; sum = T();
      for (j = 0; j < Int(C.nrows()); ++j)
        sum += MTL_ABS(C(j,i));
      norm = sum;
      ++i;
    }
    
    for (; i < Int(C.ncols()); ++i) {
      sum = T();
      for (j = 0; j < Int(C.nrows()); ++j)
        sum += MTL_ABS(C(j,i));
      norm = MTL_MAX(MTL_ABS(norm), MTL_ABS(sum));
    }
    
    if (s != norm) {
      cout << test_name.c_str() << " failed one_norm(A) ****" << endl;
      success = false;      
    }
  }

  {
    // infinity_norm
    Int i, j;
    s = mtl::infinity_norm(A);
    
    if (C.nrows() > 0) {
      i = 0; sum = T();
      for (j = 0; j < Int(C.ncols()); ++j)
        sum += MTL_ABS(C(i,j));
      norm = sum;
      ++i;
    }
    
    for (; i < Int(C.nrows()); ++i) {
      sum = T();
      for (j = 0; j < Int(C.ncols()); ++j)
        sum += MTL_ABS(C(i,j));
      norm = MTL_MAX(MTL_ABS(norm), MTL_ABS(sum));
    }
    
    if (s != norm) {
      cout << test_name.c_str() << " failed infinity_norm(A) ****" << endl;
      success = false;      
    }
  }

  set_diagonal_test(test_name, A, success);

  if (success)
    cout << test_name.c_str() << " passed single matrix algorithms test" << endl;

  return success;
}

template <class Matrix>
void
set_diagonal_test(std::string test_name, Matrix& A, bool& success)
{
  if (A.is_unit()) {
    cout << test_name.c_str() << " skipping set_diagonal" << endl;
    return;
  }
  typedef typename mtl::matrix_traits<Matrix>::value_type T;
  typedef typename mtl::matrix_traits<Matrix>::size_type Int;
  // set_diagnal
  Int i;
  mtl::set_diagonal(A, T(1));
  
  for (i = 0; i < A.nrows() && i < A.ncols(); ++i)
    if (A(i,i) != T(1)) {
      cout << test_name.c_str() << " failed set_diagnal(A,a) *****" << endl;
      success = false;
      break;
    }
}

template <class Mat>
inline void
fill_matrix(Mat& A, int /*sub*/, int /*super*/, mtl::rectangle_tag)
{
  typedef typename mtl::matrix_traits<Mat>::size_type Int;
  typedef typename mtl::matrix_traits<Mat>::value_type T;
  Int i, j;
  for (i = 0; i < A.nrows(); ++i)
    for (j = 0; j < A.ncols(); ++j)
      A(i,j) = T(i * A.ncols() + j);
}

template <class Mat>
inline void
fill_matrix(Mat& A, int sub, int super, mtl::banded_tag)
{
  typedef typename mtl::matrix_traits<Mat>::size_type Int;
  typedef typename mtl::matrix_traits<Mat>::value_type T;
  Int i, j;
  for (i = 0; i < A.nrows(); ++i) {
    Int first = MTL_MAX(0, int(i) - sub);
    Int last = MTL_MIN(int(A.ncols()), int(i) + super + 1);
    for (j = 0; j < A.ncols(); ++j)
      if (j >= first && j < last)
        A(i,j) = T(i * A.ncols() + j);
  }
}

template <class Mat>
inline void
fill_matrix(Mat& A, int sub, int super, mtl::symmetric_tag)
{
  typedef typename mtl::matrix_traits<Mat>::size_type Int;
  typedef typename mtl::matrix_traits<Mat>::value_type T;
  Int i, j;
  for (i = 0; i < A.nrows(); ++i) {
    Int first = MTL_MAX(0, int(i) - sub);
    Int last = MTL_MIN(int(A.ncols()), int(i) + super + 1);
    for (j = 0; j < A.ncols(); ++j)
      if (j >= first && j < last)
        A(i,j) = T(i + j);
  }
}



template <class Matrix>
void test_matvec_mult(std::string test, Matrix& A)
{

  typedef typename mtl::matrix_traits<Matrix>::value_type T;
  typedef typename mtl::matrix_traits<Matrix>::size_type Int;
  typedef typename mtl::matrix_traits<Matrix>::shape Shape;
  Int M = A.nrows();
  Int N = A.ncols();
  Int i, j;

  dense1D<T> x(N), y(M);

  mtl::matrix<T>::type AA(M, N);
  mtl::dense1D<T> yy(M);
  mtl::dense1D<T> z(M);

  bool pass;

  //
  // mult test    y = A x + y
  //

  mtl::set_value(AA, T());
  mtl::set_value(A, T());

  fill_matrix(A, A.sub(), A.super(), Shape());
  mtl::copy(A, AA);
  //  fill_matrix(AA, A.sub(), A.super(), Shape());

  mtl::set_value(x, T(1));
  mtl::set_value(y, T(1));
  mtl::set_value(yy, T());

  mtl::mult(A, x, y);

  for (i = 0; i < M; ++i)
    for (j = 0; j < N; ++j)
      yy[i] += AA(i,j) * x[j];

  // compare y's
  pass = std::equal(y.begin(), y.end(), yy.begin());

  if (pass)
    cout << test.c_str() << " passed matvec mult1" << endl;
  else {
    cout << "*** matvec mult1 " << test.c_str() << " failed" << endl;
#if !defined(_MSVCPP_)
    cout << "result vector ";
    mtl::print_vector(y.begin(), y.end());
    cout << "correct vector ";
    mtl::print_vector(yy.begin(), yy.end());
#endif
  }

  //
  //  mult2
  //  
  mtl::set_value(y, T(1));
  mtl::set_value(yy, T(1));
  mtl::set_value(z, T(9));
  
  mtl::mult(A, x, y, z);

  for (i = 0; i < M; ++i) {
    T tmp = y[i];
    for (j = 0; j < N; ++j)
      tmp += AA(i,j) * x[j];
    yy[i] = tmp;
  }
  pass = std::equal(z.begin(), z.end(), yy.begin());

  if (pass)
    cout << test.c_str() << " passed matvec mult2 " << endl;
  else {
    cout << "*** matvec mult2 " << test.c_str() << " failed" << endl;
#if !defined(_MSVCPP_)
    cout << "result vector ";
    mtl::print_vector(z.begin(), z.end());
    cout << "correct vector ";
    mtl::print_vector(yy.begin(), yy.end());
#endif
  }
}

template <class Matrix>
void test_matvec_rankone(std::string test, Matrix&, banded_tag)
{
  cout << test.c_str() << " skipping rank one update" << endl;
}


template <class Matrix>
void test_matvec_rankone(std::string test, Matrix& A, rectangle_tag)
{
  typedef typename mtl::matrix_traits<Matrix>::value_type T;
  typedef typename mtl::matrix_traits<Matrix>::size_type Int;
  typedef typename mtl::matrix_traits<Matrix>::shape Shape;
  Int M = A.nrows();
  Int N = A.ncols();

  mtl::matrix<T>::type AA(M, N);

  dense1D<T> x(M), y(N);
  
  //  mtl::set_value(x, T(1));
  bool pass;

  mtl::set_value(A, T());
  mtl::set_value(AA, T());

  fill_matrix(A, A.sub(), A.super(), Shape());
  mtl::copy(A, AA);
  //  fill_matrix(AA, A.sub(), A.super(), Shape());

  Int i, j;

  //
  // rank one update test    A += x * y'
  //
  pass = true;

  for (i = 0; i < M; ++i)
    x[i] = T(i);

  for (i = 0; i < N; ++i)
    y[i] = T(i);

  // mtl version
  mtl::rank_one_update(A, x, y);

  for (i = 0; i < M; ++i)
    for (j = 0; j < N; ++j)
      AA(i,j) += x[i] * MTL_CONJ(y[j]);

  // compare A's

  pass = mtl::matrix_equal(A, AA);

  if (pass)
    cout << test.c_str() << " passed matvec rank one" << endl;
  else {
    cout << "*** matvec rank one " << test.c_str() << " failed" << endl;
#if !defined(_MSVCPP_)
    cout << "result" << endl;
    mtl::print_all_matrix(A);
    cout << "correct" << endl;
    mtl::print_all_matrix(AA);
#endif
  }
}


template <class Matrix>
void test_matvec_rankone(std::string test, Matrix& A, symmetric_tag)
{
  if (A.super() == int(A.nrows()) - 1)
    test_matvec_rankone(test, A, rectangle_tag());
  else
    cout << test.c_str() << " skipping rank one (banded symm)" << endl;
}


template <class Matrix>
void test_matvec_rankone(std::string test, Matrix& A) {
  typedef typename mtl::matrix_traits<Matrix>::shape Shape;
  test_matvec_rankone(test, A, Shape());
}


template <class Matrix>
void test_ranktwo(std::string test, Matrix& A, rectangle_tag)
{
  if (A.nrows() != A.ncols()) {
    cout << test.c_str() << " skipping rank two (A must be N x N)" << endl;
    return;
  }

  typedef typename mtl::matrix_traits<Matrix>::value_type T;
  typedef typename mtl::matrix_traits<Matrix>::size_type Int;

  Int M = A.nrows();
  Int N = A.ncols();

  mtl::matrix<T>::type AA(M, N); // column oriented

  dense1D<T> x(M), y(M);

  bool pass;

  mtl::set_value(A, T());
  mtl::set_value(AA, T());

  // start with symmetric matrices
  Int i, j;
  for (i = 0; i < A.nrows(); ++++i) {
    Int first = MTL_MAX(0, int(i) - int(A.sub()));
    Int last = MTL_MIN(int(A.ncols()), int(i) + int(A.super()) + 1);
    for (j = 0; j < A.ncols(); ++++j)
      if (j >= first && j < last) {
        A(i,j) = T(i + j);
        AA(i, j) = T(i + j);
      }
  }

  //
  // rank two update test    A += x * y' + x' * y
  //
  // doing symmetric version where y == x
  // perhaps need to do two versions of the test,
  // one of them not symmetric
  //
  pass = true;

  for (i = 0; i < M; ++i)
    x[i] = T(i);

  for (i = 0; i < N; ++i)
    y[i] = T(i);

  mtl::rank_two_update(A, x, y);

  // a correct version
  for (i = 0; i < M; ++i) {
    for (j = 0; j < i; ++j) {
      T tmp = x[i] * MTL_CONJ(y[j]) + MTL_CONJ(x[j]) * y[i];
      AA(i, j) += tmp;
      AA(j, i) += tmp;