Subject: Re: [Boost-bugs] [Boost C++ Libraries] #11295: Matrix Memory problem after using lu_factorize
From: Boost C++ Libraries (noreply_at_[hidden])
Date: 2017-06-16 11:49:54
#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 anonymous):
Thanks for the explanation. Didn't realize that both L and U are combined
into a single matrix.
Thanks again, Michael
Replying to [comment:1 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.
>
>
> Replying to [ticket:11295 michael.cortis@…]:
> > 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/trac10/boost/ticket/11295#comment:2> 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-06-16 11:53:38 UTC