# Re: [Boost-bugs] [Boost C++ Libraries] #11295: Matrix Memory problem after using lu_factorize

Subject: Re: [Boost-bugs] [Boost C++ Libraries] #11295: Matrix Memory problem after using lu_factorize
Date: 2017-03-31 10:45:26

#11295: Matrix Memory problem after using lu_factorize
-------------------------------+---------------------
Reporter: michael.cortis@â€¦ | Owner: guwi17
Type: Bugs | Status: new
Milestone: To Be Determined | Component: uBLAS
Version: Boost 1.55.0 | Severity: Problem
Resolution: | Keywords:
-------------------------------+---------------------

Comment (by 2015csb1032@â€¦):

You interpreted the function in the wrong way, that's what the
lu_factorize function does.
The return value, 0 or non 0 value which you see, only specifies whether
the matrix is singular or not.
The return value is 0 if the matix is singular, and non 0 if the matrix is
non-singular.

Now intuitively, you would think that lu_factorize should give two
matrices L and U on decomposition.
That's what is done in a way.

The matrix which you pass in lu_factorize(X), here X is passed by
reference to lu_factorize function,which will '''decompose X and now, X
will contain the information of the two factorized matrix L and U. Hence X
is changed.'''[[BR]]
If you don't want X to be changed. Then, create another matrix Y, do Y=X
and then pass Y in lu_factorize : lu_factorize(Y)

It's a typical approach in LU decomposition, that the diagonal elements of
L contains 1. So, the two matrices
L and U are fitted into 1 matrix Y, putting L in the lower part omitting
the diagonal elements, and
putting U in the upper part. L and U can be extracted from Y easily (shown
in code).

If, Y (after function ends) = [[BR]]
y1 y2 y3[[BR]]
y4 y5 y6[[BR]]
y7 y8 y9[[BR]]

Then, L=[[BR]]
1 0 0[[BR]]
y4 1 0[[BR]]
y7 y7 1[[BR]]

and u =[[BR]]
y1 y2 y3[[BR]]
0 y5 y6[[BR]]
0 0 y9[[BR]]

You can run the following code, it will show what I'm trying to say above.
{{{
ublas::matrix<double> X(3,3); X.clear();
X(0,0) = 0.995182407377577; X(0,1) =-0.006473367705848; X(0,2)
=-0.002032391957706;
X(1,0) =-0.006473367705848; X(1,1) = 0.995182407377577; X(1,2)
=-0.002032391957706;
X(2,0) =-0.002032391957706; X(2,1) =-0.002032391957706; X(2,2) =
0.936175146339137;

ublas::matrix<double> Y(3,3);
Y=X;
ublas::lu_factorize(Y);

ublas::matrix<double> L(3,3); L.clear();
ublas::matrix<double> U(3,3); U.clear();

for (int i=0; i<3; i++){
for (int j=0; j<3; j++){
if (i<j){
U(i,j) = Y(i,j);
L(i,j) = 0;
}
else if (i>j){
L(i,j) = Y(i,j);
U(i,j) = 0;
}
else if (i==j){
L(i,j) = 1;
U(i,j) = Y(i,j);
}
}
}

cout << "The original matrix X :"<<endl;
cout<<setprecision(16)<<X<<endl<<endl;

cout << "LU_decomposed matrix in 1 matrix :"<<endl;
cout<<setprecision(16)<<Y<<endl<<endl;

cout << "LU_decomposed matrices in 2 diff matrix, L :"<<endl;
cout<<setprecision(16)<<L<<endl<<endl;

cout << "LU_decomposed matrices in 2 diff matrix, U :"<<endl;
cout<<setprecision(16)<<U<<endl<<endl;

ublas::matrix<double> Z(3,3);
cout << "The product Z=L*U, which should be equal to X (Z=L*U=X),
Z:"<<endl;
axpy_prod(L, U, Z, true);
cout<<setprecision(16)<<Z<<endl<<endl;
}}}

Also, run the code with
{{{
X(0,0) = 1.0; X(0,1) = 2.0; X(0,2) = 3.0;
X(1,0) = 4.0; X(1,1) = 5.0; X(1,2) = 6.0;
X(2,0) = 7.0; X(2,1) = 8.0; X(2,2) = 9.0;
}}}
instead of X(0,0) = 0.995182407377577;... for better understanding.

> Discovered that after using lu_factorize the values of the matrix have
changed!
>
>
> {{{
> #include <iostream>
>
> #include <boost/date_time/posix_time/posix_time.hpp>
> #include <boost/numeric/ublas/matrix.hpp>
> #include <boost/numeric/ublas/matrix_proxy.hpp>
> #include <boost/numeric/ublas/vector.hpp>
> #include <boost/numeric/ublas/io.hpp>
> #include <boost/numeric/ublas/lu.hpp>
>
> using namespace std;
> using namespace boost::numeric;
>
> int main()
> {
> ublas::matrix<double> X(3,3); X.clear();
> X(0,0) = 0.995182407377577; X(0,1) =-0.006473367705848; X(0,2)
=-0.002032391957706;
> X(1,0) =-0.006473367705848; X(1,1) = 0.995182407377577; X(1,2)
=-0.002032391957706;
> X(2,0) =-0.002032391957706; X(2,1) =-0.002032391957706; X(2,2) =
0.936175146339137;
> cout<<setprecision(16)<<X<<endl;
> cout<<ublas::lu_factorize(X)<<endl;
> cout<<setprecision(16)<<X<<endl;
> cout<<ublas::lu_factorize(X)<<endl;
> cout<<setprecision(16)<<X<<endl;
> }
> }}}
>
> Output:
>
>
> {{{
>
[3,3]((0.995182407377577,-0.006473367705848,-0.002032391957706),(-0.006473367705848,0.995182407377577,-0.002032391957706),(-0.002032391957706,-0.002032391957706,0.936175146339137))
> 0
>
[3,3]((0.995182407377577,-0.006473367705848,-0.002032391957706),(-0.006504704723334175,0.9951403000320849,-0.002045612067272957),(-0.002042230592742885,-0.002055601674665375,0.9361667907625133))
> 0
>
[3,3]((0.995182407377577,-0.006473367705848,-0.002032391957706),(-0.006536193440632496,0.9950979888485471,-0.002058896174255709),(-0.002052116855767581,-0.002079077442455779,0.9361583394521271))
>
> }}}
>
> Is this fixed in newer versions?
>
> Thanks
>
> Michael Cortis
>
> RA, Durham University

```--
Ticket URL: <https://svn.boost.org/trac/boost/ticket/11295#comment:1>
Boost C++ Libraries <http://www.boost.org/>
Boost provides free peer-reviewed portable C++ source libraries.
```

This archive was generated by hypermail 2.1.7 : 2017-03-31 10:50:23 UTC