text: // This file is part of Eigen, a lightweight C++ template library
  // for linear algebra.
  //
  // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
  //
  // 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/.
  
  #ifndef EIGEN_SPARSETRIANGULARSOLVER_H
  #define EIGEN_SPARSETRIANGULARSOLVER_H
  
  namespace Eigen { 
  
  namespace internal {
  
  template<typename Lhs, typename Rhs, int Mode,
    int UpLo = (Mode & Lower)
             ? Lower
             : (Mode & Upper)
             ? Upper
             : -1,
    int StorageOrder = int(traits<Lhs>::Flags) & RowMajorBit>
  struct sparse_solve_triangular_selector;
  
  // forward substitution, row-major
  template<typename Lhs, typename Rhs, int Mode>
  struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Lower,RowMajor>
  {
    typedef typename Rhs::Scalar Scalar;
    typedef evaluator<Lhs> LhsEval;
    typedef typename evaluator<Lhs>::InnerIterator LhsIterator;
    static void run(const Lhs& lhs, Rhs& other)
    {
      LhsEval lhsEval(lhs);
      for(Index col=0 ; col<other.cols() ; ++col)
      {
        for(Index i=0; i<lhs.rows(); ++i)
        {
          Scalar tmp = other.coeff(i,col);
          Scalar lastVal(0);
          Index lastIndex = 0;
          for(LhsIterator it(lhsEval, i); it; ++it)
          {
            lastVal = it.value();
            lastIndex = it.index();
            if(lastIndex==i)
              break;
            tmp -= lastVal * other.coeff(lastIndex,col);
          }
          if (Mode & UnitDiag)
            other.coeffRef(i,col) = tmp;
          else
          {
            eigen_assert(lastIndex==i);
            other.coeffRef(i,col) = tmp/lastVal;
          }
        }
      }
    }
  };
  
  // backward substitution, row-major
  template<typename Lhs, typename Rhs, int Mode>
  struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Upper,RowMajor>
  {
    typedef typename Rhs::Scalar Scalar;
    typedef evaluator<Lhs> LhsEval;
    typedef typename evaluator<Lhs>::InnerIterator LhsIterator;
    static void run(const Lhs& lhs, Rhs& other)
    {
      LhsEval lhsEval(lhs);
      for(Index col=0 ; col<other.cols() ; ++col)
      {
        for(Index i=lhs.rows()-1 ; i>=0 ; --i)
        {
          Scalar tmp = other.coeff(i,col);
          Scalar l_ii(0);
          LhsIterator it(lhsEval, i);
          while(it && it.index()<i)
            ++it;
          if(!(Mode & UnitDiag))
          {
            eigen_assert(it && it.index()==i);
            l_ii = it.value();
            ++it;
          }
          else if (it && it.index() == i)
            ++it;
          for(; it; ++it)
          {
            tmp -= it.value() * other.coeff(it.index(),col);
          }
  
          if (Mode & UnitDiag)  other.coeffRef(i,col) = tmp;
          else                  other.coeffRef(i,col) = tmp/l_ii;
        }
      }
    }
  };
  
  // forward substitution, col-major
  template<typename Lhs, typename Rhs, int Mode>
  struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Lower,ColMajor>
  {
    typedef typename Rhs::Scalar Scalar;
    typedef evaluator<Lhs> LhsEval;
    typedef typename evaluator<Lhs>::InnerIterator LhsIterator;
    static void run(const Lhs& lhs, Rhs& other)
    {
      LhsEval lhsEval(lhs);
      for(Index col=0 ; col<other.cols() ; ++col)
      {
        for(Index i=0; i<lhs.cols(); ++i)
        {
          Scalar& tmp = other.coeffRef(i,col);
          if (tmp!=Scalar(0)) // optimization when other is actually sparse
          {
            LhsIterator it(lhsEval, i);
            while(it && it.index()<i)
              ++it;
            if(!(Mode & UnitDiag))
            {
              eigen_assert(it && it.index()==i);
              tmp /= it.value();
            }
            if (it && it.index()==i)
              ++it;
            for(; it; ++it)
              other.coeffRef(it.index(), col) -= tmp * it.value();
          }
        }
      }
    }
  };
  
  // backward substitution, col-major
  template<typename Lhs, typename Rhs, int Mode>
  struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Upper,ColMajor>
  {
    typedef typename Rhs::Scalar Scalar;
    typedef evaluator<Lhs> LhsEval;
    typedef typename evaluator<Lhs>::InnerIterator LhsIterator;
    static void run(const Lhs& lhs, Rhs& other)
    {
      LhsEval lhsEval(lhs);
      for(Index col=0 ; col<other.cols() ; ++col)
      {
        for(Index i=lhs.cols()-1; i>=0; --i)
        {
          Scalar& tmp = other.coeffRef(i,col);
          if (tmp!=Scalar(0)) // optimization when other is actually sparse
          {
            if(!(Mode & UnitDiag))
            {
              // TODO replace this by a binary search. make sure the binary search is safe for partially sorted elements
              LhsIterator it(lhsEval, i);
              while(it && it.index()!=i)
                ++it;
              eigen_assert(it && it.index()==i);
              other.coeffRef(i,col) /= it.value();
            }
            LhsIterator it(lhsEval, i);
            for(; it && it.index()<i; ++it)
              other.coeffRef(it.index(), col) -= tmp * it.value();
          }
        }
      }
    }
  };
  
  } // end namespace internal
  
  template<typename ExpressionType,unsigned int Mode>
  template<typename OtherDerived>
  void TriangularViewImpl<ExpressionType,Mode,Sparse>::solveInPlace(MatrixBase<OtherDerived>& other) const
  {
    eigen_assert(derived().cols() == derived().rows() && derived().cols() == other.rows());
    eigen_assert((!(Mode & ZeroDiag)) && bool(Mode & (Upper|Lower)));
  
    enum { copy = internal::traits<OtherDerived>::Flags & RowMajorBit };
  
    typedef typename internal::conditional<copy,
      typename internal::plain_matrix_type_column_major<OtherDerived>::type, OtherDerived&>::type OtherCopy;
    OtherCopy otherCopy(other.derived());
  
    internal::sparse_solve_triangular_selector<ExpressionType, typename internal::remove_reference<OtherCopy>::type, Mode>::run(derived().nestedExpression(), otherCopy);
  
    if (copy)
      other = otherCopy;
  }
  
  // pure sparse path
  
  namespace internal {
  
  template<typename Lhs, typename Rhs, int Mode,
    int UpLo = (Mode & Lower)
             ? Lower
             : (Mode & Upper)
             ? Upper
             : -1,
    int StorageOrder = int(Lhs::Flags) & (RowMajorBit)>
  struct sparse_solve_triangular_sparse_selector;
  
  // forward substitution, col-major
  template<typename Lhs, typename Rhs, int Mode, int UpLo>
  struct sparse_solve_triangular_sparse_selector<Lhs,Rhs,Mode,UpLo,ColMajor>
  {
    typedef typename Rhs::Scalar Scalar;
    typedef typename promote_index_type<typename traits<Lhs>::StorageIndex,
                                        typename traits<Rhs>::StorageIndex>::type StorageIndex;
    static void run(const Lhs& lhs, Rhs& other)
    {
      const bool IsLower = (UpLo==Lower);
      AmbiVector<Scalar,StorageIndex> tempVector(other.rows()*2);
      tempVector.setBounds(0,other.rows());
  
      Rhs res(other.rows(), other.cols());
      res.reserve(other.nonZeros());
  
      for(Index col=0 ; col<other.cols() ; ++col)
      {
        // FIXME estimate number of non zeros
        tempVector.init(.99/*float(other.col(col).nonZeros())/float(other.rows())*/);
        tempVector.setZero();
        tempVector.restart();
        for (typename Rhs::InnerIterator rhsIt(other, col); rhsIt; ++rhsIt)
        {
          tempVector.coeffRef(rhsIt.index()) = rhsIt.value();
        }
  
        for(Index i=IsLower?0:lhs.cols()-1;
            IsLower?i<lhs.cols():i>=0;
            i+=IsLower?1:-1)
        {
          tempVector.restart();
          Scalar& ci = tempVector.coeffRef(i);
          if (ci!=Scalar(0))
          {
            // find
            typename Lhs::InnerIterator it(lhs, i);
            if(!(Mode & UnitDiag))
            {
              if (IsLower)
              {
                eigen_assert(it.index()==i);
                ci /= it.value();
              }
              else
                ci /= lhs.coeff(i,i);
            }
            tempVector.restart();
            if (IsLower)
            {
              if (it.index()==i)
                ++it;
              for(; it; ++it)
                tempVector.coeffRef(it.index()) -= ci * it.value();
            }
            else
            {
              for(; it && it.index()<i; ++it)
                tempVector.coeffRef(it.index()) -= ci * it.value();
            }
          }
        }
  
  
        Index count = 0;
        // FIXME compute a reference value to filter zeros
        for (typename AmbiVector<Scalar,StorageIndex>::Iterator it(tempVector/*,1e-12*/); it; ++it)
        {
          ++ count;
  //         std::cerr << "fill " << it.index() << ", " << col << "\n";
  //         std::cout << it.value() << "  ";
          // FIXME use insertBack
          res.insert(it.index(), col) = it.value();
        }
  //       std::cout << "tempVector.nonZeros() == " << int(count) << " / " << (other.rows()) << "\n";
      }
      res.finalize();
      other = res.markAsRValue();
    }
  };
  
  } // end namespace internal
  
  template<typename ExpressionType,unsigned int Mode>
  template<typename OtherDerived>
  void TriangularViewImpl<ExpressionType,Mode,Sparse>::solveInPlace(SparseMatrixBase<OtherDerived>& other) const
  {
    eigen_assert(derived().cols() == derived().rows() && derived().cols() == other.rows());
    eigen_assert( (!(Mode & ZeroDiag)) && bool(Mode & (Upper|Lower)));
  
  //   enum { copy = internal::traits<OtherDerived>::Flags & RowMajorBit };
  
  //   typedef typename internal::conditional<copy,
  //     typename internal::plain_matrix_type_column_major<OtherDerived>::type, OtherDerived&>::type OtherCopy;
  //   OtherCopy otherCopy(other.derived());
  
    internal::sparse_solve_triangular_sparse_selector<ExpressionType, OtherDerived, Mode>::run(derived().nestedExpression(), other.derived());
  
  //   if (copy)
  //     other = otherCopy;
  }
  
  } // end namespace Eigen
  
  #endif // EIGEN_SPARSETRIANGULARSOLVER_H
  