help-gsl
[Top][All Lists]
Advanced

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

[Help-gsl] segmentation fault while using gsl_linalg_complex_LU_decompos


From: Владимир Дрынкин
Subject: [Help-gsl] segmentation fault while using gsl_linalg_complex_LU_decomposition
Date: Tue, 9 Aug 2011 02:04:02 +0400

Hello, all!!

Tonight i was trying to make a c++ class for working with gsl_complex
matrices. But when i compiled my program, i got segmentation fault
error :(
The code was built on the mac os x lion and debian 6 stable, but the
result was the same. Can you help me, what's wrong with my code?

Thanks a lot, Vladimir.

[CODE]
#include <iostream>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_linalg.h>
using namespace std;

class matrix
{
    gsl_matrix_complex *mtrx;
public:
    size_t rows;
    size_t cols;
    matrix();
    matrix(const size_t r, const size_t c);
    ~matrix();
    matrix(const matrix &rhs);
    matrix operator=(const matrix &rhs);
    matrix operator+(const matrix &rhs);
    matrix operator-(const matrix &rhs);
    matrix operator*(const matrix &rhs);
    void print();
    void randomize();
    matrix inverse();
    gsl_complex &at(size_t r, size_t c);
};

matrix::matrix()
{
    mtrx = gsl_matrix_complex_alloc(1, 1);
    rows = cols = 1;
    cout << "Created uninitialised default 1x1 matrix\n";
}

matrix::matrix(const size_t r, const size_t c)
{
    mtrx = gsl_matrix_complex_alloc(r, c);
    rows = r;
    cols = c;
    cout << "Created uninitialised " << rows << "x" << cols << " matrix\n";
}

matrix::~matrix()
{
    gsl_matrix_complex_free(mtrx);
    cout << "The matrix has been deallocated\n";
}

matrix::matrix(const matrix &rhs)
{
    mtrx = gsl_matrix_complex_alloc(rhs.rows, rhs.cols);
    rows=rhs.rows;
    cols=rhs.cols;
    int a = gsl_matrix_complex_memcpy(mtrx, rhs.mtrx);
    cout << a << " - the code, returned by copy constructor (0=Success)\n";
}

matrix matrix::operator=(const matrix &rhs)
{
    gsl_matrix_complex_memcpy(mtrx, rhs.mtrx);
    return *this;
}

matrix matrix::operator+(const matrix &rhs)
{
    matrix temp(rows, cols);
    temp = *this;
    gsl_matrix_complex_add(temp.mtrx, rhs.mtrx);
    return temp;
}

matrix matrix::operator-(const matrix &rhs)
{
    matrix temp(rows, cols);
    temp = *this;
    gsl_matrix_complex_sub(temp.mtrx, rhs.mtrx);
    return temp;
}

matrix matrix::operator*(const matrix &rhs)
{
    matrix temp(rows, rhs.cols);
    gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, gsl_complex_rect(1.0,
0.0), mtrx, rhs.mtrx, gsl_complex_rect(0.0, 0.0), temp.mtrx);
    return temp;
}

void matrix::print()
{
    gsl_matrix_complex_fprintf(stdout, mtrx, "%f");
}

void matrix::randomize()
{
    const gsl_rng_type *T;
    gsl_rng *r;
    T = gsl_rng_taus2;
    r = gsl_rng_alloc(T);
    for (size_t n=0; n<rows; n++)
    {
        for(size_t m=0; m<cols; m++)
        {
            this->at(n,m)=gsl_complex_rect((double)gsl_rng_get(r)*0.0000001,
(double)gsl_rng_get(r)*0.0000001);
        }
    }
    gsl_rng_free(r);
}

matrix matrix::inverse()
{
    int *signum;
    gsl_permutation *p;
    p = gsl_permutation_alloc(rows);
    matrix temp = *this;
    matrix temp_inv(rows, cols);
    gsl_linalg_complex_LU_decomp(temp.mtrx, p, signum);
    temp.print();
    gsl_linalg_complex_LU_invert(temp.mtrx, p, temp_inv.mtrx);
    return temp_inv;
}

gsl_complex &matrix::at(size_t r, size_t c)
{
    gsl_complex *p;
    p = gsl_matrix_complex_ptr(mtrx, r, c);
    if(p!=0) return *p;
}

int main (int argc, const char * argv[])
{
    matrix b(3,3);
    b.randomize();
    b.inverse();
    return 0;
}
[/CODE]



reply via email to

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