nomalloc.cpp
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 // Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
00006 //
00007 // Eigen is free software; you can redistribute it and/or
00008 // modify it under the terms of the GNU Lesser General Public
00009 // License as published by the Free Software Foundation; either
00010 // version 3 of the License, or (at your option) any later version.
00011 //
00012 // Alternatively, you can redistribute it and/or
00013 // modify it under the terms of the GNU General Public License as
00014 // published by the Free Software Foundation; either version 2 of
00015 // the License, or (at your option) any later version.
00016 //
00017 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00018 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00019 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00020 // GNU General Public License for more details.
00021 //
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License and a copy of the GNU General Public License along with
00024 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00025 
00026 // this hack is needed to make this file compiles with -pedantic (gcc)
00027 #ifdef __GNUC__
00028 #define throw(X)
00029 #endif
00030 // discard stack allocation as that too bypasses malloc
00031 #define EIGEN_STACK_ALLOCATION_LIMIT 0
00032 // any heap allocation will raise an assert
00033 #define EIGEN_NO_MALLOC
00034 
00035 #include "main.h"
00036 #include <Eigen/Cholesky>
00037 #include <Eigen/Eigenvalues>
00038 #include <Eigen/LU>
00039 #include <Eigen/QR>
00040 #include <Eigen/SVD>
00041 
00042 template<typename MatrixType> void nomalloc(const MatrixType& m)
00043 {
00044   /* this test check no dynamic memory allocation are issued with fixed-size matrices
00045   */
00046   typedef typename MatrixType::Index Index;
00047   typedef typename MatrixType::Scalar Scalar;
00048   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
00049 
00050   Index rows = m.rows();
00051   Index cols = m.cols();
00052 
00053   MatrixType m1 = MatrixType::Random(rows, cols),
00054              m2 = MatrixType::Random(rows, cols),
00055              m3(rows, cols),
00056              mzero = MatrixType::Zero(rows, cols),
00057              identity = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime>
00058                               ::Identity(rows, rows),
00059              square = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime>
00060                               ::Random(rows, rows);
00061   VectorType v1 = VectorType::Random(rows),
00062              v2 = VectorType::Random(rows),
00063              vzero = VectorType::Zero(rows);
00064 
00065   Scalar s1 = internal::random<Scalar>();
00066 
00067   Index r = internal::random<Index>(0, rows-1),
00068         c = internal::random<Index>(0, cols-1);
00069 
00070   VERIFY_IS_APPROX((m1+m2)*s1,              s1*m1+s1*m2);
00071   VERIFY_IS_APPROX((m1+m2)(r,c), (m1(r,c))+(m2(r,c)));
00072   VERIFY_IS_APPROX(m1.cwiseProduct(m1.block(0,0,rows,cols)), (m1.array()*m1.array()).matrix());
00073   VERIFY_IS_APPROX((m1*m1.transpose())*m2,  m1*(m1.transpose()*m2));
00074   
00075   m2.col(0).noalias() = m1 * m1.col(0);
00076   m2.col(0).noalias() -= m1.adjoint() * m1.col(0);
00077   m2.col(0).noalias() -= m1 * m1.row(0).adjoint();
00078   m2.col(0).noalias() -= m1.adjoint() * m1.row(0).adjoint();
00079 
00080   m2.row(0).noalias() = m1.row(0) * m1;
00081   m2.row(0).noalias() -= m1.row(0) * m1.adjoint();
00082   m2.row(0).noalias() -= m1.col(0).adjoint() * m1;
00083   m2.row(0).noalias() -= m1.col(0).adjoint() * m1.adjoint();
00084   VERIFY_IS_APPROX(m2,m2);
00085   
00086   m2.col(0).noalias() = m1.template triangularView<Upper>() * m1.col(0);
00087   m2.col(0).noalias() -= m1.adjoint().template triangularView<Upper>() * m1.col(0);
00088   m2.col(0).noalias() -= m1.template triangularView<Upper>() * m1.row(0).adjoint();
00089   m2.col(0).noalias() -= m1.adjoint().template triangularView<Upper>() * m1.row(0).adjoint();
00090 
00091   m2.row(0).noalias() = m1.row(0) * m1.template triangularView<Upper>();
00092   m2.row(0).noalias() -= m1.row(0) * m1.adjoint().template triangularView<Upper>();
00093   m2.row(0).noalias() -= m1.col(0).adjoint() * m1.template triangularView<Upper>();
00094   m2.row(0).noalias() -= m1.col(0).adjoint() * m1.adjoint().template triangularView<Upper>();
00095   VERIFY_IS_APPROX(m2,m2);
00096   
00097   m2.col(0).noalias() = m1.template selfadjointView<Upper>() * m1.col(0);
00098   m2.col(0).noalias() -= m1.adjoint().template selfadjointView<Upper>() * m1.col(0);
00099   m2.col(0).noalias() -= m1.template selfadjointView<Upper>() * m1.row(0).adjoint();
00100   m2.col(0).noalias() -= m1.adjoint().template selfadjointView<Upper>() * m1.row(0).adjoint();
00101 
00102   m2.row(0).noalias() = m1.row(0) * m1.template selfadjointView<Upper>();
00103   m2.row(0).noalias() -= m1.row(0) * m1.adjoint().template selfadjointView<Upper>();
00104   m2.row(0).noalias() -= m1.col(0).adjoint() * m1.template selfadjointView<Upper>();
00105   m2.row(0).noalias() -= m1.col(0).adjoint() * m1.adjoint().template selfadjointView<Upper>();
00106   VERIFY_IS_APPROX(m2,m2);
00107   
00108   m2.template selfadjointView<Lower>().rankUpdate(m1.col(0),-1);
00109   m2.template selfadjointView<Lower>().rankUpdate(m1.row(0),-1);
00110 
00111   // The following fancy matrix-matrix products are not safe yet regarding static allocation
00112 //   m1 += m1.template triangularView<Upper>() * m2.col(;
00113 //   m1.template selfadjointView<Lower>().rankUpdate(m2);
00114 //   m1 += m1.template triangularView<Upper>() * m2;
00115 //   m1 += m1.template selfadjointView<Lower>() * m2;
00116 //   VERIFY_IS_APPROX(m1,m1);
00117 }
00118 
00119 template<typename Scalar>
00120 void ctms_decompositions()
00121 {
00122   const int maxSize = 16;
00123   const int size    = 12;
00124 
00125   typedef Eigen::Matrix<Scalar,
00126                         Eigen::Dynamic, Eigen::Dynamic,
00127                         0,
00128                         maxSize, maxSize> Matrix;
00129 
00130   typedef Eigen::Matrix<Scalar,
00131                         Eigen::Dynamic, 1,
00132                         0,
00133                         maxSize, 1> Vector;
00134 
00135   typedef Eigen::Matrix<std::complex<Scalar>,
00136                         Eigen::Dynamic, Eigen::Dynamic,
00137                         0,
00138                         maxSize, maxSize> ComplexMatrix;
00139 
00140   const Matrix A(Matrix::Random(size, size));
00141   const ComplexMatrix complexA(ComplexMatrix::Random(size, size));
00142   const Matrix saA = A.adjoint() * A;
00143 
00144   // Cholesky module
00145   Eigen::LLT<Matrix>  LLT;  LLT.compute(A);
00146   Eigen::LDLT<Matrix> LDLT; LDLT.compute(A);
00147 
00148   // Eigenvalues module
00149   Eigen::HessenbergDecomposition<ComplexMatrix> hessDecomp;        hessDecomp.compute(complexA);
00150   Eigen::ComplexSchur<ComplexMatrix>            cSchur(size);      cSchur.compute(complexA);
00151   Eigen::ComplexEigenSolver<ComplexMatrix>      cEigSolver;        cEigSolver.compute(complexA);
00152   Eigen::EigenSolver<Matrix>                    eigSolver;         eigSolver.compute(A);
00153   Eigen::SelfAdjointEigenSolver<Matrix>         saEigSolver(size); saEigSolver.compute(saA);
00154   Eigen::Tridiagonalization<Matrix>             tridiag;           tridiag.compute(saA);
00155 
00156   // LU module
00157   Eigen::PartialPivLU<Matrix> ppLU; ppLU.compute(A);
00158   Eigen::FullPivLU<Matrix>    fpLU; fpLU.compute(A);
00159 
00160   // QR module
00161   Eigen::HouseholderQR<Matrix>        hQR;  hQR.compute(A);
00162   Eigen::ColPivHouseholderQR<Matrix>  cpQR; cpQR.compute(A);
00163   Eigen::FullPivHouseholderQR<Matrix> fpQR; fpQR.compute(A);
00164 
00165   // SVD module
00166   Eigen::JacobiSVD<Matrix> jSVD; jSVD.compute(A, ComputeFullU | ComputeFullV);
00167 }
00168 
00169 void test_nomalloc()
00170 {
00171   // check that our operator new is indeed called:
00172   VERIFY_RAISES_ASSERT(MatrixXd dummy(MatrixXd::Random(3,3)));
00173   CALL_SUBTEST_1(nomalloc(Matrix<float, 1, 1>()) );
00174   CALL_SUBTEST_2(nomalloc(Matrix4d()) );
00175   CALL_SUBTEST_3(nomalloc(Matrix<float,32,32>()) );
00176   
00177   // Check decomposition modules with dynamic matrices that have a known compile-time max size (ctms)
00178   CALL_SUBTEST_4(ctms_decompositions<float>());
00179 
00180 }


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:33:07