Re: Systems of Linear Ecuations

"oswin krause" <[email protected]>
Newsgroups gmane.comp.lib.boost.ublas
Message-ID <[email protected]> 
Hi,

I think my last mail was a bit off topic as an answer to the warnign. Sorry!

This is already implemented in triangular.hpp as inplace_solve. However 
the implementation is not complete.

This reminds me that I wanted to send the improved implementation of my 
fork back to ublas as I think that it is quite complete. of course there 
is no SSE or stuff like that so there is a lot of space to improve. it 
is also not a blocked implementation but given the state of the 
matrix-matrix multiply this is hardly an efficiency issue. The 
implementation however crucially relies on the efficiency of the 
matrix-proxies, but if they are slow (in the sparse/packed case) the 
proxies should be improved. I have not benchmarked that.

the implementation is attached, namespaces need to be changed. otherwise 
everything is there.

Greetings,
Oswin

On 02.05.2013 00:01, Iulian Calciu wrote:
> Hello,
>
> I have nowhere seen a proposal for Systems of Linear Ecuations solving.
> I know they can be solved with LU descomposition for quadratic systems.
>
> Eg:
> Ax = b;
> A = LU with LU transform;
> LUx = b; note Ux = y, then
> Ly = b => y = LTRIS(L, b), then
> Ux = y => x = UTRIS(U, y);
>
> Note: LTRIS and UTRIS are methods for solving triangular lower/upper 
> systems of ecuations.
>
> LTRIS
> x = b;
> for i = 1 : n
> for j = 1 : i - 1
> x(i) = x(i) - L(i,j)x(j)
> end
> x(i) = x(i) / L(i,i);
> end
>
> UTRIS:
> x = b;
> for i = n : -1: 1
> for j = i + 1 : n
> x(i) = x(i) - U(i,j)x(j)
> end
> x(i) = x(i) / U(i,i);
> end
>
> Can this thing be proposed in Boost uBLAS?
>
> Thank you,
> Iulian Calciu
>
>
>
>
>
> _______________________________________________
> ublas mailing list
> [email protected]
> http://lists.boost.org/mailman/listinfo.cgi/ublas
> Sent to: [email protected]

_______________________________________________
ublas mailing list
[email protected]
http://lists.boost.org/mailman/listinfo.cgi/ublas
Sent to: [email protected]
triangular.hpp (text/x-c++hdr, 8.8 KB) 
#ifndef SHARK_BLAS_TRIANGULAR_HPP
#define SHARK_BLAS_TRIANGULAR_HPP

#include <shark/LinAlg/BLAS/ublas/matrix.hpp>
#include <shark/LinAlg/BLAS/ublas/matrix_proxy.hpp>

namespace shark {
namespace blas {

//Lower triangular(row-major) - vector
template<bool Unit, class E1, class E2>
void inplace_solve(
	matrix_expression<E1> const& e1, vector_expression<E2> &e2,
        lower_tag, column_major_tag
) {
	SIZE_CHECK(e1().size1() == e1().size2());
	SIZE_CHECK(e1().size2() == e2().size());
	
	typedef typename E1::value_type value_type;
	
	std::size_t size = e2().size();
	for (std::size_t n = 0; n < size; ++ n) {
		if(!Unit){
			RANGE_CHECK(e1()(n, n) != value_type());//matrix is singular
			e2()(n) /= e1()(n, n);
		}
		if (e2()(n) != value_type/*zero*/()){
			matrix_column<E1 const> col = column(e1(),n);
			noalias(subrange(e2(),n+1,size)) -= e2()(n) * subrange(col,n+1,size);
		}
	}
}
//Lower triangular(column-major) - vector
template<bool Unit, class E1, class E2>
void inplace_solve(
	matrix_expression<E1> const& e1, vector_expression<E2> &e2,
        lower_tag, row_major_tag
) {
	SIZE_CHECK(e1().size1() == e1().size2());
	SIZE_CHECK(e1().size2() == e2().size());
	
	typedef typename E1::value_type value_type;
	
	std::size_t size = e2().size();
	for (std::size_t n = 0; n < size; ++ n) {
		matrix_row<E1 const> matRow = row(e1(),n);
		e2()(n) -= inner_prod(subrange(matRow,0,n),subrange(e2(),0,n));

		if(!Unit){
			RANGE_CHECK(e1()(n, n) != value_type());//matrix is singular
			e2()(n) /= e1()(n, n);
		}
	}
}

//upper triangular(column-major)-vector
template<bool Unit, class E1, class E2>
void inplace_solve(
	matrix_expression<E1> const& e1, vector_expression<E2> &e2,
        upper_tag, column_major_tag
) {
	SIZE_CHECK(e1().size1() == e1().size2());
	SIZE_CHECK(e1().size2() == e2().size());
	
	typedef typename E1::value_type value_type;
	
	std::size_t size = e2().size();
	for (std::size_t i = 0; i < size; ++ i) {
		std::size_t n = size-i-1;
		if(!Unit){
			RANGE_CHECK(e1()(n, n) != value_type());//matrix is singular
			e2()(n) /= e1()(n, n);
		}
		if (e2()(n) != value_type/*zero*/()) {
			matrix_column<E1 const> col = column(e1(),n);
			noalias(subrange(e2(),0,n)) -= e2()(n) * subrange(col,0,n);
		}
	}
}
//upper triangular(row-major)-vector
template<bool Unit, class E1, class E2>
void inplace_solve(
	matrix_expression<E1> const& e1, vector_expression<E2> &e2,
        upper_tag, row_major_tag
) {
	SIZE_CHECK(e1().size1() == e1().size2());
	SIZE_CHECK(e1().size2() == e2().size());
	
	typedef typename E1::value_type value_type;
	
	std::size_t size = e1().size1();
	for (std::size_t i = 0; i < size; ++ i) {
		std::size_t n = size-i-1;
		matrix_row<E1 const> matRow = row(e1(),n);
		e2()(n) -= inner_prod(subrange(matRow,n+1,size),subrange(e2(),n+1,size));
		if(!Unit){
			RANGE_CHECK(e1()(n, n) != value_type());//matrix is singular
			e2()(n) /= e1()(n, n);
		}
	}
}

//public interface

// Dispatcher for Systems of the form Ax=b
template<class E1, class E2>
void inplace_solve(matrix_expression<E1> const& e1, vector_expression<E2> &e2, lower_tag tag){
	typedef typename E1::orientation_category orientation_category;
	inplace_solve<false>(e1, e2, lower_tag(), orientation_category());
}
template<class E1, class E2>
void inplace_solve(matrix_expression<E1> const& e1, vector_expression<E2> &e2, unit_lower_tag){
	typedef typename E1::orientation_category orientation_category;
	inplace_solve<true>(e1, e2, lower_tag(), orientation_category());
}

template<class E1, class E2>
void inplace_solve(matrix_expression<E1> const& e1, vector_expression<E2> &e2, upper_tag){
	typedef typename E1::orientation_category orientation_category;
	inplace_solve<false>(e1, e2,upper_tag(), orientation_category());
}
template<class E1, class E2>
void inplace_solve(matrix_expression<E1> const& e1, vector_expression<E2> &e2, unit_upper_tag){
	typedef typename E1::orientation_category orientation_category;
	inplace_solve<true>(e1, e2,upper_tag(), orientation_category());
}

// Dispatcher for Systems of the form x^TA=b
template<class E1, class E2>
void inplace_solve(vector_expression<E1> &e1, const matrix_expression<E2> &e2, lower_tag) {
	inplace_solve(trans(e2), e1, upper_tag());
}
template<class E1, class E2>
void inplace_solve(vector_expression<E1> &e1, const matrix_expression<E2> &e2, upper_tag) {
	inplace_solve(trans(e2), e1, lower_tag());
}
template<class E1, class E2>
void inplace_solve(vector_expression<E1> &e1, const matrix_expression<E2> &e2, unit_lower_tag) {
	inplace_solve(trans(e2), e1, unit_upper_tag());
}
template<class E1, class E2>
void inplace_solve(vector_expression<E1> &e1, const matrix_expression<E2> &e2, unit_upper_tag) {
	inplace_solve(trans(e2), e1, unit_lower_tag());
}


// Operations:
//  k * n * (n - 1) / 2 + k * n = k * n * (n + 1) / 2 multiplications,
//  k * n * (n - 1) / 2 additions

// Lower triangular(column major) - matrix
template<bool Unit, class E1, class E2>
void inplace_solve(
	matrix_expression<E1> const& e1, matrix_expression<E2> &e2, 
	lower_tag, column_major_tag
) {
	SIZE_CHECK(e1().size1() == e1().size2());
	SIZE_CHECK(e1().size2() == e2().size1());
	
	typedef typename E1::value_type value_type;
	
	std::size_t size1 = e2().size1();
	std::size_t size2 = e2().size2();
	for (std::size_t n = 0; n < size1; ++ n) {
		matrix_column<E1 const> columnTriangular = column(e1(),n);
		for (std::size_t l = 0; l < size2; ++ l) {
			if(!Unit){
				RANGE_CHECK(e1()(n, n) != value_type());//matrix is singular
				e2()(n, l) /= e1()(n, n);
			}
			if (e2()(n, l) != value_type/*zero*/()) {
				matrix_column<E2> columnMatrix = column(e2(),l);
				noalias(subrange(columnMatrix,n+1,size1)) -= e2()(n,l) * subrange(columnTriangular,n+1,size1);
			}
		}
	}
}
// Lower triangular(row major) - matrix
template<bool Unit, class E1, class E2>
void inplace_solve(
	matrix_expression<E1> const& e1, matrix_expression<E2> &e2, 
	lower_tag, row_major_tag
) {
	SIZE_CHECK(e1().size1() == e1().size2());
	SIZE_CHECK(e1().size2() == e2().size1());
	
	typedef typename E1::value_type value_type;
	
	std::size_t size1 = e2().size1();
	for (std::size_t n = 0; n < size1; ++ n) {
		for (std::size_t m = 0; m < n; ++m) {
			noalias(row(e2(),n)) -= e1()(n,m)*row(e2(),m);
		}
		if(!Unit){
			RANGE_CHECK(e1()(n, n) != value_type());//matrix is singular
			row(e2(),n)/=e1()(n, n);
		}
	}
}

//Upper triangular(column major) - matrix
template<bool Unit, class E1, class E2>
void inplace_solve(
	matrix_expression<E1> const& e1, matrix_expression<E2> &e2,
        upper_tag, column_major_tag
) {
	SIZE_CHECK(e1().size1() == e1().size2());
	SIZE_CHECK(e1().size2() == e2().size1());
	
	typedef typename E1::value_type value_type;
	
	std::size_t size1 = e2().size1();
	std::size_t size2 = e2().size2();
	for (std::size_t i = 0; i < size1; ++ i) {
		std::size_t n = size1-i-1;
		matrix_column<E1 const> columnTriangular = column(e1(),n);
		if(!Unit){
			RANGE_CHECK(e1()(n, n) != value_type());//matrix is singular
			row(e2(),n) /= e1()(n, n);
		}
		for (std::size_t l = 0; l < size2; ++ l) {
			if (e2()(n, l) != value_type/*zero*/()) {
				matrix_column<E2> columnMatrix = column(e2(),l);
				noalias(subrange(columnMatrix,0,n)) -= e2()(n,l) * subrange(columnTriangular,0,n);
			}
		}
	}
}

//Upper triangular(row major) - matrix
template<bool Unit, class E1, class E2>
void inplace_solve(
	matrix_expression<E1> const& e1, matrix_expression<E2> &e2,
        upper_tag, row_major_tag
) {
	SIZE_CHECK(e1().size1() == e1().size2());
	SIZE_CHECK(e1().size2() == e2().size1());
	
	typedef typename E1::value_type value_type;
	
	std::size_t size1 = e2().size1();
	for (std::size_t i = 0; i < size1; ++ i) {
		std::size_t n = size1-i-1;
		for (std::size_t m = n+1; m < size1; ++m) {
			noalias(row(e2(),n)) -= e1()(n,m)*row(e2(),m);
		}
		if(!Unit){
			RANGE_CHECK(e1()(n, n) != value_type());//matrix is singular
			row(e2(),n)/=e1()(n, n);
		}
	}
}


// Dispatcher
template<class E1, class E2>
void inplace_solve(matrix_expression<E1> const& e1, matrix_expression<E2> &e2, lower_tag) {
	typedef typename E1::orientation_category orientation_category;
	inplace_solve<false>(e1, e2, lower_tag(), orientation_category());
}
template<class E1, class E2>
void inplace_solve(matrix_expression<E1> const& e1, matrix_expression<E2> &e2, upper_tag) {
	typedef typename E1::orientation_category orientation_category;
	inplace_solve<false>(e1, e2, upper_tag(), orientation_category());
}
template<class E1, class E2>
void inplace_solve(matrix_expression<E1> const& e1, matrix_expression<E2> &e2, unit_lower_tag) {
	typedef typename E1::orientation_category orientation_category;
	inplace_solve<true>(e1, e2, lower_tag(), orientation_category());
}
template<class E1, class E2>
void inplace_solve(matrix_expression<E1> const& e1, matrix_expression<E2> &e2, unit_upper_tag) {
	typedef typename E1::orientation_category orientation_category;
	inplace_solve<true>(e1, e2, upper_tag(), orientation_category());
}


}
}

#endif
lmpx.com only provides a reader for public news (NNTP) servers. It is not affiliated with the servers or forums shown here and is not responsible for the content of articles, which is written by their respective authors.