From: Karl Meerbergen (Karl.Meerbergen_at_[hidden])
Date: 2006-12-15 08:59:24
Gunter,
I have a suggestion: instead of returning a bool you could return a
X::size_type to indicate how far the triangular solve is carried out.
So, if x is a vector and the return value is x.size(), then the solve is
fully carried out. Otherwise the return value tells you how many
elements of x have been solved for.
Karl
Gunter Winkler wrote:
>Hallo,
>
>I finally finished the new triangular solvers. A new header file and an
>example is attached. Some documentation is inside the header file. These
>solves are defined as a replacement of the solvers in triangular.hpp (which
>must be disabled in order to compile the example.)
>
>I found it to be convenient to gain a triangular type tag from a triangular
>layout type. So I added a few lines to functional.hpp (see patch file). Now
>you can use
>typedef lower TRI;
>typedef TRI::triangular_type TAG;
>where TAG is lower_tag (which can be used as an function argument of
>inplace_solve.)
>This is similar to orientation_category in row_major and column_major.
>
>The syntax of the solvers is the same:
>
>bool singular;
>singular = inplace_solve(L, x, lower_tag); // solves Lx=b
>singular = inplace_solve(L, x, unit_lower_tag); // solves (I+L)x=b
>singular = inplace_solve(U, x, upper_tag); // solves Ux=b
>singular = inplace_solve(U, x, unit_upper_tag); // solves (I+U)x=b
>
>(Most of) The solvers require the matrix to have the given structure and they
>may return wrong results if this condition is violated. The also return
>whether the matrix is singular. (Be careful: The solve was successful if the
>result is _false_.) "x" can be a dense vector or a dense matrix.
>
>
>Please have a look at the files and tell me if they work reliably.
>
>mfg
>Gunter
>
>
>
>------------------------------------------------------------------------
>
>/** \file ex_trisolve.cpp \brief main routines. -*- c++ -*- */
>/***************************************************************************
> - begin : 2006-11-03
> - copyright : (C) 2006 by Gunter Winkler
> - email : guwi17_at_[hidden]
> ***************************************************************************/
>
>// Distributed under the Boost Software License, Version 1.0. (See
>// accompanying file LICENSE_1_0.txt or copy at
>// http://www.boost.org/LICENSE_1_0.txt)
>
>// misc. boost includes
>#include <boost/timer.hpp>
>
>#include "trisolve.h"
>
>// boost ublas includes
>#include <boost/numeric/ublas/vector.hpp>
>#include <boost/numeric/ublas/matrix.hpp>
>#include <boost/numeric/ublas/matrix_sparse.hpp>
>#include <boost/numeric/ublas/matrix_proxy.hpp>
>#include <boost/numeric/ublas/triangular.hpp>
>#include <boost/numeric/ublas/io.hpp>
>
>// system includes
>#include <iostream>
>#include <limits>
>
>#ifndef NDEBUG
>const size_t default_size = 10;
>#else
>const size_t default_size = 100;
>#endif
>
>
>template<class M>
>void setup_L(M& A, bool unit=false)
>{
> using namespace boost::numeric::ublas;
> if (unit)
> A = triangular_adaptor<const scalar_matrix<double>,unit_lower>( scalar_matrix<double>(A.size1(), A.size2(), -0.2) );
> else
> A = triangular_adaptor<const scalar_matrix<double>,strict_lower>( scalar_matrix<double>(A.size1(), A.size2(), -0.2) );
>}
>
>template<class M>
>void setup_U(M& A, bool unit=false)
>{
> using namespace boost::numeric::ublas;
> if (unit)
> A = triangular_adaptor<const scalar_matrix<double>,unit_upper>( scalar_matrix<double>(A.size1(), A.size2(), 0.2) );
> else
> A = triangular_adaptor<const scalar_matrix<double>,strict_upper>( scalar_matrix<double>(A.size1(), A.size2(), 0.2) );
>}
>
>template<class V>
>void setup_x(boost::numeric::ublas::vector_expression<V>& x)
>{
> for (size_t i=0; i<x().size(); ++i) {
> x()(i) = 1.0 / sqrt(x().size());
> }
>}
>
>template<class V>
>void setup_x(boost::numeric::ublas::matrix_expression<V>& x)
>{
> double t = 1.0 / sqrt(x().size1());
> for (size_t i=0; i<x().size1(); ++i) {
> row(x(),i) = boost::numeric::ublas::scalar_vector<double>(x().size2(), t );
> }
>}
>
>
>// this is no real norm -> its only useful for this test
>template<class MATRIX>
>double
>norm_2(const boost::numeric::ublas::matrix_expression<MATRIX>& A)
>{
> double result=0.0;
> for (size_t k=0; k<A().size2(); ++k)
> result += norm_2(column( A(), k ));
> return result;
>}
>
>template<class MATRIX>
>double
>norm_inf_2(const boost::numeric::ublas::matrix_expression<MATRIX>& A)
>{
> double result=0.0;
> for (size_t k=0; k<A().size2(); ++k)
> result += norm_inf(column( A(), k ));
> return result;
>}
>
>template<class VECTOR>
>double
>norm_inf_2(const boost::numeric::ublas::vector_expression<VECTOR>& A)
>{
> return norm_inf(A());
>}
>
>template <class TRI, class MATRIX, class VECTOR>
>bool do_test(const boost::numeric::ublas::matrix_expression<MATRIX>& A,
> VECTOR& x, const VECTOR& rhs)
>{
> using namespace boost::numeric::ublas;
>
> double eps = std::numeric_limits<double>::epsilon();
>
> x=rhs;
> boost::timer t;
> t.restart();
> inplace_solve(A(), x, typename TRI::triangular_type());
> double time = t.elapsed();
>
> VECTOR resid = (rhs - prod(A(),x));
> if ( TRI::one(1,1) ) resid -= x;
> double error = norm_2(resid)/norm_2(rhs);
>
>#ifndef NDEBUG
> std::cerr << "||b-Ax|| = " << norm_2(resid) << ", " << norm_inf(resid) << std::endl;
> std::cerr << "||b|| = " << norm_2(rhs) << ", " << norm_inf(rhs) << std::endl;
> std::cerr << "||x|| = " << norm_2(x) << ", " << norm_inf(x) << std::endl;
>#endif
>
> double limit = (A().size1()*eps);
> if (error <= limit) {
> std::cout << " (" << time << " sec)"
> << " relative error " << error
> << std::endl;
> return true;
> } else {
> std::cout << " (" << time << " sec)"
> << " relative error " << error
> << " greater than threshold " << limit
> << std::endl;
>#ifndef NDEBUG
> std::cerr << " rhs = " << rhs << "\n";
> std::cerr << " x = " << x << "\n";
> std::cerr << " res = " << resid << "\n";
>#endif
> return false;
> }
>}
>
>
>template < class MATRIX >
>void run_solver_tests(size_t size, size_t nrhs=5)
>{
> using namespace boost::numeric::ublas;
>
> typedef MATRIX LMAT;
> typedef MATRIX UMAT;
>
> LMAT L(size,size);
> UMAT U(size,size);
>
> L.clear(); U.clear();
>
> vector<double> b(size);
> vector<double> x(size);
> vector<double> x2(size);
>
> matrix<double> X(size,nrhs);
> matrix<double> B(size,nrhs);
>
> setup_L(L,false);
> setup_U(U,false);
>
>
> setup_x(x);
> std::cout << "vec unit lower: ";
> b = x + prod(L, x);
> std::cerr << do_test<unit_lower>(L,x,b);
>
> setup_x(x);
> std::cout << "vec unit upper: ";
> b = x + prod(U, x);
> std::cerr << do_test<unit_upper>(U,x,b);
>
> setup_x(X);
> std::cout << "mat unit lower: ";
> B = X + prod(L, X);
> std::cerr << do_test<unit_lower>(L,X,B);
>
> setup_x(X);
> std::cout << "mat unit upper: ";
> B = X + prod(U, X);
> std::cerr << do_test<unit_upper>(U,X,B);
>
>
> L.clear(); U.clear();
> setup_L(L,true);
> setup_U(U,true);
>
> setup_x(x);
> std::cout << "vec lower: ";
> b = prod(L,x);
> std::cerr << do_test<lower>(L,x,b);
>
> setup_x(x);
> std::cout << "vec upper: ";
> b = prod(U,x);
> std::cerr << do_test<upper>(U,x,b);
>
> setup_x(X);
> std::cout << "mat lower: ";
> B = prod(L, X);
> std::cerr << do_test<lower>(L,X,B);
>
> setup_x(X);
> std::cout << "mat upper: ";
> B = prod(U, X);
> std::cerr << do_test<upper>(U,X,B);
>
>}
>
>
>int main(int argc, char *argv[])
>{
>
> typedef size_t size_type;
> size_type size = default_size;
> if (argc > 1)
> size = ::atoi (argv [1]);
> std::cout << "size " << size << std::endl;
>
> using namespace boost::numeric::ublas;
>
> std::cout << "dense column major\n";
> run_solver_tests< matrix<double, column_major> >(size);
> std::cout << "dense row major\n";
> run_solver_tests< matrix<double, row_major> >(size);
>
> std::cout << "compressed column_major major\n";
> run_solver_tests< compressed_matrix<double, column_major> >(size);
> std::cout << "compressed row major\n";
> run_solver_tests< compressed_matrix<double, row_major> >(size);
>
> std::cout << "mapped column_major major\n";
> run_solver_tests< mapped_matrix<double, column_major> >(size);
> std::cout << "mapped row major\n";
> run_solver_tests< mapped_matrix<double, row_major> >(size);
>
> return EXIT_SUCCESS;
>}
>
>
>------------------------------------------------------------------------
>
>/* trisolve.h - triangular solver */ // -*- c++ -*-
>/***************************************************************************
> begin : 2006-12-11
> copyright : (C) 2006 by Gunter Winkler
> email : guwi17_at_[hidden]
>***************************************************************************/
>
>// Distributed under the Boost Software License, Version 1.0. (See
>// accompanying file LICENSE_1_0.txt or copy at
>// http://www.boost.org/LICENSE_1_0.txt)
>
>#include <boost/numeric/ublas/matrix_sparse.hpp>
>#include <boost/numeric/ublas/matrix_proxy.hpp>
>#include <boost/numeric/ublas/vector_proxy.hpp>
>
>namespace boost { namespace numeric { namespace ublas {
>
> /*********************************
> * PART 0: general documentation *
> *********************************/
>
> /** \brief solve lower triangular system Lx=b or (I+L)x=b.
> * - \param unit controls whether L is strict lower triangular
> (\c true) or lower triangular (\c false)
> - \param MAT is the type of the matrix L
> - \param VEC is the type of the vector x
>
> On input \c x contains the right hand side, on output it
> contains the solution.
>
> The function returns \c true if the matrix is
> singular. Otherwise it returns \c false.
>
> In order to have an efficient implementation, the functions make
> no check whether L is really (strict) lower triangular. Thus the
> result is undefined if this requirement is not met.
>
> You can force the necessary structure by using the \c triangular_adaptor
> with \c lower or \c strict_lower.
>
> \note The solution of (triangular) systems of equations is
> numerically unstable. You should expect large absolute errors if
> the matrix has a large condition number. The main reason is that
> small errors introduced in the first computed components of the
> solution add up to large errors in the final components if a)
> the matrix has many rows/columns and b) the sum of the
> off-diagonal entries of a row is larger than the value at the
> diagonal. You will find a detailed error analysis in any good
> text book about numerics, see, e.g., Golub/van Loan, "Matrix
> Computations."
> */
> template <bool unit, class MAT, class VEC>
> bool inplace_solve_lower(const matrix_expression<MAT>& L, vector_expression<VEC>& x);
>
> /** \brief solve lower triangular system LX=B or (I+L)X=B.
> * - \param unit controls whether L is expected to be strict lower
> triangular (true) or lower triangular (false)
> - \param MAT is the type of the matrix L
> - \param MATX is the type of the rhs-matrix X
>
> On input \c X contains the right hand sides, on output it
> contains the solutions.
>
> \sa inplace_solve_lower(const matrix_expression<MAT>& L, vector_expression<VEC>& x)
> */
> template <bool unit, class MAT, class MATX>
> bool inplace_solve_lower(const matrix_expression<MAT>& L, matrix_expression<MATX>& X);
>
> /** \brief solve upper triangular system Ux=b or (I+U)x=b.
> * - \param unit controls whether U is expected to be strict upper
> triangular (true) or upper triangular (false)
> - \param MAT is the type of the matrix U
> - \param VEC is the type of the vector x
>
> On input \c x contains the right hand side, on output it
> contains the solution.
>
> The function returns \c true if the matrix is
> singular. Otherwise it returns \c false.
>
> In order to have an efficient implementation, the functions make
> no check whether L is really (strict) upper triangular. Thus the
> result is undefined if this requirement is not met.
>
> You can force the necessary structure by using the \c triangular_adaptor
> with \c upper or \c strict_upper.
> */
> template <bool unit, class MAT, class VEC>
> bool inplace_solve_upper(const matrix_expression<MAT>& U, vector_expression<VEC>& x);
>
> /** \brief solve upper triangular system UX=B or (I+U)X=B.
> * - \param unit controls whether U is expected to be strict upper
> triangular (true) or upper triangular (false)
> - \param MAT is the type of the matrix U
> - \param MATX is the type of the rhs-matrix X
>
> On input \c X contains the right hand sides, on output it
> contains the solutions.
>
> \sa inplace_solve_upper(const matrix_expression<MAT>& U, vector_expression<VEC>& x)
> */
> template <bool unit, class MAT, class MATX>
> bool inplace_solve_upper(const matrix_expression<MAT>& U, matrix_expression<MATX>& X);
>
> /*********************************************
> * PART 1: optimized lower and upper solvers *
> *********************************************/
>
> /* if (unit) assume L is strict lower triangular
> * if (!unit) assume L is lower triangular and invertible
> */
> template < bool unit, class T, class Z, class D, size_t IB, class IA, class TA, class VEC >
> bool
> inplace_solve_lower(const compressed_matrix<T, basic_row_major<Z,D>, IB, IA, TA>& L, vector_expression<VEC>& ve)
> {
> if ( (!unit) && L.filled1() <= L.size1() ) return true;
> typedef typename IA::value_type size_type;
> size_type imax = L.filled1()-1;
> for (size_type i = 0; i < imax; ++i ) {
> size_type iptr_a = L.index1_data()[i];
> size_type iptr_b = L.index1_data()[i+1];
> if (!unit) { -- iptr_b; }
> T value = ve()(i);
> for (size_type k = iptr_a; k < iptr_b; ++k) {
> BOOST_UBLAS_CHECK (L.index2_data()[k] <= i, bad_index ());
> value -= L.value_data()[k] * ve()( L.index2_data()[k] );
> }
> if (unit) {
> ve()(i) = value;
> } else {
> ve()(i) = value / L.value_data()[iptr_b];
> }
> }
> return false;
> }
>
> /* if (unit) assume U is strict upper triangular
> * if (!unit) assume U is upper triangular and invertible
> */
> template < bool unit, class T, class Z, class D, size_t IB, class IA, class TA, class VEC >
> bool
> inplace_solve_upper(const compressed_matrix<T, basic_row_major<Z,D>, IB, IA, TA>& U, vector_expression<VEC>& ve)
> {
> if ( (!unit) && U.filled1() <= U.size1() ) return true;
> typedef typename IA::value_type size_type;
> size_type imax = U.filled1()-1;
> for (size_type i = imax; i > 0; --i ) {
> size_type iptr_a = U.index1_data()[i-1];
> size_type iptr_b = U.index1_data()[i];
> if (!unit) { ++ iptr_a; }
> T value = ve()(i-1);
> for (size_type k = iptr_a; k < iptr_b; ++k) {
> BOOST_UBLAS_CHECK (U.index2_data()[k] >= (i-1), bad_index ());
> value -= U.value_data()[k] * ve()( U.index2_data()[k] );
> }
> if (unit) {
> ve()(i-1) = value;
> } else {
> ve()(i-1) = value / U.value_data()[iptr_a-1];
> }
> }
> return false;
> }
>
> /*******************************************
> * PART 2: general lower and upper solvers *
> *******************************************/
>
> /* only the lower left or upper right triangle of A is accessed,
> * the number of rows of A determines the (used) size of x,
> *
> * syntax 1:
> * inplace_solve_lower<unit>(A,x)
> * inplace_solve_upper<unit>(A,x)
> *
> * syntax 2:
> * inplace_solve_lower<unit>(A,x,ORI)
> * inplace_solve_upper<unit>(A,x,ORI)
> * where ORI is row_major_tag or column_major_tag
> * (here you can override the default choice of the orientation of the algorithm)
> */
>
> /** \brief Solve (column-by-column) Lx=b or (I+L)x=b where L is the (strict) lower triangle of A.
>
> This version uses the generic algorithm for dense matrices and
> thus only accesses the (strict) lower triangle of A. The result
> is even correct when the matrix A has non-zero value (at or)
> above the diagonal.
> */
> template <bool unit, class MAT, class VEC>
> bool inplace_solve_lower(const matrix_expression<MAT>& A, vector_expression<VEC>& x, column_major_tag)
> {
> typedef typename VEC::size_type size_type;
> typedef typename VEC::value_type value_type;
>
> BOOST_UBLAS_CHECK (A ().size2 () >= x ().size (), bad_size ());
> BOOST_UBLAS_CHECK (A ().size1 () <= x ().size (), bad_size ());
>
> size_type n = A().size1();
> if (n==0) return false;
>
> for (size_type j = 0; j < (n-1); ++j) {
> double t;
> if (!unit) {
> if (A()(j,j)==0.0) return true;
> t = ( x()(j) /= A()(j,j) );
> } else {
> t = x()(j);
> }
> if (t != value_type()) {
> noalias( project(x(), range(j+1,n)) ) -= t*project( column(A(), j), range(j+1,n) );
> }
> }
> if (!unit) {
> if (A()(n-1,n-1)==0.0) return true;
> x()(n-1) /= A()(n-1,n-1);
> }
> return false;
> }
>
> /** \brief Solve (row-by-row) Lx=b or (I+L)x=b where L is the (strict) lower triangle of A.
>
> This version uses the generic algorithm for dense matrices and
> thus only accesses the (strict) lower triangle of A. The result
> is even correct when the matrix A has non-zero value (at or)
> above the diagonal.
> */
> template <bool unit, class MAT, class VEC>
> bool inplace_solve_lower(const matrix_expression<MAT>& A, vector_expression<VEC>& x, row_major_tag)
> {
> typedef typename VEC::size_type size_type;
> typedef typename VEC::value_type value_type;
>
> BOOST_UBLAS_CHECK (A ().size2 () >= x ().size (), bad_size ());
> BOOST_UBLAS_CHECK (A ().size1 () <= x ().size (), bad_size ());
>
> size_type n = A().size1();
> if (n==0) return false;
> if (!unit) {
> if (A()(0,0)==0.0) return true;
> x()(0) /= A()(0,0);
> }
>
> for (size_type j = 1; j < n; ++j) {
> double t;
> t = inner_prod( project(x(), range(0,j)), project( row(A(), j), range(0,j) ) );
> if (!unit) {
> if (A()(j,j)==0.0) return true;
> x()(j) = (x()(j) - t)/A()(j,j);
> } else {
> x()(j) -= t;
> }
> }
> return false;
> }
>
> /** \brief Solve (column-by-column) Ux=b or (I+U)x=b where U is the (strict) upper triangle of A.
>
> This version uses the generic algorithm for dense matrices and
> thus only accesses the (strict) upper triangle of A. The result
> is even correct when the matrix A has non-zero value (at or)
> below the diagonal.
> */
> template <bool unit, class MAT, class VEC>
> bool inplace_solve_upper(const matrix_expression<MAT>& A, vector_expression<VEC>& x, column_major_tag)
> {
> typedef typename VEC::size_type size_type;
> typedef typename VEC::value_type value_type;
>
> BOOST_UBLAS_CHECK (A ().size2 () >= x ().size (), bad_size ());
> BOOST_UBLAS_CHECK (A ().size1 () <= x ().size (), bad_size ());
>
> size_type n = A().size1();
> if (n==0) return false;
>
> for (size_type j = n-1; j > 0; --j) {
> double t;
> if (!unit) {
> if (A()(j,j)==0.0) return true;
> t = ( x()(j) /= A()(j,j) );
> } else {
> t = x()(j);
> }
> if (t != value_type()) {
> noalias( project(x(), range(0,j)) ) -= t*project( column(A(), j), range(0,j) );
> }
> }
> if (!unit) {
> if (A()(0,0)==0.0) return true;
> x()(0) /= A()(0,0);
> }
> return false;
> }
>
> /** \brief Solve (row-by-row) Ux=b or (I+U)x=b where U is the (strict) upper triangle of A.
>
> This version uses the generic algorithm for dense matrices and
> thus only accesses the (strict) upper triangle of A. The result
> is even correct when the matrix A has non-zero value (at or)
> below the diagonal.
> */
> template <bool unit, class MAT, class VEC>
> bool inplace_solve_upper(const matrix_expression<MAT>& A, vector_expression<VEC>& x, row_major_tag)
> {
> typedef typename VEC::size_type size_type;
> typedef typename VEC::value_type value_type;
>
> BOOST_UBLAS_CHECK (A ().size2 () >= x ().size (), bad_size ());
> BOOST_UBLAS_CHECK (A ().size1 () <= x ().size (), bad_size ());
>
> size_type n = A().size1();
> if (n==0) return false;
> if (!unit) {
> x()(n-1) /= A()(n-1,n-1);
> }
>
> for (size_type j = n-2; /* at end */ ; --j) {
> double t;
> t = inner_prod( project(x(), range(j+1,n)), project( row(A(), j), range(j+1,n) ) );
> if (!unit) {
> if (A()(j,j)==0.0) return true;
> x()(j) = (x()(j) - t)/A()(j,j);
> } else {
> x()(j) -= t;
> }
> if (j==0) break;
> }
> return false;
> }
>
> /* only the lower left or upper right triangle of A is accessed,
> * the number of rows of A determines the (used) size of X,
> *
> * syntaX 1:
> * inplace_solve_lower<unit>(A,X)
> * inplace_solve_upper<unit>(A,X)
> *
> * syntaX 2:
> * inplace_solve_lower<unit>(A,X,ORI)
> * inplace_solve_upper<unit>(A,X,ORI)
> * where ORI is row_major_tag or column_major_tag
> * (here you can override the default choice of the orientation of the algorithm)
> */
>
> /** \brief Solve (column-by-column) LX=b or (I+L)X=b where L is the (strict) lower triangle of A.
>
> This version uses the generic algorithm for dense matrices and
> thus only accesses the (strict) lower triangle of A. The result
> is even correct when the matriX A has non-zero value (at or)
> above the diagonal.
> */
> template <bool unit, class MAT, class VEC>
> bool inplace_solve_lower(const matrix_expression<MAT>& A, matrix_expression<VEC>& X, column_major_tag)
> {
> typedef typename VEC::size_type size_type;
> typedef typename VEC::value_type value_type;
>
> BOOST_UBLAS_CHECK (A ().size2 () >= X ().size1 (), bad_size ());
> BOOST_UBLAS_CHECK (A ().size1 () <= X ().size1 (), bad_size ());
>
> size_type nrhs = X().size2();
>
> size_type n = A().size1();
> if (n==0) return false;
>
> for (size_type j = 0; j < (n-1); ++j) {
> if (!unit) {
> if (A()(j,j)==0.0) return true;
> row(X(),j) /= A()(j,j) ;
> }
> noalias( project(X(), range(j+1,n), range(0, nrhs) ) )
> -= outer_prod(project( column(A(), j), range(j+1,n) ), row(X(),j));
> }
> if (!unit) {
> if (A()(n-1,n-1)==0.0) return true;
> row(X(),n-1) /= A()(n-1,n-1);
> }
> return false;
> }
>
> /** \brief Solve (row-by-row) LX=b or (I+L)X=b where L is the (strict) lower triangle of A.
>
> This version uses the generic algorithm for dense matrices and
> thus only accesses the (strict) lower triangle of A. The result
> is even correct when the matriX A has non-zero value (at or)
> above the diagonal.
> */
> template <bool unit, class MAT, class VEC>
> bool inplace_solve_lower(const matrix_expression<MAT>& A, matrix_expression<VEC>& X, row_major_tag)
> {
> typedef typename VEC::size_type size_type;
> typedef typename VEC::value_type value_type;
>
> BOOST_UBLAS_CHECK (A ().size2 () >= X ().size1 (), bad_size ());
> BOOST_UBLAS_CHECK (A ().size1 () <= X ().size1 (), bad_size ());
>
> size_type nrhs = X().size2();
>
> size_type n = A().size1();
> if (n==0) return false;
> if (!unit) {
> if (A()(0,0)==0.0) return true;
> row(X(),0) /= A()(0,0);
> }
>
> for (size_type j = 1; j < n; ++j) {
> if (!unit) {
> if (A()(j,j)==0.0) return true;
> noalias( row(X(),j) )
> = ( row(X(),j) - prod( project( row(A(), j), range(0,j) ),
> project(X(), range(0,j), range(0,nrhs)) ) )
> / A()(j,j);
> } else {
> noalias( row(X(),j) )
> -= prod( project( row(A(), j), range(0,j) ),
> project(X(), range(0,j), range(0,nrhs)) );
> }
> }
> return false;
> }
>
> /** \brief Solve (column-by-column) UX=b or (I+U)X=b where U is the (strict) upper triangle of A.
>
> This version uses the generic algorithm for dense matrices and
> thus only accesses the (strict) upper triangle of A. The result
> is even correct when the matriX A has non-zero value (at or)
> below the diagonal.
> */
> template <bool unit, class MAT, class VEC>
> bool inplace_solve_upper(const matrix_expression<MAT>& A, matrix_expression<VEC>& X, column_major_tag)
> {
> typedef typename VEC::size_type size_type;
> typedef typename VEC::value_type value_type;
>
> BOOST_UBLAS_CHECK (A ().size2 () >= X ().size1 (), bad_size ());
> BOOST_UBLAS_CHECK (A ().size1 () <= X ().size1 (), bad_size ());
>
> size_type nrhs = X().size2();
>
> size_type n = A().size1();
> if (n==0) return false;
>
> for (size_type j = n-1; j > 0; --j) {
> if (!unit) {
> if (A()(j,j)==0.0) return true;
> row(X(),j) /= A()(j,j);
> }
> noalias( project(X(), range(0,j), range(0,nrhs) ) )
> -= outer_prod( project( column(A(), j), range(0,j) ),
> row(X(),j) );
> }
> if (!unit) {
> if (A()(0,0)==0.0) return true;
> row(X(),0) /= A()(0,0);
> }
> return false;
> }
>
> /** \brief Solve (row-by-row) UX=b or (I+U)X=b where U is the (strict) upper triangle of A.
>
> This version uses the generic algorithm for dense matrices and
> thus only accesses the (strict) upper triangle of A. The result
> is even correct when the matriX A has non-zero value (at or)
> below the diagonal.
> */
> template <bool unit, class MAT, class VEC>
> bool inplace_solve_upper(const matrix_expression<MAT>& A, matrix_expression<VEC>& X, row_major_tag)
> {
> typedef typename VEC::size_type size_type;
> typedef typename VEC::value_type value_type;
>
> BOOST_UBLAS_CHECK (A ().size2 () >= X ().size1 (), bad_size ());
> BOOST_UBLAS_CHECK (A ().size1 () <= X ().size1 (), bad_size ());
>
> size_type nrhs = X().size2();
>
> size_type n = A().size1();
> if (n==0) return false;
> if (!unit) {
> row(X(),n-1) /= A()(n-1,n-1);
> }
>
> for (size_type j = n-2; /* at end */ ; --j) {
> if (!unit) {
> if (A()(j,j)==0.0) return true;
> noalias( row(X(),j) )
> = ( row(X(),j) - prod( project( row(A(), j), range(j+1,n) ),
> project(X(), range(j+1,n), range(0,nrhs)) ) )
> / A()(j,j);
> } else {
> noalias( row(X(),j) )
> -= prod( project( row(A(), j), range(j+1,n) ),
> project(X(), range(j+1,n), range(0,nrhs)) );
> }
> if (j==0) break;
> }
> return false;
> }
>
> // Dispatcher for default orientation
>
> template <bool unit, class MAT, class VEC>
> bool inplace_solve_lower(const matrix_expression<MAT>& A, vector_expression<VEC>& x)
> {
> return inplace_solve_lower<unit>(A, x, typename MAT::orientation_category() );
> }
>
> template <bool unit, class MAT, class VEC>
> bool inplace_solve_upper(const matrix_expression<MAT>& A, vector_expression<VEC>& x)
> {
> return inplace_solve_upper<unit>(A, x, typename MAT::orientation_category() );
> }
>
> template <bool unit, class MAT, class MATX>
> bool inplace_solve_lower(const matrix_expression<MAT>& A, matrix_expression<MATX>& X)
> {
> return inplace_solve_lower<unit>(A, X, typename MAT::orientation_category() );
> }
>
> template <bool unit, class MAT, class MATX>
> bool inplace_solve_upper(const matrix_expression<MAT>& A, matrix_expression<MATX>& X)
> {
> return inplace_solve_upper<unit>(A, X, typename MAT::orientation_category() );
> }
>
> /***********************************
> * PART 3: upper/lower dispatchers *
> ***********************************/
>
> // note: the dispatcher even work when VEC is a matrix expression.
>
> template < class MAT, class VEC >
> bool
> inplace_solve( const MAT& L, VEC& x, unit_lower_tag )
> {
> return inplace_solve_lower<true>(L, x);
> }
>
> template < class MAT, class VEC >
> bool
> inplace_solve( const MAT& L, VEC& x, lower_tag )
> {
> return inplace_solve_lower<false>(L, x);
> }
>
> template < class MAT, class VEC >
> bool
> inplace_solve( const MAT& U, VEC& x, unit_upper_tag )
> {
> return inplace_solve_upper<true>(U, x);
> }
>
> template < class MAT, class VEC >
> bool
> inplace_solve( const MAT& U, VEC& x, upper_tag )
> {
> return inplace_solve_upper<false>(U, x);
> }
>
>}}}
>
>
>------------------------------------------------------------------------
>
>Index: functional.hpp
>===================================================================
>RCS file: /cvsroot/boost/boost/boost/numeric/ublas/functional.hpp,v
>retrieving revision 1.39
>diff -u -p -r1.39 functional.hpp
>--- functional.hpp 27 Jul 2006 10:27:37 -0000 1.39
>+++ functional.hpp 15 Dec 2006 13:30:57 -0000
>@@ -1774,6 +1786,7 @@ namespace boost { namespace numeric { na
> template <class Z>
> struct basic_lower {
> typedef Z size_type;
>+ typedef lower_tag triangular_type;
>
> template<class L>
> static
>@@ -1828,6 +1841,7 @@ namespace boost { namespace numeric { na
> template <class Z>
> struct basic_upper {
> typedef Z size_type;
>+ typedef upper_tag triangular_type;
>
> template<class L>
> static
>@@ -1882,6 +1896,7 @@ namespace boost { namespace numeric { na
> template <class Z>
> struct basic_unit_lower : public basic_lower<Z> {
> typedef Z size_type;
>+ typedef unit_lower_tag triangular_type;
>
> template<class L>
> static
>@@ -1925,6 +1940,7 @@ namespace boost { namespace numeric { na
> template <class Z>
> struct basic_unit_upper : public basic_upper<Z> {
> typedef Z size_type;
>+ typedef unit_upper_tag triangular_type;
>
> template<class L>
> static
>@@ -1968,6 +1984,7 @@ namespace boost { namespace numeric { na
> template <class Z>
> struct basic_strict_lower : public basic_lower<Z> {
> typedef Z size_type;
>+ typedef strict_lower_tag triangular_type;
>
> template<class L>
> static
>@@ -2026,6 +2043,7 @@ namespace boost { namespace numeric { na
> template <class Z>
> struct basic_strict_upper : public basic_upper<Z> {
> typedef Z size_type;
>+ typedef strict_upper_tag triangular_type;
>
> template<class L>
> static
>
>
>------------------------------------------------------------------------
>
>_______________________________________________
>ublas mailing list
>ublas_at_[hidden]
>http://listarchives.boost.org/mailman/listinfo.cgi/ublas
>
>
Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm