* LU Decomposition error
@ 2002-11-13 12:11 Peter Roche
2002-11-14 3:11 ` Slaven Peles
2002-11-16 11:05 ` Brian Gough
0 siblings, 2 replies; 4+ messages in thread
From: Peter Roche @ 2002-11-13 12:11 UTC (permalink / raw)
To: gsl-discuss
Dear All,
I hope this isn't a repeat of an already dealt with problem, but I
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
------------------------------------------------------------------------
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_Lower,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
*********************************************************
^ permalink raw reply [flat|nested] 4+ messages in thread* Re: LU Decomposition error 2002-11-13 12:11 LU Decomposition error Peter Roche @ 2002-11-14 3:11 ` Slaven Peles 2002-11-18 8:55 ` Peter Roche 2002-11-16 11:05 ` Brian Gough 1 sibling, 1 reply; 4+ messages in thread From: Slaven Peles @ 2002-11-14 3:11 UTC (permalink / raw) To: Peter Roche, gsl-discuss 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 > > ********************************************************* ^ permalink raw reply [flat|nested] 4+ messages in thread
* Re: LU Decomposition error 2002-11-14 3:11 ` Slaven Peles @ 2002-11-18 8:55 ` Peter Roche 0 siblings, 0 replies; 4+ messages in thread From: Peter Roche @ 2002-11-18 8:55 UTC (permalink / raw) To: gsl-discuss 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 > > > > ********************************************************* > ^ permalink raw reply [flat|nested] 4+ messages in thread
* Re: LU Decomposition error 2002-11-13 12:11 LU Decomposition error Peter Roche 2002-11-14 3:11 ` Slaven Peles @ 2002-11-16 11:05 ` Brian Gough 1 sibling, 0 replies; 4+ messages in thread From: Brian Gough @ 2002-11-16 11:05 UTC (permalink / raw) To: Peter Roche; +Cc: gsl-discuss, bug-gsl Peter Roche writes: > I hope this isn't a repeat of an already dealt with problem, but I > 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.... Please send a complete example program which reproduces the problem. Thanks. ^ permalink raw reply [flat|nested] 4+ messages in thread
end of thread, other threads:[~2002-11-16 22:18 UTC | newest] Thread overview: 4+ messages (download: mbox.gz / follow: Atom feed) -- links below jump to the message on this page -- 2002-11-13 12:11 LU Decomposition error Peter Roche 2002-11-14 3:11 ` Slaven Peles 2002-11-18 8:55 ` Peter Roche 2002-11-16 11:05 ` Brian Gough
This is a public inbox, see mirroring instructions for how to clone and mirror all data and code used for this inbox; as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).