ground/gcs/src/libs/eigen/test/sparse_product.cpp

   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@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.noalias()+=mUp.template selfadjointView<Upper>()*b, refX+=refS*b);
 301     VERIFY_IS_APPROX(x.noalias()-=mLo.template selfadjointView<Lower>()*b, refX-=refS*b);
 302     VERIFY_IS_APPROX(x.noalias()+=mS.template selfadjointView<Upper|Lower>()*b, refX+=refS*b);
 303
 304     // sparse selfadjointView with sparse matrices
 305     SparseMatrixType mSres(rows,rows);
 306     VERIFY_IS_APPROX(mSres = mLo.template selfadjointView<Lower>()*mS,
 307                      refX = refLo.template selfadjointView<Lower>()*refS);
 308     VERIFY_IS_APPROX(mSres = mS * mLo.template selfadjointView<Lower>(),
 309                      refX = refS * refLo.template selfadjointView<Lower>());
 310
 311     // sparse triangularView with dense matrices
 312     VERIFY_IS_APPROX(x=mA.template triangularView<Upper>()*b, refX=refA.template triangularView<Upper>()*b);
 313     VERIFY_IS_APPROX(x=mA.template triangularView<Lower>()*b, refX=refA.template triangularView<Lower>()*b);
 314     VERIFY_IS_APPROX(x=b*mA.template triangularView<Upper>(), refX=b*refA.template triangularView<Upper>());
 315     VERIFY_IS_APPROX(x=b*mA.template triangularView<Lower>(), refX=b*refA.template triangularView<Lower>());
 316
 317     // sparse triangularView with sparse matrices
 318     VERIFY_IS_APPROX(mSres = mA.template triangularView<Lower>()*mS,   refX = refA.template triangularView<Lower>()*refS);
 319     VERIFY_IS_APPROX(mSres = mS * mA.template triangularView<Lower>(), refX = refS * refA.template triangularView<Lower>());
 320     VERIFY_IS_APPROX(mSres = mA.template triangularView<Upper>()*mS,   refX = refA.template triangularView<Upper>()*refS);
 321     VERIFY_IS_APPROX(mSres = mS * mA.template triangularView<Upper>(), refX = refS * refA.template triangularView<Upper>());
 322   }
 323 }
 324
 325 // New test for Bug in SparseTimeDenseProduct
 326 template<typename SparseMatrixType, typename DenseMatrixType> void sparse_product_regression_test()
 327 {
 328   // This code does not compile with afflicted versions of the bug
 329   SparseMatrixType sm1(3,2);
 330   DenseMatrixType m2(2,2);
 331   sm1.setZero();
 332   m2.setZero();
 333
 334   DenseMatrixType m3 = sm1*m2;
 335
 336
 337   // This code produces a segfault with afflicted versions of another SparseTimeDenseProduct
 338   // bug
 339
 340   SparseMatrixType sm2(20000,2);
 341   sm2.setZero();
 342   DenseMatrixType m4(sm2*m2);
 343
 344   VERIFY_IS_APPROX( m4(0,0), 0.0 );
 345 }
 346
 347 template<typename Scalar>
 348 void bug_942()
 349 {
 350   typedef Matrix<Scalar, Dynamic, 1>     Vector;
 351   typedef SparseMatrix<Scalar, ColMajor> ColSpMat;
 352   typedef SparseMatrix<Scalar, RowMajor> RowSpMat;
 353   ColSpMat cmA(1,1);
 354   cmA.insert(0,0) = 1;
 355
 356   RowSpMat rmA(1,1);
 357   rmA.insert(0,0) = 1;
 358
 359   Vector d(1);
 360   d[0] = 2;
 361
 362   double res = 2;
 363
 364   VERIFY_IS_APPROX( ( cmA*d.asDiagonal() ).eval().coeff(0,0), res );
 365   VERIFY_IS_APPROX( ( d.asDiagonal()*rmA ).eval().coeff(0,0), res );
 366   VERIFY_IS_APPROX( ( rmA*d.asDiagonal() ).eval().coeff(0,0), res );
 367   VERIFY_IS_APPROX( ( d.asDiagonal()*cmA ).eval().coeff(0,0), res );
 368 }
 369
 370 void test_sparse_product()
 371 {
 372   for(int i = 0; i < g_repeat; i++) {
 373     CALL_SUBTEST_1( (sparse_product<SparseMatrix<double,ColMajor> >()) );
 374     CALL_SUBTEST_1( (sparse_product<SparseMatrix<double,RowMajor> >()) );
 375     CALL_SUBTEST_1( (bug_942<double>()) );
 376     CALL_SUBTEST_2( (sparse_product<SparseMatrix<std::complex<double>, ColMajor > >()) );
 377     CALL_SUBTEST_2( (sparse_product<SparseMatrix<std::complex<double>, RowMajor > >()) );
 378     CALL_SUBTEST_3( (sparse_product<SparseMatrix<float,ColMajor,long int> >()) );
 379     CALL_SUBTEST_4( (sparse_product_regression_test<SparseMatrix<double,RowMajor>, Matrix<double, Dynamic, Dynamic, RowMajor> >()) );
 380   }
 381 }