Loading...

help-gsl@gnu.org

[Prev] Thread [Next]  |  [Prev] Date [Next]

[Help-gsl] Help with gsl_linalg_complex_LU_det() function Ruben Dario Guerrero Mancilla Fri Feb 17 16:00:54 2012

Hello there,

Thanks in advance for reading.

I'm newbie using GSL. Currently I am trying to make a function to
which I pass a flatten matrix of doubles in form of a vector x[],
within the function I build again the matrix and assing it to the
object  "gsl_matrix * M = gsl_matrix_alloc (4, 4);", next to it I am
interested in calculate the quantity sqrt(Det(eye+i*M)) where "eye" is
the identity 4x4, To at future work with bigger matrices I use the
function  gsl_linalg_complex_LU_det(&C.matrix_complex,k). The function
must return the complex value  "sqrt(Det(eye+i*M)) "

The problem:

When I pass a vector with data X[] to the function and print the
result I find that always the function returns norm zero complex
numbers.

I'm sure tht I incurre in a missconception in the code. Some body can
tell my  where is my error in the function? and how can I correct it?


I attach the source code  of the function  here:
//****************************************************************************************
#include "/usr/include/boost/random.hpp"
#include "/usr/include/boost/math/distributions/normal.hpp"
#include <malloc.h>
#include <complex>
#include <math.h>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <cmath>
#include <stdio.h>
#include <stdlib.h>
#include "fftw3.h"

using namespace std;
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_cblas.h>

std::complex<double> CFFVV(double *x)
{
  double RC_aux, IC_aux,RC,IC;
  double c0=C0;
  double c1=C1;
  double g=G;
  gsl_matrix_complex *C = gsl_matrix_complex_alloc (4, 4);
  gsl_matrix * M = gsl_matrix_alloc (4, 4);
  gsl_matrix * one = gsl_matrix_alloc (4, 4);
  gsl_matrix_set_identity(one);
  for(int i=0; i<4;i++)
    {
      for(int j=0; j<4; j++)
        {
          gsl_matrix_set (M, i, j, x[5+i+j]);
        }
    }
  for(int i=0; i<4;i++)
    {
      for(int j=0; j<4; j++)
        {
          double ra=gsl_matrix_get(one,i,j);
          double ia= gsl_matrix_get(M,i,j);
          gsl_complex z;
          GSL_SET_COMPLEX(&z,ra,ia);
          gsl_matrix_complex_set (C, i, j,z );
        }
    }
  gsl_permutation * p = gsl_permutation_alloc (4);
  int k;
  gsl_complex det;
  det=gsl_linalg_complex_LU_det(&C.matrix_complex,k);
  cout<< k<<endl;
  return std::sqrt(std::complex<double>(GSL_REAL (det), GSL_IMAG (det)));
}

//*************************************************************************************************

-- 
Rubén Darío Guerrero M. Sc.
Estudiante Doctorado en Ciencias-Física (Universidad Nacional Bogota - Colombia)

Unidad Camilo Torres
Bloque 9 Oficina 2B2
Calle 44 #45-67
Bogotá - Colombia