This is the mail archive of the gsl-discuss@sources.redhat.com mailing list for the GSL project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

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
> >
> > *********************************************************
>


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]