Re: LU Decomposition error

[Slaven Peles <peles@cns.physics.gatech.edu>]

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