This is the mail archive of the
gsl-discuss@sources.redhat.com
mailing list for the GSL project.
Re: LU Decomposition error
- From: Peter Roche <P dot J dot P dot Roche at damtp dot cam dot ac dot uk>
- To: gsl-discuss at sources dot redhat dot com
- Date: Sat, 16 Nov 2002 22:18:18 +0000 (GMT)
- Subject: Re: LU Decomposition error
On more careful reflection the LU decomposition seems OK. I had simply
not taken into consideration that the marix product of the Upper and Lower
triangular matrices gives a permutation of the original matrix, that is
P.A = L.U, where A is the original matrix, L the lower triangular matrix,
U the Upper triangular and P a vector containing a record of the row
permutations carried out in the decomposition.
Cheers, Peter
*********************************************************
Peter Roche
Department of Applied Mathematics and Theoretical Physics,
University of Cambridge, Silver Street, Cambridge, CB3 9EW
email: p.j.p.roche@damtp.cam.ac.uk
*********************************************************
On Wed, 13 Nov 2002, Slaven Peles wrote:
> On Wednesday 13 November 2002 02:39 pm, Peter Roche wrote:
> > couldn't see anything in the archives. I have encountered some problems
> > with the gsl_linalg_LU_decomp rountine. The routine ran without exiting
> > with an error message, but I found that if I multiply the decomposed
> > Lower and Upper triangular matrices I do not get the original matrix I
> > started with before decomposing. I have included a section of my code
> > and also a section of the original matrix (sym_denom) and the result of
> > matrix multiplying the decomposed lower and upper triangular matrices
> > (LU check).
> >
> > Some of the rows are the same as the original but are not in the correct
> > row. For example row 4 in "LU check" has the same elements as row 2
> > "sym_denom"!
> >
> > Has anyone else encountered similar problems or can help me with this
> > one.
> >
> > Cheers, Peter
> >
>
> Have you checked if your matrix is singular? It was reported before that LU
> decomposition does not exit with error when the input matrix is singular
> matrix. I've just had a look at function gsl_linalg_LU_decomp in lu.c in GSL
> 1.2. It does check for zeros at the main diagonal after the decomposition,
> and proceeds with calculation in case there are none. But it seems to me that
> it does nothing if it finds a zero, i.e. when the input matrix is singular. I
> copy that piece of gsl_linalg_LU_decomp code below.
>
> ======== from gsl_linalg_LU_decomp in lu.c ==================
>
> ajj = gsl_matrix_get (A, j, j);
>
> if (ajj != 0.0)
> {
> for (i = j + 1; i < N; i++)
> {
> REAL aij = gsl_matrix_get (A, i, j) / ajj;
> gsl_matrix_set (A, i, j, aij);
>
> for (k = j + 1; k < N; k++)
> {
> REAL aik = gsl_matrix_get (A, i, k);
> REAL ajk = gsl_matrix_get (A, j, k);
> gsl_matrix_set (A, i, k, aik - aij * ajk);
> }
> }
> }
> }
>
> return GSL_SUCCESS;
> }
> }
> ==========================================
>
> I guess there should be an else statement at the end calling GSL_ERROR.
>
> Cheers,
> Slaven
>
>
> > sym_denom_LU = gsl_matrix_calloc(num_ind,num_ind);
> > sym_denom_Lower = gsl_matrix_calloc(num_ind,num_ind);
> > sym_denom_Upper = gsl_matrix_calloc(num_ind,num_ind);
> >
> > perm = gsl_permutation_calloc(num_ind);
> > gsl_matrix_memcpy(sym_denom_LU,sym_denom);
> > gsl_linalg_LU_decomp(sym_denom_LU,perm,&sign);
> >
> > gsl_matrix_memcpy(sym_denom_Lower,sym_denom_LU);
> > gsl_matrix_memcpy(sym_denom_Upper,sym_denom_LU);
> >
> > /* check on LU decomp */
> > for(ina=0;ina<num_ind;ina++){
> > for(inb=ina;inb<num_ind;inb++){
> > if(inb > ina ) gsl_matrix_set(sym_denom_Lower,ina,inb,0.0);
> > }
> > for(inb=ina+1;inb<num_ind;inb++){
> > gsl_matrix_set(sym_denom_Upper,inb,ina,0.0);
> > }
> > }
> >
> > /* Multiply lower/upper matrices to get original sym_denom */
> >
> > gsl_blas_dtrmm(CblasLeft,CblasLower,CblasNoTrans,CblasUnit,1.0,sym_denom_Lo
> >wer,sym_denom_Upper);
> >
> > ------------------------------------------------------------------------
> >
> >
> > sym_denom
> > 1 -7.920023e-19 -1.652853e-19 -1.560040e-19 -6.114344e-20
> > -3.066652e-19 4.391243e-19 2 -1.652853e-19 -3.923216e-20
> > -3.608558e-20 -1.378546e-20 -7.772215e-20 9.616245e-20 3
> > -1.560040e-19 -3.608558e-20 -8.714679e-20 -2.635082e-20 -3.128613e-20
> > 7.866459e-20 4 -6.114344e-20 -1.378546e-20 -2.635082e-20
> > -9.419878e-21 -1.627012e-20 3.311913e-20 5 -3.066652e-19
> > -7.772215e-20 -3.128613e-20 -1.627012e-20 -2.676652e-19 2.252427e-19 6
> > 4.391243e-19 9.616245e-20 7.866459e-20 3.311913e-20 2.252427e-19
> > -3.018037e-19 7 -1.748050e-20 -3.759400e-21 -2.671663e-20
> > -8.956211e-21 2.023082e-20 4.257951e-21 8 2.556328e-19 6.230603e-20
> > 1.710786e-19 5.691356e-20 3.972437e-20 -1.485747e-19 9
> > 4.908470e-20 1.223947e-20 3.696976e-20 1.238438e-20 3.418982e-21
> > -2.760516e-20
> >
> >
> > LU check - should be the same as sym_denom
> > 1 -7.920023e-19 -1.652853e-19 -1.560040e-19 -6.114344e-20
> > -3.066652e-19 4.391243e-19 2 3.491029e-19 8.709227e-20
> > 5.217473e-20 2.250474e-20 2.709399e-19 -2.325323e-19 3 2.556328e-19
> > 6.230603e-20 1.710786e-19 5.691356e-20 3.972437e-20 -1.485747e-19 4
> > -1.560040e-19 -3.608558e-20 -8.714679e-20 -2.635082e-20 -3.128613e-20
> > 7.866459e-20 5 -1.652853e-19 -3.923216e-20 -3.608558e-20
> > -1.378546e-20 -7.772215e-20 9.616245e-20 6 2.847849e-19
> > 6.186536e-20 4.510521e-20 2.035696e-20 1.565680e-19 -2.154436e-19 7
> > -1.748050e-20 -3.759400e-21 -2.671663e-20 -8.956211e-21 2.023082e-20
> > 4.257951e-21 8 -1.493627e-19 -4.009839e-20 -1.313729e-19
> > -4.382593e-20 -5.583881e-21 8.188773e-20 9 1.981133e-20
> > 5.659666e-21 2.816777e-20 9.600430e-21 -8.449448e-21 -8.121479e-21
> >
> > *********************************************************
> >
> > Peter Roche
> > Department of Applied Mathematics and Theoretical Physics,
> > University of Cambridge, Silver Street, Cambridge, CB3 9EW
> >
> > email: p.j.p.roche@damtp.cam.ac.uk
> >
> > *********************************************************
>