Help with lu_factorize

[Marco Guazzone]

Hello!

I don't understand how the function lu_factorize(A) works (note: only
the version with one parameter, which I think is for LU without
pivoting).
I didn't find any documentation of this function so I've decided to
experiment a bit with it.

I've attached an example.
The input matrix is:
---[snip]---
Input matrix
        matrix:
[4,4]((0.55595,0.27469,0.540605,0.798938),(0.108929,0.830123,0.891726,0.895283),(0.948014,0.973234,0.216504,0.883152),(0.023787,0.675382,0.231751,0.450332))
---[/snip]---

In this example I use both lu_factorize(A) and lu_factorize(A,P)
(i.e.,with the permutation matrix).
The last one works perfectly. After call it, A is what I expected
(also verified with Mathematica and Octave).
---[snip]---
LU *with* permutation matrix
        res: 0
        matrix:
[4,4]((0.948014,0.973234,0.216504,0.883152),(0.114902,0.718296,0.866849,0.793807),(0.586436,-0.412156,0.770916,0.608198),(0.0250914,0.906259,-0.725463,0.150003))
---[/snip]---

The first one instead gives a result that I'm not able to understand.
---[snip]---
LU *without* permutation matrix
        res: 0
        matrix:
[4,4]((0.55595,0.27469,0.540605,0.798938),(0.195933,0.776302,0.785804,0.738745),(1.70521,0.650299,-1.21635,-0.959614),(0.0427862,0.854859,0.380754,0.150003))
---[/snip]---

NOTE: I'm using Boost 1.40 and GCC 4.4.1. I've compiled with:
$ g++ -o lu -Wall -Wextra -pedantic -ansi lu.cpp

Can you give me an help, please?

Best,

-- Marco

[lu.cpp]

---

#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <iostream>

int main()
{
        typedef double value_type;
    typedef boost::numeric::ublas::matrix<value_type> matrix_type;
    typedef boost::numeric::ublas::permutation_matrix<value_type> pmatrix_type;

    matrix_type A(4,4);

    A(0,0) = 0.555950; A(0,1) = 0.274690; A(0,2) = 0.540605; A(0,3) = 0.798938;
    A(1,0) = 0.108929; A(1,1) = 0.830123; A(1,2) = 0.891726; A(1,3) = 0.895283;
    A(2,0) = 0.948014; A(2,1) = 0.973234; A(2,2) = 0.216504; A(2,3) = 0.883152;
    A(3,0) = 0.023787; A(3,1) = 0.675382; A(3,2) = 0.231751; A(3,3) = 0.450332;

        int res;

        // LU without pivoting
        matrix_type A1(A);

        res = boost::numeric::ublas::lu_factorize(A1);

        // LU with partial pivoting
        matrix_type A2(A);

        pmatrix_type P(A2.size1());

        res = boost::numeric::ublas::lu_factorize(A2, P);

        std::cout << "Input matrix" << std::endl;
        std::cout << "\tmatrix: " << A << std::endl;
        std::cout << "LU *without* permutation matrix" << std::endl;
        std::cout << "\tres: " << res << std::endl;
        std::cout << "\tmatrix: " << A1 << std::endl;
        std::cout << "LU *with* permutation matrix" << std::endl;
        std::cout << "\tres: " << res << std::endl;
        std::cout << "\tmatrix: " << A2 << std::endl;
}

---

_______________________________________________
ublas mailing list
ublas@...
http://lists.boost.org/mailman/listinfo.cgi/ublas
Sent to: lists@...

[Nasos Iliopoulos]

Some parts of this message have been removed. Learn more about Nabble's security policy.

Marco,
lu_factorize(A) : A == L*U
lu_factorize(A2, P) : P*A2 ==L*U
P is a permutation matrix that interchanges the rows of A2. In some cases a matrix cannot be factored, so using pivoting with the permutation matrix can give you a factorisation. Also solving a system with and without pivoting has some slight differences.

http://math.fullerton.edu/mathews/n2003/LUFactorMod.html

I hope this helps.

Best
Nasos

> Date: Tue, 15 Sep 2009 17:46:15 +0200
> From: marco.guazzone@...
> To: ublas@...
> Subject: [ublas] Help with lu_factorize
>
> Hello!
>
> I don't understand how the function lu_factorize(A) works (note: only
> the version with one parameter, which I think is for LU without
> pivoting).
> I didn't find any documentation of this function so I've decided to
> experiment a bit with it.
>
> I've attached an example.
> The input matrix is:
> ---[snip]---
> Input matrix
> matrix:
> [4,4]((0.55595,0.27469,0.540605,0.798938),(0.108929,0.830123,0.891726,0.895283),(0.948014,0.973234,0.216504,0.883152),(0.023787,0.675382,0.231751,0.450332))
> ---[/snip]---
>
> In this example I use both lu_factorize(A) and lu_factorize(A,P)
> (i.e.,with the permutation matrix).
> The last one works perfectly. After call it, A is what I expected
> (also verified with Mathematica and Octave).
> ---[snip]---
> LU *with* permutation matrix
> res: 0
> matrix:
> [4,4]((0.948014,0.973234,0.216504,0.883152),(0.114902,0.718296,0.866849,0.793807),(0.586436,-0.412156,0.770916,0.608198),(0.0250914,0.906259,-0.725463,0.150003))
> ---[/snip]---
>
> The first one instead gives a result that I'm not able to understand.
> ---[snip]---
> LU *without* permutation matrix
> res: 0
> matrix:
> [4,4]((0.55595,0.27469,0.540605,0.798938),(0.195933,0.776302,0.785804,0.738745),(1.70521,0.650299,-1.21635,-0.959614),(0.0427862,0.854859,0.380754,0.150003))
> ---[/snip]---
>
> NOTE: I'm using Boost 1.40 and GCC 4.4.1. I've compiled with:
> $ g++ -o lu -Wall -Wextra -pedantic -ansi lu.cpp
>
> Can you give me an help, please?
>
> Best,
>
> -- Marco

---

Your E-mail and More On-the-Go. Get Windows Live Hotmail Free. Sign up now.
_______________________________________________
ublas mailing list
ublas@...
http://lists.boost.org/mailman/listinfo.cgi/ublas
Sent to: lists@...

[Marco Guazzone]

On Tue, Sep 15, 2009 at 6:12 PM, Nasos Iliopoulos <nasos_i@...> wrote:

> Marco,
> lu_factorize(A) : A == L*U
> lu_factorize(A2, P) : P*A2 ==L*U
> P is a permutation matrix that interchanges the rows of A2. In some cases a
> matrix cannot be factored, so using pivoting with the permutation matrix can
> give you a factorisation. Also solving a system with and without pivoting
> has some slight differences.
>
> http://math.fullerton.edu/mathews/n2003/LUFactorMod.html
>
> I hope this helps.
>

Thanks Nasos!
In effect it seems that both Octave and Mathematica always use LU
decomposition with partial pivoting.
I've found a little Octave program luguass.m which implements LU
without pivoting.
It comes from the book:
  Quarteroni et al. "Scientific Computing with MATLAB and Octave" 2nd
edition (2006)
  http://mox.polimi.it/qs/
The result is the same one I get from uBlas. :)

For obtaining L and U matrices I would proceed as follows (sorry this
is Matlab/Octave code, hope do you understand):
L = tril(LU,-1) + eye(4)
U = triu(LU)
That is L is the lower-triangular matrix (without the main diagonal)
plus the identity matrix (of order 4) and U is the upper-triangular
matrix.
And in effect A2==L*U

Is there a better way to obtain them?

Thank you very much!

Cheers!

-- Marco
_______________________________________________
ublas mailing list
ublas@...
http://lists.boost.org/mailman/listinfo.cgi/ublas
Sent to: lists@...