help-gsl
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

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


From: Joakim Hove
Subject: Re: [Help-gsl] Determinant of a matrix:
Date: Fri, 11 Feb 2005 09:55:01 +0100
User-agent: Gnus/5.1006 (Gnus v5.10.6) Emacs/21.2 (gnu/linux)

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



-- 
Joakim Hove
hove AT ift uib no
Tlf: +47 (55 5)8 27 90 
Fax: +47 (55 5)8 94 40
http://www.ift.uib.no/~hove/





reply via email to

[Prev in Thread] Current Thread [Next in Thread]