Home | History | Annotate | Download | only in test 
      1 // This file is part of Eigen, a lightweight C++ template library
      2 // for linear algebra.
      3 //
      4 // Copyright (C) 2008-2011 Gael Guennebaud <gael.guennebaud (at) inria.fr>
      5 //
      6 // This Source Code Form is subject to the terms of the Mozilla
      7 // Public License v. 2.0. If a copy of the MPL was not distributed
      8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
      9 
     10 static long int nb_temporaries;
     11 
     12 inline void on_temporary_creation() {
     13   // here's a great place to set a breakpoint when debugging failures in this test!
     14   nb_temporaries++;
     15 }
     16 
     17 #define EIGEN_SPARSE_CREATE_TEMPORARY_PLUGIN { on_temporary_creation(); }
     18 
     19 #include "sparse.h"
     20 
     21 #define VERIFY_EVALUATION_COUNT(XPR,N) {\
     22     nb_temporaries = 0; \
     23     CALL_SUBTEST( XPR ); \
     24     if(nb_temporaries!=N) std::cerr << "nb_temporaries == " << nb_temporaries << "\n"; \
     25     VERIFY( (#XPR) && nb_temporaries==N ); \
     26   }
     27 
     28 
     29 
     30 template<typename SparseMatrixType> void sparse_product()
     31 {
     32   typedef typename SparseMatrixType::StorageIndex StorageIndex;
     33   Index n = 100;
     34   const Index rows  = internal::random<Index>(1,n);
     35   const Index cols  = internal::random<Index>(1,n);
     36   const Index depth = internal::random<Index>(1,n);
     37   typedef typename SparseMatrixType::Scalar Scalar;
     38   enum { Flags = SparseMatrixType::Flags };
     39 
     40   double density = (std::max)(8./(rows*cols), 0.2);
     41   typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
     42   typedef Matrix<Scalar,Dynamic,1> DenseVector;
     43   typedef Matrix<Scalar,1,Dynamic> RowDenseVector;
     44   typedef SparseVector<Scalar,0,StorageIndex> ColSpVector;
     45   typedef SparseVector<Scalar,RowMajor,StorageIndex> RowSpVector;
     46 
     47   Scalar s1 = internal::random<Scalar>();
     48   Scalar s2 = internal::random<Scalar>();
     49 
     50   // test matrix-matrix product
     51   {
     52     DenseMatrix refMat2  = DenseMatrix::Zero(rows, depth);
     53     DenseMatrix refMat2t = DenseMatrix::Zero(depth, rows);
     54     DenseMatrix refMat3  = DenseMatrix::Zero(depth, cols);
     55     DenseMatrix refMat3t = DenseMatrix::Zero(cols, depth);
     56     DenseMatrix refMat4  = DenseMatrix::Zero(rows, cols);
     57     DenseMatrix refMat4t = DenseMatrix::Zero(cols, rows);
     58     DenseMatrix refMat5  = DenseMatrix::Random(depth, cols);
     59     DenseMatrix refMat6  = DenseMatrix::Random(rows, rows);
     60     DenseMatrix dm4 = DenseMatrix::Zero(rows, rows);
     61 //     DenseVector dv1 = DenseVector::Random(rows);
     62     SparseMatrixType m2 (rows, depth);
     63     SparseMatrixType m2t(depth, rows);
     64     SparseMatrixType m3 (depth, cols);
     65     SparseMatrixType m3t(cols, depth);
     66     SparseMatrixType m4 (rows, cols);
     67     SparseMatrixType m4t(cols, rows);
     68     SparseMatrixType m6(rows, rows);
     69     initSparse(density, refMat2,  m2);
     70     initSparse(density, refMat2t, m2t);
     71     initSparse(density, refMat3,  m3);
     72     initSparse(density, refMat3t, m3t);
     73     initSparse(density, refMat4,  m4);
     74     initSparse(density, refMat4t, m4t);
     75     initSparse(density, refMat6, m6);
     76 
     77 //     int c = internal::random<int>(0,depth-1);
     78 
     79     // sparse * sparse
     80     VERIFY_IS_APPROX(m4=m2*m3, refMat4=refMat2*refMat3);
     81     VERIFY_IS_APPROX(m4=m2t.transpose()*m3, refMat4=refMat2t.transpose()*refMat3);
     82     VERIFY_IS_APPROX(m4=m2t.transpose()*m3t.transpose(), refMat4=refMat2t.transpose()*refMat3t.transpose());
     83     VERIFY_IS_APPROX(m4=m2*m3t.transpose(), refMat4=refMat2*refMat3t.transpose());
     84 
     85     VERIFY_IS_APPROX(m4 = m2*m3/s1, refMat4 = refMat2*refMat3/s1);
     86     VERIFY_IS_APPROX(m4 = m2*m3*s1, refMat4 = refMat2*refMat3*s1);
     87     VERIFY_IS_APPROX(m4 = s2*m2*m3*s1, refMat4 = s2*refMat2*refMat3*s1);
     88     VERIFY_IS_APPROX(m4 = (m2+m2)*m3, refMat4 = (refMat2+refMat2)*refMat3);
     89     VERIFY_IS_APPROX(m4 = m2*m3.leftCols(cols/2), refMat4 = refMat2*refMat3.leftCols(cols/2));
     90     VERIFY_IS_APPROX(m4 = m2*(m3+m3).leftCols(cols/2), refMat4 = refMat2*(refMat3+refMat3).leftCols(cols/2));
     91 
     92     VERIFY_IS_APPROX(m4=(m2*m3).pruned(0), refMat4=refMat2*refMat3);
     93     VERIFY_IS_APPROX(m4=(m2t.transpose()*m3).pruned(0), refMat4=refMat2t.transpose()*refMat3);
     94     VERIFY_IS_APPROX(m4=(m2t.transpose()*m3t.transpose()).pruned(0), refMat4=refMat2t.transpose()*refMat3t.transpose());
     95     VERIFY_IS_APPROX(m4=(m2*m3t.transpose()).pruned(0), refMat4=refMat2*refMat3t.transpose());
     96 
     97     // make sure the right product implementation is called:
     98     if((!SparseMatrixType::IsRowMajor) && m2.rows()<=m3.cols())
     99     {
    100       VERIFY_EVALUATION_COUNT(m4 = m2*m3, 3); // 1 temp for the result + 2 for transposing and get a sorted result.
    101       VERIFY_EVALUATION_COUNT(m4 = (m2*m3).pruned(0), 1);
    102       VERIFY_EVALUATION_COUNT(m4 = (m2*m3).eval().pruned(0), 4);
    103     }
    104 
    105     // and that pruning is effective:
    106     {
    107       DenseMatrix Ad(2,2);
    108       Ad << -1, 1, 1, 1;
    109       SparseMatrixType As(Ad.sparseView()), B(2,2);
    110       VERIFY_IS_EQUAL( (As*As.transpose()).eval().nonZeros(), 4);
    111       VERIFY_IS_EQUAL( (Ad*Ad.transpose()).eval().sparseView().eval().nonZeros(), 2);
    112       VERIFY_IS_EQUAL( (As*As.transpose()).pruned(1e-6).eval().nonZeros(), 2);
    113     }
    114 
    115     // dense ?= sparse * sparse
    116     VERIFY_IS_APPROX(dm4 =m2*m3, refMat4 =refMat2*refMat3);
    117     VERIFY_IS_APPROX(dm4+=m2*m3, refMat4+=refMat2*refMat3);
    118     VERIFY_IS_APPROX(dm4-=m2*m3, refMat4-=refMat2*refMat3);
    119     VERIFY_IS_APPROX(dm4 =m2t.transpose()*m3, refMat4 =refMat2t.transpose()*refMat3);
    120     VERIFY_IS_APPROX(dm4+=m2t.transpose()*m3, refMat4+=refMat2t.transpose()*refMat3);
    121     VERIFY_IS_APPROX(dm4-=m2t.transpose()*m3, refMat4-=refMat2t.transpose()*refMat3);
    122     VERIFY_IS_APPROX(dm4 =m2t.transpose()*m3t.transpose(), refMat4 =refMat2t.transpose()*refMat3t.transpose());
    123     VERIFY_IS_APPROX(dm4+=m2t.transpose()*m3t.transpose(), refMat4+=refMat2t.transpose()*refMat3t.transpose());
    124     VERIFY_IS_APPROX(dm4-=m2t.transpose()*m3t.transpose(), refMat4-=refMat2t.transpose()*refMat3t.transpose());
    125     VERIFY_IS_APPROX(dm4 =m2*m3t.transpose(), refMat4 =refMat2*refMat3t.transpose());
    126     VERIFY_IS_APPROX(dm4+=m2*m3t.transpose(), refMat4+=refMat2*refMat3t.transpose());
    127     VERIFY_IS_APPROX(dm4-=m2*m3t.transpose(), refMat4-=refMat2*refMat3t.transpose());
    128     VERIFY_IS_APPROX(dm4 = m2*m3*s1, refMat4 = refMat2*refMat3*s1);
    129 
    130     // test aliasing
    131     m4 = m2; refMat4 = refMat2;
    132     VERIFY_IS_APPROX(m4=m4*m3, refMat4=refMat4*refMat3);
    133 
    134     // sparse * dense matrix
    135     VERIFY_IS_APPROX(dm4=m2*refMat3, refMat4=refMat2*refMat3);
    136     VERIFY_IS_APPROX(dm4=m2*refMat3t.transpose(), refMat4=refMat2*refMat3t.transpose());
    137     VERIFY_IS_APPROX(dm4=m2t.transpose()*refMat3, refMat4=refMat2t.transpose()*refMat3);
    138     VERIFY_IS_APPROX(dm4=m2t.transpose()*refMat3t.transpose(), refMat4=refMat2t.transpose()*refMat3t.transpose());
    139 
    140     VERIFY_IS_APPROX(dm4=m2*refMat3, refMat4=refMat2*refMat3);
    141     VERIFY_IS_APPROX(dm4=dm4+m2*refMat3, refMat4=refMat4+refMat2*refMat3);
    142     VERIFY_IS_APPROX(dm4+=m2*refMat3, refMat4+=refMat2*refMat3);
    143     VERIFY_IS_APPROX(dm4-=m2*refMat3, refMat4-=refMat2*refMat3);
    144     VERIFY_IS_APPROX(dm4.noalias()+=m2*refMat3, refMat4+=refMat2*refMat3);
    145     VERIFY_IS_APPROX(dm4.noalias()-=m2*refMat3, refMat4-=refMat2*refMat3);
    146     VERIFY_IS_APPROX(dm4=m2*(refMat3+refMat3), refMat4=refMat2*(refMat3+refMat3));
    147     VERIFY_IS_APPROX(dm4=m2t.transpose()*(refMat3+refMat5)*0.5, refMat4=refMat2t.transpose()*(refMat3+refMat5)*0.5);
    148 
    149     // sparse * dense vector
    150     VERIFY_IS_APPROX(dm4.col(0)=m2*refMat3.col(0), refMat4.col(0)=refMat2*refMat3.col(0));
    151     VERIFY_IS_APPROX(dm4.col(0)=m2*refMat3t.transpose().col(0), refMat4.col(0)=refMat2*refMat3t.transpose().col(0));
    152     VERIFY_IS_APPROX(dm4.col(0)=m2t.transpose()*refMat3.col(0), refMat4.col(0)=refMat2t.transpose()*refMat3.col(0));
    153     VERIFY_IS_APPROX(dm4.col(0)=m2t.transpose()*refMat3t.transpose().col(0), refMat4.col(0)=refMat2t.transpose()*refMat3t.transpose().col(0));
    154 
    155     // dense * sparse
    156     VERIFY_IS_APPROX(dm4=refMat2*m3, refMat4=refMat2*refMat3);
    157     VERIFY_IS_APPROX(dm4=dm4+refMat2*m3, refMat4=refMat4+refMat2*refMat3);
    158     VERIFY_IS_APPROX(dm4+=refMat2*m3, refMat4+=refMat2*refMat3);
    159     VERIFY_IS_APPROX(dm4-=refMat2*m3, refMat4-=refMat2*refMat3);
    160     VERIFY_IS_APPROX(dm4.noalias()+=refMat2*m3, refMat4+=refMat2*refMat3);
    161     VERIFY_IS_APPROX(dm4.noalias()-=refMat2*m3, refMat4-=refMat2*refMat3);
    162     VERIFY_IS_APPROX(dm4=refMat2*m3t.transpose(), refMat4=refMat2*refMat3t.transpose());
    163     VERIFY_IS_APPROX(dm4=refMat2t.transpose()*m3, refMat4=refMat2t.transpose()*refMat3);
    164     VERIFY_IS_APPROX(dm4=refMat2t.transpose()*m3t.transpose(), refMat4=refMat2t.transpose()*refMat3t.transpose());
    165 
    166     // sparse * dense and dense * sparse outer product
    167     {
    168       Index c  = internal::random<Index>(0,depth-1);
    169       Index r  = internal::random<Index>(0,rows-1);
    170       Index c1 = internal::random<Index>(0,cols-1);
    171       Index r1 = internal::random<Index>(0,depth-1);
    172       DenseMatrix dm5  = DenseMatrix::Random(depth, cols);
    173 
    174       VERIFY_IS_APPROX( m4=m2.col(c)*dm5.col(c1).transpose(), refMat4=refMat2.col(c)*dm5.col(c1).transpose());
    175       VERIFY_IS_EQUAL(m4.nonZeros(), (refMat4.array()!=0).count());
    176       VERIFY_IS_APPROX( m4=m2.middleCols(c,1)*dm5.col(c1).transpose(), refMat4=refMat2.col(c)*dm5.col(c1).transpose());
    177       VERIFY_IS_EQUAL(m4.nonZeros(), (refMat4.array()!=0).count());
    178       VERIFY_IS_APPROX(dm4=m2.col(c)*dm5.col(c1).transpose(), refMat4=refMat2.col(c)*dm5.col(c1).transpose());
    179 
    180       VERIFY_IS_APPROX(m4=dm5.col(c1)*m2.col(c).transpose(), refMat4=dm5.col(c1)*refMat2.col(c).transpose());
    181       VERIFY_IS_EQUAL(m4.nonZeros(), (refMat4.array()!=0).count());
    182       VERIFY_IS_APPROX(m4=dm5.col(c1)*m2.middleCols(c,1).transpose(), refMat4=dm5.col(c1)*refMat2.col(c).transpose());
    183       VERIFY_IS_EQUAL(m4.nonZeros(), (refMat4.array()!=0).count());
    184       VERIFY_IS_APPROX(dm4=dm5.col(c1)*m2.col(c).transpose(), refMat4=dm5.col(c1)*refMat2.col(c).transpose());
    185 
    186       VERIFY_IS_APPROX( m4=dm5.row(r1).transpose()*m2.col(c).transpose(), refMat4=dm5.row(r1).transpose()*refMat2.col(c).transpose());
    187       VERIFY_IS_EQUAL(m4.nonZeros(), (refMat4.array()!=0).count());
    188       VERIFY_IS_APPROX(dm4=dm5.row(r1).transpose()*m2.col(c).transpose(), refMat4=dm5.row(r1).transpose()*refMat2.col(c).transpose());
    189 
    190       VERIFY_IS_APPROX( m4=m2.row(r).transpose()*dm5.col(c1).transpose(), refMat4=refMat2.row(r).transpose()*dm5.col(c1).transpose());
    191       VERIFY_IS_EQUAL(m4.nonZeros(), (refMat4.array()!=0).count());
    192       VERIFY_IS_APPROX( m4=m2.middleRows(r,1).transpose()*dm5.col(c1).transpose(), refMat4=refMat2.row(r).transpose()*dm5.col(c1).transpose());
    193       VERIFY_IS_EQUAL(m4.nonZeros(), (refMat4.array()!=0).count());
    194       VERIFY_IS_APPROX(dm4=m2.row(r).transpose()*dm5.col(c1).transpose(), refMat4=refMat2.row(r).transpose()*dm5.col(c1).transpose());
    195 
    196       VERIFY_IS_APPROX( m4=dm5.col(c1)*m2.row(r), refMat4=dm5.col(c1)*refMat2.row(r));
    197       VERIFY_IS_EQUAL(m4.nonZeros(), (refMat4.array()!=0).count());
    198       VERIFY_IS_APPROX( m4=dm5.col(c1)*m2.middleRows(r,1), refMat4=dm5.col(c1)*refMat2.row(r));
    199       VERIFY_IS_EQUAL(m4.nonZeros(), (refMat4.array()!=0).count());
    200       VERIFY_IS_APPROX(dm4=dm5.col(c1)*m2.row(r), refMat4=dm5.col(c1)*refMat2.row(r));
    201 
    202       VERIFY_IS_APPROX( m4=dm5.row(r1).transpose()*m2.row(r), refMat4=dm5.row(r1).transpose()*refMat2.row(r));
    203       VERIFY_IS_EQUAL(m4.nonZeros(), (refMat4.array()!=0).count());
    204       VERIFY_IS_APPROX(dm4=dm5.row(r1).transpose()*m2.row(r), refMat4=dm5.row(r1).transpose()*refMat2.row(r));
    205     }
    206 
    207     VERIFY_IS_APPROX(m6=m6*m6, refMat6=refMat6*refMat6);
    208 
    209     // sparse matrix * sparse vector
    210     ColSpVector cv0(cols), cv1;
    211     DenseVector dcv0(cols), dcv1;
    212     initSparse(2*density,dcv0, cv0);
    213 
    214     RowSpVector rv0(depth), rv1;
    215     RowDenseVector drv0(depth), drv1(rv1);
    216     initSparse(2*density,drv0, rv0);
    217 
    218     VERIFY_IS_APPROX(cv1=m3*cv0, dcv1=refMat3*dcv0);
    219     VERIFY_IS_APPROX(rv1=rv0*m3, drv1=drv0*refMat3);
    220     VERIFY_IS_APPROX(cv1=m3t.adjoint()*cv0, dcv1=refMat3t.adjoint()*dcv0);
    221     VERIFY_IS_APPROX(cv1=rv0*m3, dcv1=drv0*refMat3);
    222     VERIFY_IS_APPROX(rv1=m3*cv0, drv1=refMat3*dcv0);
    223   }
    224 
    225   // test matrix - diagonal product
    226   {
    227     DenseMatrix refM2 = DenseMatrix::Zero(rows, cols);
    228     DenseMatrix refM3 = DenseMatrix::Zero(rows, cols);
    229     DenseMatrix d3 = DenseMatrix::Zero(rows, cols);
    230     DiagonalMatrix<Scalar,Dynamic> d1(DenseVector::Random(cols));
    231     DiagonalMatrix<Scalar,Dynamic> d2(DenseVector::Random(rows));
    232     SparseMatrixType m2(rows, cols);
    233     SparseMatrixType m3(rows, cols);
    234     initSparse<Scalar>(density, refM2, m2);
    235     initSparse<Scalar>(density, refM3, m3);
    236     VERIFY_IS_APPROX(m3=m2*d1, refM3=refM2*d1);
    237     VERIFY_IS_APPROX(m3=m2.transpose()*d2, refM3=refM2.transpose()*d2);
    238     VERIFY_IS_APPROX(m3=d2*m2, refM3=d2*refM2);
    239     VERIFY_IS_APPROX(m3=d1*m2.transpose(), refM3=d1*refM2.transpose());
    240 
    241     // also check with a SparseWrapper:
    242     DenseVector v1 = DenseVector::Random(cols);
    243     DenseVector v2 = DenseVector::Random(rows);
    244     DenseVector v3 = DenseVector::Random(rows);
    245     VERIFY_IS_APPROX(m3=m2*v1.asDiagonal(), refM3=refM2*v1.asDiagonal());
    246     VERIFY_IS_APPROX(m3=m2.transpose()*v2.asDiagonal(), refM3=refM2.transpose()*v2.asDiagonal());
    247     VERIFY_IS_APPROX(m3=v2.asDiagonal()*m2, refM3=v2.asDiagonal()*refM2);
    248     VERIFY_IS_APPROX(m3=v1.asDiagonal()*m2.transpose(), refM3=v1.asDiagonal()*refM2.transpose());
    249 
    250     VERIFY_IS_APPROX(m3=v2.asDiagonal()*m2*v1.asDiagonal(), refM3=v2.asDiagonal()*refM2*v1.asDiagonal());
    251 
    252     VERIFY_IS_APPROX(v2=m2*v1.asDiagonal()*v1, refM2*v1.asDiagonal()*v1);
    253     VERIFY_IS_APPROX(v3=v2.asDiagonal()*m2*v1, v2.asDiagonal()*refM2*v1);
    254 
    255     // evaluate to a dense matrix to check the .row() and .col() iterator functions
    256     VERIFY_IS_APPROX(d3=m2*d1, refM3=refM2*d1);
    257     VERIFY_IS_APPROX(d3=m2.transpose()*d2, refM3=refM2.transpose()*d2);
    258     VERIFY_IS_APPROX(d3=d2*m2, refM3=d2*refM2);
    259     VERIFY_IS_APPROX(d3=d1*m2.transpose(), refM3=d1*refM2.transpose());
    260   }
    261 
    262   // test self-adjoint and triangular-view products
    263   {
    264     DenseMatrix b = DenseMatrix::Random(rows, rows);
    265     DenseMatrix x = DenseMatrix::Random(rows, rows);
    266     DenseMatrix refX = DenseMatrix::Random(rows, rows);
    267     DenseMatrix refUp = DenseMatrix::Zero(rows, rows);
    268     DenseMatrix refLo = DenseMatrix::Zero(rows, rows);
    269     DenseMatrix refS = DenseMatrix::Zero(rows, rows);
    270     DenseMatrix refA = DenseMatrix::Zero(rows, rows);
    271     SparseMatrixType mUp(rows, rows);
    272     SparseMatrixType mLo(rows, rows);
    273     SparseMatrixType mS(rows, rows);
    274     SparseMatrixType mA(rows, rows);
    275     initSparse<Scalar>(density, refA, mA);
    276     do {
    277       initSparse<Scalar>(density, refUp, mUp, ForceRealDiag|/*ForceNonZeroDiag|*/MakeUpperTriangular);
    278     } while (refUp.isZero());
    279     refLo = refUp.adjoint();
    280     mLo = mUp.adjoint();
    281     refS = refUp + refLo;
    282     refS.diagonal() *= 0.5;
    283     mS = mUp + mLo;
    284     // TODO be able to address the diagonal....
    285     for (int k=0; k<mS.outerSize(); ++k)
    286       for (typename SparseMatrixType::InnerIterator it(mS,k); it; ++it)
    287         if (it.index() == k)
    288           it.valueRef() *= Scalar(0.5);
    289 
    290     VERIFY_IS_APPROX(refS.adjoint(), refS);
    291     VERIFY_IS_APPROX(mS.adjoint(), mS);
    292     VERIFY_IS_APPROX(mS, refS);
    293     VERIFY_IS_APPROX(x=mS*b, refX=refS*b);
    294 
    295     // sparse selfadjointView with dense matrices
    296     VERIFY_IS_APPROX(x=mUp.template selfadjointView<Upper>()*b, refX=refS*b);
    297     VERIFY_IS_APPROX(x=mLo.template selfadjointView<Lower>()*b, refX=refS*b);
    298     VERIFY_IS_APPROX(x=mS.template selfadjointView<Upper|Lower>()*b, refX=refS*b);
    299 
    300     VERIFY_IS_APPROX(x=b * mUp.template selfadjointView<Upper>(),       refX=b*refS);
    301     VERIFY_IS_APPROX(x=b * mLo.template selfadjointView<Lower>(),       refX=b*refS);
    302     VERIFY_IS_APPROX(x=b * mS.template selfadjointView<Upper|Lower>(),  refX=b*refS);
    303 
    304     VERIFY_IS_APPROX(x.noalias()+=mUp.template selfadjointView<Upper>()*b, refX+=refS*b);
    305     VERIFY_IS_APPROX(x.noalias()-=mLo.template selfadjointView<Lower>()*b, refX-=refS*b);
    306     VERIFY_IS_APPROX(x.noalias()+=mS.template selfadjointView<Upper|Lower>()*b, refX+=refS*b);
    307 
    308     // sparse selfadjointView with sparse matrices
    309     SparseMatrixType mSres(rows,rows);
    310     VERIFY_IS_APPROX(mSres = mLo.template selfadjointView<Lower>()*mS,
    311                      refX = refLo.template selfadjointView<Lower>()*refS);
    312     VERIFY_IS_APPROX(mSres = mS * mLo.template selfadjointView<Lower>(),
    313                      refX = refS * refLo.template selfadjointView<Lower>());
    314 
    315     // sparse triangularView with dense matrices
    316     VERIFY_IS_APPROX(x=mA.template triangularView<Upper>()*b, refX=refA.template triangularView<Upper>()*b);
    317     VERIFY_IS_APPROX(x=mA.template triangularView<Lower>()*b, refX=refA.template triangularView<Lower>()*b);
    318     VERIFY_IS_APPROX(x=b*mA.template triangularView<Upper>(), refX=b*refA.template triangularView<Upper>());
    319     VERIFY_IS_APPROX(x=b*mA.template triangularView<Lower>(), refX=b*refA.template triangularView<Lower>());
    320 
    321     // sparse triangularView with sparse matrices
    322     VERIFY_IS_APPROX(mSres = mA.template triangularView<Lower>()*mS,   refX = refA.template triangularView<Lower>()*refS);
    323     VERIFY_IS_APPROX(mSres = mS * mA.template triangularView<Lower>(), refX = refS * refA.template triangularView<Lower>());
    324     VERIFY_IS_APPROX(mSres = mA.template triangularView<Upper>()*mS,   refX = refA.template triangularView<Upper>()*refS);
    325     VERIFY_IS_APPROX(mSres = mS * mA.template triangularView<Upper>(), refX = refS * refA.template triangularView<Upper>());
    326   }
    327 }
    328 
    329 // New test for Bug in SparseTimeDenseProduct
    330 template<typename SparseMatrixType, typename DenseMatrixType> void sparse_product_regression_test()
    331 {
    332   // This code does not compile with afflicted versions of the bug
    333   SparseMatrixType sm1(3,2);
    334   DenseMatrixType m2(2,2);
    335   sm1.setZero();
    336   m2.setZero();
    337 
    338   DenseMatrixType m3 = sm1*m2;
    339 
    340 
    341   // This code produces a segfault with afflicted versions of another SparseTimeDenseProduct
    342   // bug
    343 
    344   SparseMatrixType sm2(20000,2);
    345   sm2.setZero();
    346   DenseMatrixType m4(sm2*m2);
    347 
    348   VERIFY_IS_APPROX( m4(0,0), 0.0 );
    349 }
    350 
    351 template<typename Scalar>
    352 void bug_942()
    353 {
    354   typedef Matrix<Scalar, Dynamic, 1>     Vector;
    355   typedef SparseMatrix<Scalar, ColMajor> ColSpMat;
    356   typedef SparseMatrix<Scalar, RowMajor> RowSpMat;
    357   ColSpMat cmA(1,1);
    358   cmA.insert(0,0) = 1;
    359 
    360   RowSpMat rmA(1,1);
    361   rmA.insert(0,0) = 1;
    362 
    363   Vector d(1);
    364   d[0] = 2;
    365 
    366   double res = 2;
    367 
    368   VERIFY_IS_APPROX( ( cmA*d.asDiagonal() ).eval().coeff(0,0), res );
    369   VERIFY_IS_APPROX( ( d.asDiagonal()*rmA ).eval().coeff(0,0), res );
    370   VERIFY_IS_APPROX( ( rmA*d.asDiagonal() ).eval().coeff(0,0), res );
    371   VERIFY_IS_APPROX( ( d.asDiagonal()*cmA ).eval().coeff(0,0), res );
    372 }
    373 
    374 void test_sparse_product()
    375 {
    376   for(int i = 0; i < g_repeat; i++) {
    377     CALL_SUBTEST_1( (sparse_product<SparseMatrix<double,ColMajor> >()) );
    378     CALL_SUBTEST_1( (sparse_product<SparseMatrix<double,RowMajor> >()) );
    379     CALL_SUBTEST_1( (bug_942<double>()) );
    380     CALL_SUBTEST_2( (sparse_product<SparseMatrix<std::complex<double>, ColMajor > >()) );
    381     CALL_SUBTEST_2( (sparse_product<SparseMatrix<std::complex<double>, RowMajor > >()) );
    382     CALL_SUBTEST_3( (sparse_product<SparseMatrix<float,ColMajor,long int> >()) );
    383     CALL_SUBTEST_4( (sparse_product_regression_test<SparseMatrix<double,RowMajor>, Matrix<double, Dynamic, Dynamic, RowMajor> >()) );
    384   }
    385 }
    386