vlsvrs 1.0.0

VLSV reader with C and Python bindings
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// 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 "product.h"
#include <Eigen/LU>

template<typename T>
void test_aliasing()
{
  int rows = internal::random<int>(1,12);
  int cols = internal::random<int>(1,12);
  typedef Matrix<T,Dynamic,Dynamic> MatrixType;
  typedef Matrix<T,Dynamic,1> VectorType;
  VectorType x(cols); x.setRandom();
  VectorType z(x);
  VectorType y(rows); y.setZero();
  MatrixType A(rows,cols); A.setRandom();
  // CwiseBinaryOp
  VERIFY_IS_APPROX(x = y + A*x, A*z);     // OK because "y + A*x" is marked as "assume-aliasing"
  x = z;
  // CwiseUnaryOp
  VERIFY_IS_APPROX(x = T(1.)*(A*x), A*z); // OK because 1*(A*x) is replaced by (1*A*x) which is a Product<> expression
  x = z;
  // VERIFY_IS_APPROX(x = y-A*x, -A*z);   // Not OK in 3.3 because x is resized before A*x gets evaluated
  x = z;
}

template<int>
void product_large_regressions()
{
  {
    // test a specific issue in DiagonalProduct
    int N = 1000000;
    VectorXf v = VectorXf::Ones(N);
    MatrixXf m = MatrixXf::Ones(N,3);
    m = (v+v).asDiagonal() * m;
    VERIFY_IS_APPROX(m, MatrixXf::Constant(N,3,2));
  }

  {
    // test deferred resizing in Matrix::operator=
    MatrixXf a = MatrixXf::Random(10,4), b = MatrixXf::Random(4,10), c = a;
    VERIFY_IS_APPROX((a = a * b), (c * b).eval());
  }

  {
    // check the functions to setup blocking sizes compile and do not segfault
    // FIXME check they do what they are supposed to do !!
    std::ptrdiff_t l1 = internal::random<int>(10000,20000);
    std::ptrdiff_t l2 = internal::random<int>(100000,200000);
    std::ptrdiff_t l3 = internal::random<int>(1000000,2000000);
    setCpuCacheSizes(l1,l2,l3);
    VERIFY(l1==l1CacheSize());
    VERIFY(l2==l2CacheSize());
    std::ptrdiff_t k1 = internal::random<int>(10,100)*16;
    std::ptrdiff_t m1 = internal::random<int>(10,100)*16;
    std::ptrdiff_t n1 = internal::random<int>(10,100)*16;
    // only makes sure it compiles fine
    internal::computeProductBlockingSizes<float,float,std::ptrdiff_t>(k1,m1,n1,1);
  }

  {
    // test regression in row-vector by matrix (bad Map type)
    MatrixXf mat1(10,32); mat1.setRandom();
    MatrixXf mat2(32,32); mat2.setRandom();
    MatrixXf r1 = mat1.row(2)*mat2.transpose();
    VERIFY_IS_APPROX(r1, (mat1.row(2)*mat2.transpose()).eval());

    MatrixXf r2 = mat1.row(2)*mat2;
    VERIFY_IS_APPROX(r2, (mat1.row(2)*mat2).eval());
  }

  {
    Eigen::MatrixXd A(10,10), B, C;
    A.setRandom();
    C = A;
    for(int k=0; k<79; ++k)
      C = C * A;
    B.noalias() = (((A*A)*(A*A))*((A*A)*(A*A))*((A*A)*(A*A))*((A*A)*(A*A))*((A*A)*(A*A)) * ((A*A)*(A*A))*((A*A)*(A*A))*((A*A)*(A*A))*((A*A)*(A*A))*((A*A)*(A*A)))
                * (((A*A)*(A*A))*((A*A)*(A*A))*((A*A)*(A*A))*((A*A)*(A*A))*((A*A)*(A*A)) * ((A*A)*(A*A))*((A*A)*(A*A))*((A*A)*(A*A))*((A*A)*(A*A))*((A*A)*(A*A)));
    VERIFY_IS_APPROX(B,C);
  }
}

template<int>
void bug_1622() {
  typedef Matrix<double, 2, -1, 0, 2, -1> Mat2X;
  Mat2X x(2,2); x.setRandom();
  MatrixXd y(2,2); y.setRandom();
  const Mat2X K1 = x * y.inverse();
  const Matrix2d K2 = x * y.inverse();
  VERIFY_IS_APPROX(K1,K2);
}

EIGEN_DECLARE_TEST(product_large)
{
  for(int i = 0; i < g_repeat; i++) {
    CALL_SUBTEST_1( product(MatrixXf(internal::random<int>(1,EIGEN_TEST_MAX_SIZE), internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
    CALL_SUBTEST_2( product(MatrixXd(internal::random<int>(1,EIGEN_TEST_MAX_SIZE), internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
    CALL_SUBTEST_2( product(MatrixXd(internal::random<int>(1,10), internal::random<int>(1,10))) );

    CALL_SUBTEST_3( product(MatrixXi(internal::random<int>(1,EIGEN_TEST_MAX_SIZE), internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
    CALL_SUBTEST_4( product(MatrixXcf(internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2), internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2))) );
    CALL_SUBTEST_5( product(Matrix<float,Dynamic,Dynamic,RowMajor>(internal::random<int>(1,EIGEN_TEST_MAX_SIZE), internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );

    CALL_SUBTEST_1( test_aliasing<float>() );

    CALL_SUBTEST_6( bug_1622<1>() );

    CALL_SUBTEST_7( product(MatrixXcd(internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2), internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2))) );
    CALL_SUBTEST_8( product(Matrix<double,Dynamic,Dynamic,RowMajor>(internal::random<int>(1,EIGEN_TEST_MAX_SIZE), internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
    CALL_SUBTEST_9( product(Matrix<std::complex<float>,Dynamic,Dynamic,RowMajor>(internal::random<int>(1,EIGEN_TEST_MAX_SIZE), internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
    CALL_SUBTEST_10( product(Matrix<std::complex<double>,Dynamic,Dynamic,RowMajor>(internal::random<int>(1,EIGEN_TEST_MAX_SIZE), internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
  }

  CALL_SUBTEST_6( product_large_regressions<0>() );

  // Regression test for bug 714:
#if defined EIGEN_HAS_OPENMP
  omp_set_dynamic(1);
  for(int i = 0; i < g_repeat; i++) {
    CALL_SUBTEST_6( product(Matrix<float,Dynamic,Dynamic>(internal::random<int>(1,EIGEN_TEST_MAX_SIZE), internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
  }
#endif
}