Re: [Help-gsl] Determinant of a matrix:

[Nomesh Bolia]

Thanks a lot for the reply.

I am sorry, I did copy the matrix A into tmp_ptr. So, the code is like: 
(exactly as you suggested in the reply):

double get_det(gsl_matrix * A) {
  int sign=0; double det=0.0; int row_sq = A->size1;
  gsl_permutation * p = gsl_permutation_calloc(row_sq);
  gsl_matrix * tmp_ptr = gsl_matrix_calloc(row_sq, row_sq);
  int * signum = &sign;
  gsl_matrix_memcpy(tmp_ptr, A);
  gsl_linalg_LU_decomp(tmp_ptr, p, signum);
  det = gsl_linalg_LU_det(tmp_ptr, *signum);
  gsl_permutation_free(p);
  gsl_matrix_free(tmp_ptr);
  return det;
}

It doesnt work, even for simple 2X2 or 3X3 matrices. 

-Nomesh

On Fri, 11 Feb 2005, 
Joakim Hove wrote:

> Nomesh Bolia <address@hidden> writes:
> 
> > Hi,
> >
> > I am using the following program to compute the determinant of a matrix A:
> >
> > double get_det(gsl_matrix * A) {
> >   int sign=0; double det=0.0; int row_sq = mat_ptr->size1;
> >   gsl_permutation * p = gsl_permutation_calloc(row_sq);
> >   gsl_matrix * tmp_ptr = gsl_matrix_calloc(row_sq, row_sq);
> >   int * signum = &sign;
> >   gsl_matrix_memcpy(tmp_ptr, mat_ptr);
> >   gsl_linalg_LU_decomp(tmp_ptr, p, signum);
> >   det = gsl_linalg_LU_det(tmp_ptr, *signum);
> >   gsl_permutation_free(p);
> >   gsl_matrix_free(tmp_ptr);
> >   return det;
> > }
> >
> > Basically, the functions:
> >   
> > gsl_linalg_LU_decomp(tmp_ptr, p, signum); and 
> > det = gsl_linalg_LU_det(tmp_ptr, *signum);
> >
> > have been used to compute the determinant. However, I get completely 
> > arbitrary answers. 
> 
> 
> Well,
> 
> that seems reasonable. You are determining the LU decomposition and det of
> the newly allocated matrix pointed to by tmp_ptr, which is seemingly
> not related to A - are the arguments of the gsl_matrix_memcpy() call
> correct?
> 
> Why don't you just use A - I assume that is the matrix you are
> interested in. Alternative code (not tested):
> 
> 
> double det_get(gsl_matrix *A, int inPlace) {
> 
> /*
>   inPlace = 1 => A is replaced with the LU decomposed copy.
>   inPlace = 0 => A is retained, and a copy is used for LU.
> */
> 
>    double det;
>    int signum;
>    gsl_permutation *p = gsl_perumtation_alloc(A->size1);
>    gsl_matrix *tmpA;
> 
>    if (inPlace)
>       tmpA = A;
>    else {
>      gsl_matrix *tmpA = gsl_matrix_alloc(A->size1, A->size2);
>      gsl_matrix_memcpy(tmpA , A);
>    }
> 
> 
>    gsl_linalg_LU_decomp(tmpA , p , &signum);
>    det = gsl_linalg_LU_det(tmpA , signum);
>    gsl_permutation_free(p);
>    if (! inPlace)
>       gsl_matrix_free(tmpA);
> 
>    
>    return det;
> }
> 
> 
> HTH - Joakim
> 
> 
> 
> 

-- 
Regards,
Nomesh.

The world is a dangerous place to live; not because of the people who are 
evil, but because of the people who don't do anything about it. -- Albert 
Einstein.