1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr> 5c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 6c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This Source Code Form is subject to the terms of the Mozilla 7c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Public License v. 2.0. If a copy of the MPL was not distributed 8c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 10c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#ifndef EIGEN_SPARSETRIANGULARSOLVER_H 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_SPARSETRIANGULARSOLVER_H 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace Eigen { 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Lhs, typename Rhs, int Mode, 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int UpLo = (Mode & Lower) 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ? Lower 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : (Mode & Upper) 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ? Upper 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : -1, 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int StorageOrder = int(traits<Lhs>::Flags) & RowMajorBit> 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct sparse_solve_triangular_selector; 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// forward substitution, row-major 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Lhs, typename Rhs, int Mode> 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Lower,RowMajor> 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Rhs::Scalar Scalar; 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static void run(const Lhs& lhs, Rhs& other) 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int col=0 ; col<other.cols() ; ++col) 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<lhs.rows(); ++i) 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar tmp = other.coeff(i,col); 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar lastVal(0); 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int lastIndex = 0; 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(typename Lhs::InnerIterator it(lhs, i); it; ++it) 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lastVal = it.value(); 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lastIndex = it.index(); 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(lastIndex==i) 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath break; 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp -= lastVal * other.coeff(lastIndex,col); 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (Mode & UnitDiag) 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath other.coeffRef(i,col) = tmp; 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(lastIndex==i); 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath other.coeffRef(i,col) = tmp/lastVal; 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// backward substitution, row-major 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Lhs, typename Rhs, int Mode> 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Upper,RowMajor> 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Rhs::Scalar Scalar; 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static void run(const Lhs& lhs, Rhs& other) 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int col=0 ; col<other.cols() ; ++col) 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=lhs.rows()-1 ; i>=0 ; --i) 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar tmp = other.coeff(i,col); 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar l_ii = 0; 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename Lhs::InnerIterator it(lhs, i); 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath while(it && it.index()<i) 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++it; 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(!(Mode & UnitDiag)) 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(it && it.index()==i); 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath l_ii = it.value(); 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++it; 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if (it && it.index() == i) 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++it; 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(; it; ++it) 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp -= it.value() * other.coeff(it.index(),col); 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (Mode & UnitDiag) 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath other.coeffRef(i,col) = tmp; 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath other.coeffRef(i,col) = tmp/l_ii; 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// forward substitution, col-major 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Lhs, typename Rhs, int Mode> 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Lower,ColMajor> 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Rhs::Scalar Scalar; 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static void run(const Lhs& lhs, Rhs& other) 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int col=0 ; col<other.cols() ; ++col) 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<lhs.cols(); ++i) 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar& tmp = other.coeffRef(i,col); 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (tmp!=Scalar(0)) // optimization when other is actually sparse 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename Lhs::InnerIterator it(lhs, i); 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath while(it && it.index()<i) 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++it; 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(!(Mode & UnitDiag)) 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(it && it.index()==i); 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp /= it.value(); 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (it && it.index()==i) 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++it; 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(; it; ++it) 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath other.coeffRef(it.index(), col) -= tmp * it.value(); 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// backward substitution, col-major 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Lhs, typename Rhs, int Mode> 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Upper,ColMajor> 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Rhs::Scalar Scalar; 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static void run(const Lhs& lhs, Rhs& other) 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int col=0 ; col<other.cols() ; ++col) 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=lhs.cols()-1; i>=0; --i) 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar& tmp = other.coeffRef(i,col); 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (tmp!=Scalar(0)) // optimization when other is actually sparse 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(!(Mode & UnitDiag)) 145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // TODO replace this by a binary search. make sure the binary search is safe for partially sorted elements 147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename Lhs::ReverseInnerIterator it(lhs, i); 148c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath while(it && it.index()!=i) 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath --it; 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(it && it.index()==i); 151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath other.coeffRef(i,col) /= it.value(); 152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 153c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename Lhs::InnerIterator it(lhs, i); 154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(; it && it.index()<i; ++it) 155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath other.coeffRef(it.index(), col) -= tmp * it.value(); 156c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 157c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 158c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 159c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 160c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace internal 163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename ExpressionType,int Mode> 165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename OtherDerived> 166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid SparseTriangularView<ExpressionType,Mode>::solveInPlace(MatrixBase<OtherDerived>& other) const 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_matrix.cols() == m_matrix.rows() && m_matrix.cols() == other.rows()); 169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert((!(Mode & ZeroDiag)) && bool(Mode & (Upper|Lower))); 170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 171c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { copy = internal::traits<OtherDerived>::Flags & RowMajorBit }; 172c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename internal::conditional<copy, 174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename internal::plain_matrix_type_column_major<OtherDerived>::type, OtherDerived&>::type OtherCopy; 175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath OtherCopy otherCopy(other.derived()); 176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 177c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::sparse_solve_triangular_selector<ExpressionType, typename internal::remove_reference<OtherCopy>::type, Mode>::run(m_matrix, otherCopy); 178c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 179c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (copy) 180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath other = otherCopy; 181c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 182c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 183c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename ExpressionType,int Mode> 184c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename OtherDerived> 185c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtypename internal::plain_matrix_type_column_major<OtherDerived>::type 186c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathSparseTriangularView<ExpressionType,Mode>::solve(const MatrixBase<OtherDerived>& other) const 187c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 188c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename internal::plain_matrix_type_column_major<OtherDerived>::type res(other); 189c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath solveInPlace(res); 190c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return res; 191c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 192c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 193c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// pure sparse path 194c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 195c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 196c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 197c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Lhs, typename Rhs, int Mode, 198c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int UpLo = (Mode & Lower) 199c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ? Lower 200c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : (Mode & Upper) 201c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ? Upper 202c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : -1, 203c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int StorageOrder = int(Lhs::Flags) & (RowMajorBit)> 204c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct sparse_solve_triangular_sparse_selector; 205c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 206c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// forward substitution, col-major 207c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Lhs, typename Rhs, int Mode, int UpLo> 208c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct sparse_solve_triangular_sparse_selector<Lhs,Rhs,Mode,UpLo,ColMajor> 209c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 210c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Rhs::Scalar Scalar; 211c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename promote_index_type<typename traits<Lhs>::Index, 212c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename traits<Rhs>::Index>::type Index; 213c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static void run(const Lhs& lhs, Rhs& other) 214c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 215c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const bool IsLower = (UpLo==Lower); 216c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath AmbiVector<Scalar,Index> tempVector(other.rows()*2); 217c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tempVector.setBounds(0,other.rows()); 218c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 219c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Rhs res(other.rows(), other.cols()); 220c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath res.reserve(other.nonZeros()); 221c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 222c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int col=0 ; col<other.cols() ; ++col) 223c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 224c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // FIXME estimate number of non zeros 225c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tempVector.init(.99/*float(other.col(col).nonZeros())/float(other.rows())*/); 226c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tempVector.setZero(); 227c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tempVector.restart(); 228c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (typename Rhs::InnerIterator rhsIt(other, col); rhsIt; ++rhsIt) 229c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 230c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tempVector.coeffRef(rhsIt.index()) = rhsIt.value(); 231c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 232c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 233c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=IsLower?0:lhs.cols()-1; 234c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath IsLower?i<lhs.cols():i>=0; 235c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath i+=IsLower?1:-1) 236c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 237c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tempVector.restart(); 238c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar& ci = tempVector.coeffRef(i); 239c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (ci!=Scalar(0)) 240c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 241c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // find 242c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename Lhs::InnerIterator it(lhs, i); 243c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(!(Mode & UnitDiag)) 244c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 245c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (IsLower) 246c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 247c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(it.index()==i); 248c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ci /= it.value(); 249c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 250c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 251c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ci /= lhs.coeff(i,i); 252c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 253c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tempVector.restart(); 254c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (IsLower) 255c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 256c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (it.index()==i) 257c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++it; 258c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(; it; ++it) 259c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tempVector.coeffRef(it.index()) -= ci * it.value(); 260c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 261c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 262c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 263c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(; it && it.index()<i; ++it) 264c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tempVector.coeffRef(it.index()) -= ci * it.value(); 265c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 266c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 267c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 268c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 269c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 270c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int count = 0; 271c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // FIXME compute a reference value to filter zeros 272c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (typename AmbiVector<Scalar,Index>::Iterator it(tempVector/*,1e-12*/); it; ++it) 273c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 274c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++ count; 275c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// std::cerr << "fill " << it.index() << ", " << col << "\n"; 276c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// std::cout << it.value() << " "; 277c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // FIXME use insertBack 278c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath res.insert(it.index(), col) = it.value(); 279c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 280c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// std::cout << "tempVector.nonZeros() == " << int(count) << " / " << (other.rows()) << "\n"; 281c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 282c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath res.finalize(); 283c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath other = res.markAsRValue(); 284c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 285c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 286c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 287c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace internal 288c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 289c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename ExpressionType,int Mode> 290c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename OtherDerived> 291c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid SparseTriangularView<ExpressionType,Mode>::solveInPlace(SparseMatrixBase<OtherDerived>& other) const 292c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 293c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_matrix.cols() == m_matrix.rows() && m_matrix.cols() == other.rows()); 294c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert( (!(Mode & ZeroDiag)) && bool(Mode & (Upper|Lower))); 295c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 296c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// enum { copy = internal::traits<OtherDerived>::Flags & RowMajorBit }; 297c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 298c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// typedef typename internal::conditional<copy, 299c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// typename internal::plain_matrix_type_column_major<OtherDerived>::type, OtherDerived&>::type OtherCopy; 300c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// OtherCopy otherCopy(other.derived()); 301c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 302c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::sparse_solve_triangular_sparse_selector<ExpressionType, OtherDerived, Mode>::run(m_matrix, other.derived()); 303c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 304c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// if (copy) 305c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// other = otherCopy; 306c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 307c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 308c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#ifdef EIGEN2_SUPPORT 309c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 310c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// deprecated stuff: 311c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 312c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \deprecated */ 313c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Derived> 314c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename OtherDerived> 315c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid SparseMatrixBase<Derived>::solveTriangularInPlace(MatrixBase<OtherDerived>& other) const 316c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 317c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath this->template triangular<Flags&(Upper|Lower)>().solveInPlace(other); 318c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 319c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 320c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \deprecated */ 321c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Derived> 322c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename OtherDerived> 323c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtypename internal::plain_matrix_type_column_major<OtherDerived>::type 324c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathSparseMatrixBase<Derived>::solveTriangular(const MatrixBase<OtherDerived>& other) const 325c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 326c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename internal::plain_matrix_type_column_major<OtherDerived>::type res(other); 327c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath derived().solveTriangularInPlace(res); 328c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return res; 329c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 330c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN2_SUPPORT 331c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 332c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 333c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 334c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_SPARSETRIANGULARSOLVER_H 335