From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 18052 invoked by alias); 16 Nov 2002 22:18:22 -0000 Mailing-List: contact gsl-discuss-help@sources.redhat.com; run by ezmlm Precedence: bulk List-Subscribe: List-Archive: List-Post: List-Help: , Sender: gsl-discuss-owner@sources.redhat.com Received: (qmail 18045 invoked from network); 16 Nov 2002 22:18:19 -0000 Received: from unknown (HELO axon.amtp.cam.ac.uk) (131.111.16.133) by sources.redhat.com with SMTP; 16 Nov 2002 22:18:19 -0000 Received: from narmada.amtp.cam.ac.uk ([131.111.18.45] ident=root) by axon.amtp.cam.ac.uk with esmtp (Exim 3.36 #1) id 18DBGc-0000Rv-00 for gsl-discuss@sources.redhat.com; Sat, 16 Nov 2002 22:18:18 +0000 Received: from localhost by narmada.amtp.cam.ac.uk (8.11.6/DAMTP 3.1-Client (damtp.cam.ac.uk)) id gAGMIIq28944; Sat, 16 Nov 2002 22:18:18 GMT Date: Mon, 18 Nov 2002 08:55:00 -0000 From: Peter Roche X-X-Sender: pjpr2@narmada.amtp.cam.ac.uk To: gsl-discuss@sources.redhat.com Subject: Re: LU Decomposition error In-Reply-To: <20021113201143.B7FD94EC9F@acmez.gatech.edu> Message-ID: MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII X-SW-Source: 2002-q4/txt/msg00155.txt.bz2 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 > for(inb=ina;inb > if(inb > ina ) gsl_matrix_set(sym_denom_Lower,ina,inb,0.0); > > } > > for(inb=ina+1;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 > > > > ********************************************************* >