|
|
Ublas : |
From: Gunter Winkler (gunter.winkler_at_[hidden])
Date: 2006-10-09 02:41:07
On Saturday 07 October 2006 20:53, Michael Stevens wrote:
> Hi Gunter, Hi Guimond,
>
> On Thursday, 7. September 2006 17:32, Gunter Winkler wrote:
> > On Thursday 07 September 2006 16:28, Guimond, Alexandre wrote:
> > Interesting, the functions can not be distinguished. I always use the
> > version with pivoting, which works perfectly.
>
> GCC agrees the two function definitions are ambigeous. The two functions in
> question are.
>
> line 258:
> template<class M, class E>
> void lu_substitute (const M &m, vector_expression<E> &e) {
> inplace_solve (m, e, unit_lower_tag ());
> inplace_solve (m, e, upper_tag ());
>
> line 317:
> template<class E, class M>
> void lu_substitute (matrix_expression<E> &e, const M &m) {
> inplace_solve (e, m, upper_tag ());
> inplace_solve (e, m, unit_lower_tag ());
>
> > Maybe one should drop the (IMHO strange) versions of lu substitute for
> > solving X A = B .
>
> The function (defined at line 317 ) certainly looks wrong. The argument 'm'
> is const. This 'm' is passed to inplace_solver as the second argument this
> is not const and is modified to contain the result.
No, its correct. m is the factorized matrix, e contains the right hand sides.
>
> Unless anyone knows what the function is there for and how to code it
> correctly I will remove it.
The purpose is the transposed solve. Similar to the two types of prod
( prod(M,v) ... M*v and prod(v,M) ... M'*v ) there are two types of triangular
solves: solve(M,v,...) ... M*x=v and solve(v,M,...) ... M'*x=v.
The above function tries to implement the same idea for multiple right hand
sides v which are stored in the matrix_expression e. I think we can drop the
transposed solves (or rename them to ..._trans).
btw. just calling lu_substitue( trans(M), x ) does not work because the
diagonal elements only belong to the upper triangle of M.
mfg
Gunter