Cloudy & Associates

Commit 324a9a09 authored by M Chatzikos's avatar M Chatzikos
Browse files

Add solve of linear system with LU-decomposition

parent 949c1be9
......@@ -21,6 +21,8 @@ using namespace std;
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_linalg.h>
......@@ -191,6 +193,17 @@ static gsl_matrix *copy_matrix( const int naxis, const vector<double> &matrix,
return gmat;
}
static gsl_vector *copy_vector( const int naxis, const vector<double> &vec )
{
gsl_vector *gvec = gsl_vector_calloc( naxis );
for( int icol = 0; icol < naxis; icol++ )
{
gsl_vector_set( gvec, icol, vec[ icol ] );
}
return gvec;
}
static void calc_eigen( const int naxis, const vector<double> &matrix,
bool compute_Schur = false, bool balance = true )
{
......@@ -291,6 +304,58 @@ static void calc_eigen( const int naxis, const vector<double> &matrix,
gsl_eigen_nonsymm_free( ws );
}
static void solve_system( const int naxis, const vector<double> &matrix,
const vector<double> &creation,
const vector<double> &pop )
{
// Solve the system Ax=b for x
//
gsl_matrix *A_save = copy_matrix( naxis, matrix );
gsl_matrix *A = copy_matrix( naxis, matrix );
gsl_vector *b = copy_vector( naxis, creation );
gsl_permutation *p = gsl_permutation_alloc( naxis );
int n;
gsl_linalg_LU_decomp( A, p, &n );
cout << endl;
cout << "Determinant: " << gsl_linalg_LU_det( A, n ) << endl;
cout << "Determinant sign: " << gsl_linalg_LU_sgndet( A, n ) << endl;
gsl_vector *x = gsl_vector_calloc( naxis );
gsl_linalg_LU_solve( A, p, b, x );
gsl_vector *work = gsl_vector_calloc( naxis );
gsl_linalg_LU_refine( A_save, A, p, b, x, work );
bool prtHeader = false;
for( int ip = 0; ip < naxis; ip++ )
{
double np = gsl_vector_get( x, ip );
double diff = np - pop[ip];
if( fabs( diff ) > 1e-4 * pop[ip] )
{
if( not prtHeader )
{
cout << endl;
cout << "Compare GSL v Cloudy results:" << endl;
prtHeader = true;
}
cout << "\t" << np << "\t" << pop[ip] <<
"\t" << diff << endl;
}
}
gsl_matrix_free( A_save );
gsl_matrix_free( A );
gsl_vector_free( b );
gsl_permutation_free( p );
gsl_vector_free( x );
gsl_vector_free( work );
}
int main( const int argc, const char **argv )
{
string fitsFile;
......@@ -303,5 +368,7 @@ int main( const int argc, const char **argv )
calc_eigen( axes[0], matrix );
solve_system( axes[0], matrix, creation, pop );
return EXIT_SUCCESS;
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment