LU Decomposition error

[Peter Roche <P.J.P.Roche@damtp.cam.ac.uk>]

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

*********************************************************