C++ MatrixType类代码示例

template<typename MatrixType> void vectorwiseop_matrix(const MatrixType& m)
{
  typedef typename MatrixType::Index Index;
  typedef typename MatrixType::Scalar Scalar;
  typedef typename NumTraits<Scalar>::Real RealScalar;
  typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> ColVectorType;
  typedef Matrix<Scalar, 1, MatrixType::ColsAtCompileTime> RowVectorType;
  typedef Matrix<RealScalar, MatrixType::RowsAtCompileTime, 1> RealColVectorType;
  typedef Matrix<RealScalar, 1, MatrixType::ColsAtCompileTime> RealRowVectorType;

  Index rows = m.rows();
  Index cols = m.cols();
  Index r = internal::random<Index>(0, rows-1),
        c = internal::random<Index>(0, cols-1);

  MatrixType m1 = MatrixType::Random(rows, cols),
            m2(rows, cols),
            m3(rows, cols);

  ColVectorType colvec = ColVectorType::Random(rows);
  RowVectorType rowvec = RowVectorType::Random(cols);
  RealColVectorType rcres;
  RealRowVectorType rrres;

  // test addition

  m2 = m1;
  m2.colwise() += colvec;
  VERIFY_IS_APPROX(m2, m1.colwise() + colvec);
  VERIFY_IS_APPROX(m2.col(c), m1.col(c) + colvec);

  VERIFY_RAISES_ASSERT(m2.colwise() += colvec.transpose());
  VERIFY_RAISES_ASSERT(m1.colwise() + colvec.transpose());

  m2 = m1;
  m2.rowwise() += rowvec;
  VERIFY_IS_APPROX(m2, m1.rowwise() + rowvec);
  VERIFY_IS_APPROX(m2.row(r), m1.row(r) + rowvec);

  VERIFY_RAISES_ASSERT(m2.rowwise() += rowvec.transpose());
  VERIFY_RAISES_ASSERT(m1.rowwise() + rowvec.transpose());

  // test substraction

  m2 = m1;
  m2.colwise() -= colvec;
  VERIFY_IS_APPROX(m2, m1.colwise() - colvec);
  VERIFY_IS_APPROX(m2.col(c), m1.col(c) - colvec);

  VERIFY_RAISES_ASSERT(m2.colwise() -= colvec.transpose());
  VERIFY_RAISES_ASSERT(m1.colwise() - colvec.transpose());

  m2 = m1;
  m2.rowwise() -= rowvec;
  VERIFY_IS_APPROX(m2, m1.rowwise() - rowvec);
  VERIFY_IS_APPROX(m2.row(r), m1.row(r) - rowvec);

  VERIFY_RAISES_ASSERT(m2.rowwise() -= rowvec.transpose());
  VERIFY_RAISES_ASSERT(m1.rowwise() - rowvec.transpose());

  // test norm
  rrres = m1.colwise().norm();
  VERIFY_IS_APPROX(rrres(c), m1.col(c).norm());
  rcres = m1.rowwise().norm();
  VERIFY_IS_APPROX(rcres(r), m1.row(r).norm());

  // test normalized
  m2 = m1.colwise().normalized();
  VERIFY_IS_APPROX(m2.col(c), m1.col(c).normalized());
  m2 = m1.rowwise().normalized();
  VERIFY_IS_APPROX(m2.row(r), m1.row(r).normalized());

  // test normalize
  m2 = m1;
  m2.colwise().normalize();
  VERIFY_IS_APPROX(m2.col(c), m1.col(c).normalized());
  m2 = m1;
  m2.rowwise().normalize();
  VERIFY_IS_APPROX(m2.row(r), m1.row(r).normalized());
}

 // returns the # of cols of the matrix *seen* as fortran matrix
 template <typename MatrixType> int get_n_cols (MatrixType const & A) { 
  return (A.memory_layout_is_fortran() ? second_dim(A) : first_dim(A));
 }

__attribute__ ((noinline)) void benchLLT(const MatrixType& m)
{
  int rows = m.rows();
  int cols = m.cols();

  double cost = 0;
  for (int j=0; j<rows; ++j)
  {
    int r = std::max(rows - j -1,0);
    cost += 2*(r*j+r+j);
  }

  int repeats = (REPEAT*1000)/(rows*rows);

  typedef typename MatrixType::Scalar Scalar;
  typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> SquareMatrixType;

  MatrixType a = MatrixType::Random(rows,cols);
  SquareMatrixType covMat =  a * a.adjoint();

  BenchTimer timerNoSqrt, timerSqrt;

  Scalar acc = 0;
  int r = internal::random<int>(0,covMat.rows()-1);
  int c = internal::random<int>(0,covMat.cols()-1);
  for (int t=0; t<TRIES; ++t)
  {
    timerNoSqrt.start();
    for (int k=0; k<repeats; ++k)
    {
      LDLT<SquareMatrixType> cholnosqrt(covMat);
      acc += cholnosqrt.matrixL().coeff(r,c);
    }
    timerNoSqrt.stop();
  }

  for (int t=0; t<TRIES; ++t)
  {
    timerSqrt.start();
    for (int k=0; k<repeats; ++k)
    {
      LLT<SquareMatrixType> chol(covMat);
      acc += chol.matrixL().coeff(r,c);
    }
    timerSqrt.stop();
  }

  if (MatrixType::RowsAtCompileTime==Dynamic)
    std::cout << "dyn   ";
  else
    std::cout << "fixed ";
  std::cout << covMat.rows() << " \t"
            << (timerNoSqrt.best()) / repeats << "s "
            << "(" << 1e-9 * cost*repeats/timerNoSqrt.best() << " GFLOPS)\t"
            << (timerSqrt.best()) / repeats << "s "
            << "(" << 1e-9 * cost*repeats/timerSqrt.best() << " GFLOPS)\n";


  #ifdef BENCH_GSL
  if (MatrixType::RowsAtCompileTime==Dynamic)
  {
    timerSqrt.reset();

    gsl_matrix* gslCovMat = gsl_matrix_alloc(covMat.rows(),covMat.cols());
    gsl_matrix* gslCopy = gsl_matrix_alloc(covMat.rows(),covMat.cols());

    eiToGsl(covMat, &gslCovMat);
    for (int t=0; t<TRIES; ++t)
    {
      timerSqrt.start();
      for (int k=0; k<repeats; ++k)
      {
        gsl_matrix_memcpy(gslCopy,gslCovMat);
        gsl_linalg_cholesky_decomp(gslCopy);
        acc += gsl_matrix_get(gslCopy,r,c);
      }
      timerSqrt.stop();
    }

    std::cout << " | \t"
              << timerSqrt.value() * REPEAT / repeats << "s";

    gsl_matrix_free(gslCovMat);
  }
  #endif
  std::cout << "\n";
  // make sure the compiler does not optimize too much
  if (acc==123)
    std::cout << acc;
}

 /** @brief Constructor for the preconditioner
 *
 * @param mat   The system matrix
 * @param tag   A row scaling tag holding the desired norm.
 */
 row_scaling(MatrixType const & mat, row_scaling_tag const & tag) : diag_M(viennacl::traits::size1(mat))
 {
   assert(mat.size1() == mat.size2() && bool("Size mismatch"));
   init(mat, tag);
 }

template<typename PointT> template <typename MatrixType> void
pcl::search::OrganizedNeighbor<PointT>::makeSymmetric (MatrixType& matrix, bool use_upper_triangular) const
{
  if (use_upper_triangular)
  {
    matrix.coeffRef (4) = matrix.coeff (1);
    matrix.coeffRef (8) = matrix.coeff (2);
    matrix.coeffRef (9) = matrix.coeff (6);
    matrix.coeffRef (12) = matrix.coeff (3);
    matrix.coeffRef (13) = matrix.coeff (7);
    matrix.coeffRef (14) = matrix.coeff (11);
  }
  else
  {
    matrix.coeffRef (1) = matrix.coeff (4);
    matrix.coeffRef (2) = matrix.coeff (8);
    matrix.coeffRef (6) = matrix.coeff (9);
    matrix.coeffRef (3) = matrix.coeff (12);
    matrix.coeffRef (7) = matrix.coeff (13);
    matrix.coeffRef (11) = matrix.coeff (14);
  }
}

template<typename MatrixType> void product_notemporary(const MatrixType& m)
{
    /* This test checks the number of temporaries created
     * during the evaluation of a complex expression */
    typedef typename MatrixType::Index Index;
    typedef typename MatrixType::Scalar Scalar;
    typedef typename MatrixType::RealScalar RealScalar;
    typedef Matrix<Scalar, 1, Dynamic> RowVectorType;
    typedef Matrix<Scalar, Dynamic, 1> ColVectorType;
    typedef Matrix<Scalar, Dynamic, Dynamic, ColMajor> ColMajorMatrixType;
    typedef Matrix<Scalar, Dynamic, Dynamic, RowMajor> RowMajorMatrixType;

    Index rows = m.rows();
    Index cols = m.cols();

    ColMajorMatrixType m1 = MatrixType::Random(rows, cols),
                       m2 = MatrixType::Random(rows, cols),
                       m3(rows, cols);
    RowVectorType rv1 = RowVectorType::Random(rows), rvres(rows);
    ColVectorType cv1 = ColVectorType::Random(cols), cvres(cols);
    RowMajorMatrixType rm3(rows, cols);

    Scalar s1 = internal::random<Scalar>(),
           s2 = internal::random<Scalar>(),
           s3 = internal::random<Scalar>();

    Index c0 = internal::random<Index>(4,cols-8),
          c1 = internal::random<Index>(8,cols-c0),
          r0 = internal::random<Index>(4,cols-8),
          r1 = internal::random<Index>(8,rows-r0);

    VERIFY_EVALUATION_COUNT( m3 = (m1 * m2.adjoint()), 1);
    VERIFY_EVALUATION_COUNT( m3 = (m1 * m2.adjoint()).transpose(), 1);
    VERIFY_EVALUATION_COUNT( m3.noalias() = m1 * m2.adjoint(), 0);

    VERIFY_EVALUATION_COUNT( m3 = s1 * (m1 * m2.transpose()), 1);
//   VERIFY_EVALUATION_COUNT( m3 = m3 + s1 * (m1 * m2.transpose()), 1);
    VERIFY_EVALUATION_COUNT( m3.noalias() = s1 * (m1 * m2.transpose()), 0);

    VERIFY_EVALUATION_COUNT( m3 = m3 + (m1 * m2.adjoint()), 1);

    VERIFY_EVALUATION_COUNT( m3 = m3 + (m1 * m2.adjoint()).transpose(), 1);
    VERIFY_EVALUATION_COUNT( m3.noalias() = m3 + m1 * m2.transpose(), 0);
    VERIFY_EVALUATION_COUNT( m3.noalias() += m3 + m1 * m2.transpose(), 0);
    VERIFY_EVALUATION_COUNT( m3.noalias() -= m3 + m1 * m2.transpose(), 0);

    VERIFY_EVALUATION_COUNT( m3.noalias() = s1 * m1 * s2 * m2.adjoint(), 0);
    VERIFY_EVALUATION_COUNT( m3.noalias() = s1 * m1 * s2 * (m1*s3+m2*s2).adjoint(), 1);
    VERIFY_EVALUATION_COUNT( m3.noalias() = (s1 * m1).adjoint() * s2 * m2, 0);
    VERIFY_EVALUATION_COUNT( m3.noalias() += s1 * (-m1*s3).adjoint() * (s2 * m2 * s3), 0);
    VERIFY_EVALUATION_COUNT( m3.noalias() -= s1 * (m1.transpose() * m2), 0);

    VERIFY_EVALUATION_COUNT(( m3.block(r0,r0,r1,r1).noalias() += -m1.block(r0,c0,r1,c1) * (s2*m2.block(r0,c0,r1,c1)).adjoint() ), 0);
    VERIFY_EVALUATION_COUNT(( m3.block(r0,r0,r1,r1).noalias() -= s1 * m1.block(r0,c0,r1,c1) * m2.block(c0,r0,c1,r1) ), 0);

    // NOTE this is because the Block expression is not handled yet by our expression analyser
    VERIFY_EVALUATION_COUNT(( m3.block(r0,r0,r1,r1).noalias() = s1 * m1.block(r0,c0,r1,c1) * (s1*m2).block(c0,r0,c1,r1) ), 1);

    VERIFY_EVALUATION_COUNT( m3.noalias() -= (s1 * m1).template triangularView<Lower>() * m2, 0);
    VERIFY_EVALUATION_COUNT( rm3.noalias() = (s1 * m1.adjoint()).template triangularView<Upper>() * (m2+m2), 1);
    VERIFY_EVALUATION_COUNT( rm3.noalias() = (s1 * m1.adjoint()).template triangularView<UnitUpper>() * m2.adjoint(), 0);

    VERIFY_EVALUATION_COUNT( m3.template triangularView<Upper>() = (m1 * m2.adjoint()), 0);
    VERIFY_EVALUATION_COUNT( m3.template triangularView<Upper>() -= (m1 * m2.adjoint()), 0);

    // NOTE this is because the blas_traits require innerstride==1 to avoid a temporary, but that doesn't seem to be actually needed for the triangular products
    VERIFY_EVALUATION_COUNT( rm3.col(c0).noalias() = (s1 * m1.adjoint()).template triangularView<UnitUpper>() * (s2*m2.row(c0)).adjoint(), 1);

    VERIFY_EVALUATION_COUNT( m1.template triangularView<Lower>().solveInPlace(m3), 0);
    VERIFY_EVALUATION_COUNT( m1.adjoint().template triangularView<Lower>().solveInPlace(m3.transpose()), 0);

    VERIFY_EVALUATION_COUNT( m3.noalias() -= (s1 * m1).adjoint().template selfadjointView<Lower>() * (-m2*s3).adjoint(), 0);
    VERIFY_EVALUATION_COUNT( m3.noalias() = s2 * m2.adjoint() * (s1 * m1.adjoint()).template selfadjointView<Upper>(), 0);
    VERIFY_EVALUATION_COUNT( rm3.noalias() = (s1 * m1.adjoint()).template selfadjointView<Lower>() * m2.adjoint(), 0);

    // NOTE this is because the blas_traits require innerstride==1 to avoid a temporary, but that doesn't seem to be actually needed for the triangular products
    VERIFY_EVALUATION_COUNT( m3.col(c0).noalias() = (s1 * m1).adjoint().template selfadjointView<Lower>() * (-m2.row(c0)*s3).adjoint(), 1);
    VERIFY_EVALUATION_COUNT( m3.col(c0).noalias() -= (s1 * m1).adjoint().template selfadjointView<Upper>() * (-m2.row(c0)*s3).adjoint(), 1);

    VERIFY_EVALUATION_COUNT( m3.block(r0,c0,r1,c1).noalias() += m1.block(r0,r0,r1,r1).template selfadjointView<Upper>() * (s1*m2.block(r0,c0,r1,c1)), 0);
    VERIFY_EVALUATION_COUNT( m3.block(r0,c0,r1,c1).noalias() = m1.block(r0,r0,r1,r1).template selfadjointView<Upper>() * m2.block(r0,c0,r1,c1), 0);

    VERIFY_EVALUATION_COUNT( m3.template selfadjointView<Lower>().rankUpdate(m2.adjoint()), 0);

    // Here we will get 1 temporary for each resize operation of the lhs operator; resize(r1,c1) would lead to zero temporaries
    m3.resize(1,1);
    VERIFY_EVALUATION_COUNT( m3.noalias() = m1.block(r0,r0,r1,r1).template selfadjointView<Lower>() * m2.block(r0,c0,r1,c1), 1);
    m3.resize(1,1);
    VERIFY_EVALUATION_COUNT( m3.noalias() = m1.block(r0,r0,r1,r1).template triangularView<UnitUpper>()  * m2.block(r0,c0,r1,c1), 1);

    // Zero temporaries for lazy products ...
    VERIFY_EVALUATION_COUNT( Scalar tmp = 0; tmp += Scalar(RealScalar(1)) /  (m3.transpose().lazyProduct(m3)).diagonal().sum(), 0 );

    // ... and even no temporary for even deeply (>=2) nested products
    VERIFY_EVALUATION_COUNT( Scalar tmp = 0; tmp += Scalar(RealScalar(1)) /  (m3.transpose() * m3).diagonal().sum(), 0 );
    VERIFY_EVALUATION_COUNT( Scalar tmp = 0; tmp += Scalar(RealScalar(1)) /  (m3.transpose() * m3).diagonal().array().abs().sum(), 0 );

    // Zero temporaries for ... CoeffBasedProductMode
    VERIFY_EVALUATION_COUNT( m3.col(0).template head<5>() * m3.col(0).transpose() + m3.col(0).template head<5>() * m3.col(0).transpose(), 0 );

//.........这里部分代码省略.........

template<typename MatrixType> void map_class_matrix(const MatrixType& m)
{
  typedef typename MatrixType::Index Index;
  typedef typename MatrixType::Scalar Scalar;

  Index rows = m.rows(), cols = m.cols(), size = rows*cols;
  Scalar s1 = internal::random<Scalar>();

  // array1 and array2 -> aligned heap allocation
  Scalar* array1 = internal::aligned_new<Scalar>(size);
  for(int i = 0; i < size; i++) array1[i] = Scalar(1);
  Scalar* array2 = internal::aligned_new<Scalar>(size);
  for(int i = 0; i < size; i++) array2[i] = Scalar(1);
  // array3unaligned -> unaligned pointer to heap
  Scalar* array3 = new Scalar[size+1];
  for(int i = 0; i < size+1; i++) array3[i] = Scalar(1);
  Scalar* array3unaligned = internal::UIntPtr(array3)%EIGEN_MAX_ALIGN_BYTES == 0 ? array3+1 : array3;
  Scalar array4[256];
  if(size<=256)
    for(int i = 0; i < size; i++) array4[i] = Scalar(1);
  
  Map<MatrixType> map1(array1, rows, cols);
  Map<MatrixType, AlignedMax> map2(array2, rows, cols);
  Map<MatrixType> map3(array3unaligned, rows, cols);
  Map<MatrixType> map4(array4, rows, cols);
  
  VERIFY_IS_EQUAL(map1, MatrixType::Ones(rows,cols));
  VERIFY_IS_EQUAL(map2, MatrixType::Ones(rows,cols));
  VERIFY_IS_EQUAL(map3, MatrixType::Ones(rows,cols));
  map1 = MatrixType::Random(rows,cols);
  map2 = map1;
  map3 = map1;
  MatrixType ma1 = map1;
  MatrixType ma2 = map2;
  MatrixType ma3 = map3;
  VERIFY_IS_EQUAL(map1, map2);
  VERIFY_IS_EQUAL(map1, map3);
  VERIFY_IS_EQUAL(ma1, ma2);
  VERIFY_IS_EQUAL(ma1, ma3);
  VERIFY_IS_EQUAL(ma1, map3);
  
  VERIFY_IS_APPROX(s1*map1, s1*map2);
  VERIFY_IS_APPROX(s1*ma1, s1*ma2);
  VERIFY_IS_EQUAL(s1*ma1, s1*ma3);
  VERIFY_IS_APPROX(s1*map1, s1*map3);
  
  map2 *= s1;
  map3 *= s1;
  VERIFY_IS_APPROX(s1*map1, map2);
  VERIFY_IS_APPROX(s1*map1, map3);
  
  if(size<=256)
  {
    VERIFY_IS_EQUAL(map4, MatrixType::Ones(rows,cols));
    map4 = map1;
    MatrixType ma4 = map4;
    VERIFY_IS_EQUAL(map1, map4);
    VERIFY_IS_EQUAL(ma1, map4);
    VERIFY_IS_EQUAL(ma1, ma4);
    VERIFY_IS_APPROX(s1*map1, s1*map4);
    
    map4 *= s1;
    VERIFY_IS_APPROX(s1*map1, map4);
  }

  internal::aligned_delete(array1, size);
  internal::aligned_delete(array2, size);
  delete[] array3;
}

template<typename MatrixType> void inverse(const MatrixType& m)
{
  using std::abs;
  typedef typename MatrixType::Index Index;
  /* this test covers the following files:
     Inverse.h
  */
  Index rows = m.rows();
  Index cols = m.cols();

  typedef typename MatrixType::Scalar Scalar;
  typedef typename NumTraits<Scalar>::Real RealScalar;
  typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;

  MatrixType m1(rows, cols),
             m2(rows, cols),
             identity = MatrixType::Identity(rows, rows);
  createRandomPIMatrixOfRank(rows,rows,rows,m1);
  m2 = m1.inverse();
  VERIFY_IS_APPROX(m1, m2.inverse() );

  VERIFY_IS_APPROX((Scalar(2)*m2).inverse(), m2.inverse()*Scalar(0.5));

  VERIFY_IS_APPROX(identity, m1.inverse() * m1 );
  VERIFY_IS_APPROX(identity, m1 * m1.inverse() );

  VERIFY_IS_APPROX(m1, m1.inverse().inverse() );

  // since for the general case we implement separately row-major and col-major, test that
  VERIFY_IS_APPROX(MatrixType(m1.transpose().inverse()), MatrixType(m1.inverse().transpose()));

#if !defined(EIGEN_TEST_PART_5) && !defined(EIGEN_TEST_PART_6)
  //computeInverseAndDetWithCheck tests
  //First: an invertible matrix
  bool invertible;
  RealScalar det;

  m2.setZero();
  m1.computeInverseAndDetWithCheck(m2, det, invertible);
  VERIFY(invertible);
  VERIFY_IS_APPROX(identity, m1*m2);
  VERIFY_IS_APPROX(det, m1.determinant());

  m2.setZero();
  m1.computeInverseWithCheck(m2, invertible);
  VERIFY(invertible);
  VERIFY_IS_APPROX(identity, m1*m2);

  //Second: a rank one matrix (not invertible, except for 1x1 matrices)
  VectorType v3 = VectorType::Random(rows);
  MatrixType m3 = v3*v3.transpose(), m4(rows,cols);
  m3.computeInverseAndDetWithCheck(m4, det, invertible);
  VERIFY( rows==1 ? invertible : !invertible );
  VERIFY_IS_MUCH_SMALLER_THAN(abs(det-m3.determinant()), RealScalar(1));
  m3.computeInverseWithCheck(m4, invertible);
  VERIFY( rows==1 ? invertible : !invertible );
#endif

  // check in-place inversion
  if(MatrixType::RowsAtCompileTime>=2 && MatrixType::RowsAtCompileTime<=4)
  {
    // in-place is forbidden
    VERIFY_RAISES_ASSERT(m1 = m1.inverse());
  }
  else
  {
    m2 = m1.inverse();
    m1 = m1.inverse();
    VERIFY_IS_APPROX(m1,m2);
  }
}

template<typename MatrixType> void adjoint(const MatrixType& m)
{
    /* this test covers the following files:
       Transpose.h Conjugate.h Dot.h
    */

    typedef typename MatrixType::Scalar Scalar;
    typedef typename NumTraits<Scalar>::Real RealScalar;
    typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
    typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> SquareMatrixType;
    int rows = m.rows();
    int cols = m.cols();

    RealScalar largerEps = test_precision<RealScalar>();
    if (ei_is_same_type<RealScalar,float>::ret)
        largerEps = RealScalar(1e-3f);

    MatrixType m1 = MatrixType::Random(rows, cols),
               m2 = MatrixType::Random(rows, cols),
               m3(rows, cols),
               mzero = MatrixType::Zero(rows, cols),
               identity = SquareMatrixType::Identity(rows, rows),
               square = SquareMatrixType::Random(rows, rows);
    VectorType v1 = VectorType::Random(rows),
               v2 = VectorType::Random(rows),
               v3 = VectorType::Random(rows),
               vzero = VectorType::Zero(rows);

    Scalar s1 = ei_random<Scalar>(),
           s2 = ei_random<Scalar>();

    // check basic compatibility of adjoint, transpose, conjugate
    VERIFY_IS_APPROX(m1.transpose().conjugate().adjoint(),    m1);
    VERIFY_IS_APPROX(m1.adjoint().conjugate().transpose(),    m1);

    // check multiplicative behavior
    VERIFY_IS_APPROX((m1.adjoint() * m2).adjoint(),           m2.adjoint() * m1);
    VERIFY_IS_APPROX((s1 * m1).adjoint(),                     ei_conj(s1) * m1.adjoint());

    // check basic properties of dot, norm, norm2
    typedef typename NumTraits<Scalar>::Real RealScalar;
    VERIFY(ei_isApprox((s1 * v1 + s2 * v2).dot(v3),   s1 * v1.dot(v3) + s2 * v2.dot(v3), largerEps));
    VERIFY(ei_isApprox(v3.dot(s1 * v1 + s2 * v2),     ei_conj(s1)*v3.dot(v1)+ei_conj(s2)*v3.dot(v2), largerEps));
    VERIFY_IS_APPROX(ei_conj(v1.dot(v2)),               v2.dot(v1));
    VERIFY_IS_APPROX(ei_abs(v1.dot(v1)),                v1.squaredNorm());
    if(NumTraits<Scalar>::HasFloatingPoint)
        VERIFY_IS_APPROX(v1.squaredNorm(),                      v1.norm() * v1.norm());
    VERIFY_IS_MUCH_SMALLER_THAN(ei_abs(vzero.dot(v1)),  static_cast<RealScalar>(1));
    if(NumTraits<Scalar>::HasFloatingPoint)
        VERIFY_IS_MUCH_SMALLER_THAN(vzero.norm(),         static_cast<RealScalar>(1));

    // check compatibility of dot and adjoint
    VERIFY(ei_isApprox(v1.dot(square * v2), (square.adjoint() * v1).dot(v2), largerEps));

    // like in testBasicStuff, test operator() to check const-qualification
    int r = ei_random<int>(0, rows-1),
        c = ei_random<int>(0, cols-1);
    VERIFY_IS_APPROX(m1.conjugate()(r,c), ei_conj(m1(r,c)));
    VERIFY_IS_APPROX(m1.adjoint()(c,r), ei_conj(m1(r,c)));

    if(NumTraits<Scalar>::HasFloatingPoint)
    {
        // check that Random().normalized() works: tricky as the random xpr must be evaluated by
        // normalized() in order to produce a consistent result.
        VERIFY_IS_APPROX(VectorType::Random(rows).normalized().norm(), RealScalar(1));
    }

    // check inplace transpose
    m3 = m1;
    m3.transposeInPlace();
    VERIFY_IS_APPROX(m3,m1.transpose());
    m3.transposeInPlace();
    VERIFY_IS_APPROX(m3,m1);

}

//************************************************************************************
//************************************************************************************
void UpdatedLagrangianFluid3Dinc::CalculateLocalSystem(MatrixType& rLeftHandSideMatrix, VectorType& rRightHandSideVector, ProcessInfo& rCurrentProcessInfo)
{
    KRATOS_TRY

    const double& density = 0.25*(GetGeometry()[0].FastGetSolutionStepValue(DENSITY)+
                                  GetGeometry()[1].FastGetSolutionStepValue(DENSITY) +
                                  GetGeometry()[2].FastGetSolutionStepValue(DENSITY)+
                                  GetGeometry()[3].FastGetSolutionStepValue(DENSITY));

    double K = 0.25* (GetGeometry()[0].FastGetSolutionStepValue(BULK_MODULUS)+
                      GetGeometry()[1].FastGetSolutionStepValue(BULK_MODULUS) +
                      GetGeometry()[2].FastGetSolutionStepValue(BULK_MODULUS)+
                      GetGeometry()[3].FastGetSolutionStepValue(BULK_MODULUS));
    K *= density;

    //unsigned int dim = 3;
    unsigned int number_of_nodes = 4;

    if(rLeftHandSideMatrix.size1() != 12)
        rLeftHandSideMatrix.resize(12,12,false);

    if(rRightHandSideVector.size() != 12)
        rRightHandSideVector.resize(12,false);

    //calculate current area
    double current_volume;
    GeometryUtils::CalculateGeometryData(GetGeometry(), msDN_Dx, msN, current_volume);


    //writing the body force
    const array_1d<double,3>& body_force = 0.25*(GetGeometry()[0].FastGetSolutionStepValue(BODY_FORCE)+
                                           GetGeometry()[1].FastGetSolutionStepValue(BODY_FORCE) +
                                           GetGeometry()[2].FastGetSolutionStepValue(BODY_FORCE) +
                                           GetGeometry()[3].FastGetSolutionStepValue(BODY_FORCE));

    for(unsigned int i = 0; i<number_of_nodes; i++)
    {
        rRightHandSideVector[i*3] = body_force[0]* density * mA0 * 0.25;
        rRightHandSideVector[i*3+1] = body_force[1] * density * mA0 * 0.25;
        rRightHandSideVector[i*3+2] = body_force[2] * density * mA0 * 0.25;

    }
    /*
    		//VISCOUS CONTRIBUTION TO THE STIFFNESS MATRIX
    		// rLeftHandSideMatrix += Laplacian * nu;
    		//filling matrix B
    		for (unsigned int i=0;i<number_of_nodes;i++)
    		{
    			unsigned int index = dim*i;
    			msB(0,index+0)=msDN_Dx(i,0);					msB(0,index+1)= 0.0;
    			msB(1,index+0)=0.0;								msB(1,index+1)= msDN_Dx(i,1);
    			msB(2,index+0)= msDN_Dx(i,1);					msB(2,index+1)= msDN_Dx(i,0);
    		}

    		//constitutive tensor
    		ms_constitutive_matrix(0,0) = K;	ms_constitutive_matrix(0,1) = K ;	ms_constitutive_matrix(0,2) = 0.0;
    		ms_constitutive_matrix(1,0) = K; 	ms_constitutive_matrix(1,1) = K;	ms_constitutive_matrix(1,2) = 0.0;
    		ms_constitutive_matrix(2,0) = 0.0;	ms_constitutive_matrix(2,1) = 0.0;	ms_constitutive_matrix(2,2) = 0.0;

    		//calculating viscous contributions
    		ms_temp = prod( ms_constitutive_matrix , msB);
    		noalias(rLeftHandSideMatrix) = prod( trans(msB) , ms_temp);
    		rLeftHandSideMatrix *= -current_area;
    */
    noalias(rLeftHandSideMatrix) = ZeroMatrix(12,12);
    //get the nodal pressure
    double p0 = GetGeometry()[0].FastGetSolutionStepValue(PRESSURE);
    double p1 = GetGeometry()[1].FastGetSolutionStepValue(PRESSURE);
    double p2 = GetGeometry()[2].FastGetSolutionStepValue(PRESSURE);
    double p3 = GetGeometry()[3].FastGetSolutionStepValue(PRESSURE);

    //adding pressure gradient
    double pavg = p0 + p1 + p2 + p3; //calculate old pressure over the element
    pavg *= 0.25 * current_volume;


    rRightHandSideVector[0] += msDN_Dx(0,0)*pavg;
    rRightHandSideVector[1] += msDN_Dx(0,1)*pavg;
    rRightHandSideVector[2] += msDN_Dx(0,2)*pavg;

    rRightHandSideVector[3] += msDN_Dx(1,0)*pavg;
    rRightHandSideVector[4] += msDN_Dx(1,1)*pavg;
    rRightHandSideVector[5] += msDN_Dx(1,2)*pavg;

    rRightHandSideVector[6] += msDN_Dx(2,0)*pavg;
    rRightHandSideVector[7] += msDN_Dx(2,1)*pavg;
    rRightHandSideVector[8] += msDN_Dx(2,2)*pavg;

    rRightHandSideVector[9]	 += msDN_Dx(3,0)*pavg;
    rRightHandSideVector[10] += msDN_Dx(3,1)*pavg;
    rRightHandSideVector[11] += msDN_Dx(3,2)*pavg;

    KRATOS_CATCH("")
}

//************************************************************************************
//************************************************************************************
void Monolithic2DNeumann::CalculateLocalSystem(MatrixType& rLeftHandSideMatrix, VectorType& rRightHandSideVector, ProcessInfo& rCurrentProcessInfo)
{

    if(rLeftHandSideMatrix.size1() != 4)
    {
        rLeftHandSideMatrix.resize(4,4,false);
        rRightHandSideVector.resize(4,false);

    }

    noalias(rLeftHandSideMatrix) = ZeroMatrix(4,4);

    //calculate normal to element.(normal follows the cross rule)
    array_1d<double,2> An,edge;
    edge[0] = GetGeometry()[1].X() - GetGeometry()[0].X();
    edge[1] = GetGeometry()[1].Y() - GetGeometry()[0].Y();


    double norm = edge[0]*edge[0] + edge[1]*edge[1];
    norm = pow(norm,0.5);

    An[0] = -edge[1];
    An[1] = edge[0];
    //An /= norm; this is then simplified by length of element in integration so is not divided.

    double mean_ex_p = 0.0;

    for(unsigned int i = 0; i<2 ; i++)
        mean_ex_p += 0.5*GetGeometry()[i].FastGetSolutionStepValue(EXTERNAL_PRESSURE);

    double p0 = GetGeometry()[0].FastGetSolutionStepValue(EXTERNAL_PRESSURE);
    rRightHandSideVector[0] = -0.5*An[0]*p0;
    rRightHandSideVector[1] = -0.5*An[1]*p0;

    double p1 = GetGeometry()[1].FastGetSolutionStepValue(EXTERNAL_PRESSURE);
    rRightHandSideVector[2] = -0.5*An[0]*p1;
    rRightHandSideVector[3] = -0.5*An[1]*p1;

    //	if(mean_ex_p !=GetGeometry()[0].FastGetSolutionStepValue(EXTERNAL_PRESSURE))
    //		mean_ex_p = 0.0;
    //KRATOS_WATCH(mean_ex_p);

    /*			for(unsigned int ii = 0; ii< 2; ++ii)
    				{

    					int id = (2 + 1)*(ii);
    					rRightHandSideVector[id] = mean_ex_p * An[0]* 0.5;
    					rRightHandSideVector[id + 1] = mean_ex_p * An[1]* 0.5;
    					rRightHandSideVector[id + 2] = 0.0;
    				//KRATOS_WATCH(An);
    				}*/
// 			KRATOS_WATCH(An);
//KRATOS_WATCH(p0);
//KRATOS_WATCH(p1);
//KRATOS_WATCH("TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT");
    //	KRATOS_WATCH(rRightHandSideVector);




}

typename GaussianProcess<TScalarType>::MatrixType GaussianProcess<TScalarType>::InvertKernelMatrix(const typename GaussianProcess<TScalarType>::MatrixType &K,
                                                      typename GaussianProcess<TScalarType>::InversionMethod inv_method = GaussianProcess<TScalarType>::FullPivotLU,
                                                                                                   bool stable) const{
    // compute core matrix
    if(debug){
        std::cout << "GaussianProcess::InvertKernelMatrix: inverting kernel matrix... ";
        std::cout.flush();
    }

    typename GaussianProcess<TScalarType>::MatrixType core;

    switch(inv_method){
    // standard method: fast but not that accurate
    // Uses the LU decomposition with full pivoting for the inversion
    case FullPivotLU:{
        if(debug) std::cout << " (inversion method: FullPivotLU) " << std::flush;
        try{
            if(stable){
                core = K.inverse();
            }
            else{
                if(debug) std::cout << " (using lapack) " << std::flush;
                core = lapack::lu_invert<TScalarType>(K);
            }
        }
        catch(lapack::LAPACKException& e){
            core = K.inverse();
        }
    }
    break;

    // very accurate and very slow method, use it for small problems
    // Uses the two-sided Jacobi SVD decomposition
    case JacobiSVD:{
        if(debug) std::cout << " (inversion method: JacobiSVD) " << std::flush;
        Eigen::JacobiSVD<MatrixType> jacobisvd(K, Eigen::ComputeThinU | Eigen::ComputeThinV);
        if((jacobisvd.singularValues().real().array() < 0).any() && debug){
            std::cout << "GaussianProcess::InvertKernelMatrix: warning: there are negative eigenvalues.";
            std::cout.flush();
        }
        core = jacobisvd.matrixV() * VectorType(1/jacobisvd.singularValues().array()).asDiagonal() * jacobisvd.matrixU().transpose();
    }
    break;

    // accurate method and faster than Jacobi SVD.
    // Uses the bidiagonal divide and conquer SVD
    case BDCSVD:{
        if(debug) std::cout << " (inversion method: BDCSVD) " << std::flush;
#ifdef EIGEN_BDCSVD_H
        Eigen::BDCSVD<MatrixType> bdcsvd(K, Eigen::ComputeThinU | Eigen::ComputeThinV);
        if((bdcsvd.singularValues().real().array() < 0).any() && debug){
            std::cout << "GaussianProcess::InvertKernelMatrix: warning: there are negative eigenvalues.";
            std::cout.flush();
        }
        core = bdcsvd.matrixV() * VectorType(1/bdcsvd.singularValues().array()).asDiagonal() * bdcsvd.matrixU().transpose();
#else
        // this is checked, since BDCSVD is currently not in the newest release
        throw std::string("GaussianProcess::InvertKernelMatrix: BDCSVD is not supported by the provided Eigen library.");
#endif

    }
    break;

    // faster than the SVD method but less stable
    // computes the eigenvalues/eigenvectors of selfadjoint matrices
    case SelfAdjointEigenSolver:{
        if(debug) std::cout << " (inversion method: SelfAdjointEigenSolver) " << std::flush;
        try{
            core = lapack::chol_invert<TScalarType>(K);
        }
        catch(lapack::LAPACKException& e){
            Eigen::SelfAdjointEigenSolver<MatrixType> es;
            es.compute(K);
            VectorType eigenValues = es.eigenvalues().reverse();
            MatrixType eigenVectors = es.eigenvectors().rowwise().reverse();
            if((eigenValues.real().array() < 0).any() && debug){
                std::cout << "GaussianProcess::InvertKernelMatrix: warning: there are negative eigenvalues.";
                std::cout.flush();
            }
            core = eigenVectors * VectorType(1/eigenValues.array()).asDiagonal() * eigenVectors.transpose();
        }
    }
    break;
    }

    if(debug) std::cout << "[done]" << std::endl;
    return core;
}

void GaussianProcess<TScalarType>::Load(std::string prefix){
    if(debug){
        std::cout << "GaussianProcess::Load: loading gaussian process: " << std::endl;
        std::cout << "\t " << prefix+"-RegressionVectors.txt" << std::endl;
        std::cout << "\t " << prefix+"-CoreMatrix.txt" << std::endl;
        std::cout << "\t " << prefix+"-SampleVectors.txt" << std::endl;
        std::cout << "\t " << prefix+"-LabelVectors.txt" << std::endl;
        std::cout << "\t " << prefix+"-ParameterFile.txt" << std::endl;
    }

    // load regression vectors
    std::string rv_filename = prefix+"-RegressionVectors.txt";
    fs::path rv_file(rv_filename.c_str());
    if(!(fs::exists(rv_file) && !fs::is_directory(rv_file))){
        throw std::string("GaussianProcess::Load: "+rv_filename+" does not exist or is a directory.");
    }
    m_RegressionVectors = ReadMatrix<MatrixType>(rv_filename);

    // load core matrix
    std::string cm_filename = prefix+"-CoreMatrix.txt";
    fs::path cm_file(cm_filename.c_str());
    if(!(fs::exists(cm_file) && !fs::is_directory(cm_file))){
        throw std::string("GaussianProcess::Load: "+cm_filename+" does not exist or is a directory.");
    }
    m_CoreMatrix = ReadMatrix<MatrixType>(cm_filename);

    // load sample vectors
    std::string sv_filename = prefix+"-SampleVectors.txt";
    fs::path sv_file(sv_filename.c_str());
    if(!(fs::exists(sv_file) && !fs::is_directory(sv_file))){
        throw std::string("GaussianProcess::Load: "+sv_filename+" does not exist or is a directory.");
    }
    MatrixType X = ReadMatrix<MatrixType>(sv_filename);
    m_SampleVectors.clear();
    for(unsigned i=0; i<X.cols(); i++){
        m_SampleVectors.push_back(X.col(i));
    }

    // load label vectors
    std::string lv_filename = prefix+"-LabelVectors.txt";
    fs::path lv_file(lv_filename.c_str());
    if(!(fs::exists(lv_file) && !fs::is_directory(lv_file))){
        throw std::string("GaussianProcess::Load: "+lv_filename+" does not exist or is a directory.");
    }
    MatrixType Y = ReadMatrix<MatrixType>(lv_filename);
    m_LabelVectors.clear();
    for(unsigned i=0; i<Y.cols(); i++){
        m_LabelVectors.push_back(Y.col(i));
    }

    // load parameters
    std::string pf_filename = prefix+"-ParameterFile.txt";
    if(!(fs::exists(pf_filename) && !fs::is_directory(pf_filename))){
        throw std::string("GaussianProcess::Load: "+pf_filename+" does not exist or is a directory.");
    }

    std::ifstream parameter_infile;
    parameter_infile.open( pf_filename.c_str() );

    // reading parameter file
    std::string line;
    bool load = true;
    if(std::getline(parameter_infile, line)) {
        std::stringstream line_stream(line);

        if(!(line_stream >> m_Sigma &&
             line_stream >> m_InputDimension &&
             line_stream >> m_OutputDimension &&
             line_stream >> m_EfficientStorage &&
             line_stream >> debug) ||
                load == false){
            throw std::string("GaussianProcess::Load: parameter file is corrupt");
        }

        std::string kernel_string;
        line_stream >> kernel_string;

        m_Kernel = KernelFactory<TScalarType>::GetKernel(kernel_string);

    }
    parameter_infile.close();


    m_Initialized = true;
}

//----------------------
//-----  PRIVATE  ------
//----------------------
//***********************************************************************************
void FaceHeatConvection::CalculateAll( MatrixType& rLeftHandSideMatrix, VectorType& rRightHandSideVector, const ProcessInfo& rCurrentProcessInfo, bool CalculateStiffnessMatrixFlag, bool CalculateResidualVectorFlag )
{
    KRATOS_TRY

    const unsigned int number_of_nodes = GetGeometry().size();
    unsigned int MatSize = number_of_nodes;

    //resizing as needed the LHS
    if ( CalculateStiffnessMatrixFlag == true ) //calculation of the matrix is required
    {
        if ( rLeftHandSideMatrix.size1() != MatSize )
            rLeftHandSideMatrix.resize( MatSize, MatSize, false );
        noalias( rLeftHandSideMatrix ) = ZeroMatrix( MatSize, MatSize ); //resetting LHS
    }

    //resizing as needed the RHS
    if ( CalculateResidualVectorFlag == true ) //calculation of the matrix is required
    {
        if ( rRightHandSideVector.size() != MatSize )
            rRightHandSideVector.resize( MatSize );
        rRightHandSideVector = ZeroVector( MatSize ); //resetting RHS
    }

    //reading integration points and local gradients
    const GeometryType::IntegrationPointsArrayType& integration_points = GetGeometry().IntegrationPoints();
    const GeometryType::ShapeFunctionsGradientsType& DN_DeContainer = GetGeometry().ShapeFunctionsLocalGradients();
    const Matrix& Ncontainer = GetGeometry().ShapeFunctionsValues();

    //calculating actual jacobian
    GeometryType::JacobiansType J;
    J = GetGeometry().Jacobian( J );

    //auxiliary terms
    for ( unsigned int PointNumber = 0; PointNumber < integration_points.size(); PointNumber++ )
    {
        double convection_coefficient = 0.0;
        double T_ambient = 0.0;
        double T = 0.0;


        for ( unsigned int n = 0; n < GetGeometry().size(); n++ )
        {
            convection_coefficient += ( GetGeometry()[n] ).GetSolutionStepValue( CONVECTION_COEFFICIENT ) * Ncontainer( PointNumber, n );
            T_ambient += ( GetGeometry()[n] ).GetSolutionStepValue( AMBIENT_TEMPERATURE ) * Ncontainer( PointNumber, n );
            T += ( GetGeometry()[n] ).GetSolutionStepValue( TEMPERATURE ) * Ncontainer( PointNumber, n );
        }

//         if ( PointNumber == 1 )
//             std::cout << "CONDITION ### HeatConvection:  h= " << convection_coefficient << ",\t T_ambient= " << T_ambient << ",\t T_surface= " << T << std::endl;

        double IntegrationWeight = GetGeometry().IntegrationPoints()[PointNumber].Weight();
        Vector t1 = ZeroVector( 3 );//first tangential vector
        Vector t2 = ZeroVector( 3 );//second tangential vector

        for ( unsigned int n = 0; n < number_of_nodes; n++ )
        {
            t1[0] += GetGeometry().GetPoint( n ).X0() * DN_DeContainer[PointNumber]( n, 0 );
            t1[1] += GetGeometry().GetPoint( n ).Y0() * DN_DeContainer[PointNumber]( n, 0 );
            t1[2] += GetGeometry().GetPoint( n ).Z0() * DN_DeContainer[PointNumber]( n, 0 );
            t2[0] += GetGeometry().GetPoint( n ).X0() * DN_DeContainer[PointNumber]( n, 1 );
            t2[1] += GetGeometry().GetPoint( n ).Y0() * DN_DeContainer[PointNumber]( n, 1 );
            t2[2] += GetGeometry().GetPoint( n ).Z0() * DN_DeContainer[PointNumber]( n, 1 );
        }

        //calculating normal
        Vector v3 = ZeroVector( 3 );
        v3[0] = t1[1] * t2[2] - t1[2] * t2[1];
        v3[1] = t1[2] * t2[0] - t1[0] * t2[2];
        v3[2] = t1[0] * t2[1] - t1[1] * t2[0];

        double	dA = sqrt( v3[0] * v3[0] + v3[1] * v3[1] + v3[2] * v3[2] );

        if ( CalculateResidualVectorFlag == true ) //calculation of the matrix is required
        {
            for ( unsigned int n = 0; n < number_of_nodes; n++ )
                rRightHandSideVector( n ) -= Ncontainer( PointNumber, n )
                                             * convection_coefficient * ( T - T_ambient )
                                             * IntegrationWeight * dA; // W/(m^2*°C) = kg*s^-3*°C°^-1
        }

        if ( CalculateStiffnessMatrixFlag == true ) //calculation of the matrix is required
        {
            for ( unsigned int n = 0; n < number_of_nodes; n++ )
                rLeftHandSideMatrix( n, n ) += Ncontainer( PointNumber, n )
                                               * convection_coefficient
                                               * Ncontainer( PointNumber, n ) * IntegrationWeight * dA; // W/(m^2*°C) = kg*s^-3*°C°^-1
        }

    }
    KRATOS_CATCH( "" )
}

 static void run(MatrixType& result, typename MatrixType::Index size)
 {
   result.resize(size, size);
   result.template triangularView<Upper>() = MatrixType::Random(size, size);
 }