From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 6117 invoked by alias); 13 Nov 2002 20:11:48 -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 6109 invoked from network); 13 Nov 2002 20:11:48 -0000 Received: from unknown (HELO hellmouth5.gatech.edu) (130.207.165.165) by sources.redhat.com with SMTP; 13 Nov 2002 20:11:48 -0000 Received: from hellmouth5.gatech.edu (localhost [127.0.0.1]) by hellmouth5.gatech.edu (Postfix) with SMTP id 181E91D1A22; Wed, 13 Nov 2002 15:11:48 -0500 (EST) (envelope-from peles@cns.physics.gatech.edu) Received: from acmez.gatech.edu (acmez.prism.gatech.edu [130.207.171.28]) by hellmouth5.gatech.edu (Postfix) with ESMTP id E35691D1A1B; Wed, 13 Nov 2002 15:11:47 -0500 (EST) (envelope-from peles@cns.physics.gatech.edu) Received: from there (kaos.physics.gatech.edu [130.207.140.139]) by acmez.gatech.edu (Postfix) with SMTP id B7FD94EC9F; Wed, 13 Nov 2002 15:11:43 -0500 (EST) (envelope-from peles@cns.physics.gatech.edu) Content-Type: text/plain; charset="iso-8859-1" From: Slaven Peles To: Peter Roche , gsl-discuss@sources.redhat.com Subject: Re: LU Decomposition error Date: Thu, 14 Nov 2002 03:11:00 -0000 References: In-Reply-To: MIME-Version: 1.0 Content-Transfer-Encoding: 8bit Message-Id: <20021113201143.B7FD94EC9F@acmez.gatech.edu> X-SW-Source: 2002-q4/txt/msg00147.txt.bz2 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 > > *********************************************************