1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64
|
#include <stdio.h>
#include <gsl/gsl_linalg.h>
int
namespace dvp
{
template<class T, size_t N>
class matrix
{
std::array<T,N*N> _data;
public:
matrix() : _data({0}) {}//constructeur par défaut
matrix(size_t line=0, size_t column=0)//constructeur avec arguments
// copy & move constructors & operator= left to lecter
T& operator()(size_t line, size_t column) { return _data[N * column + line]; }
const T& operator()(size_t line, size_t column) const { return _data[N * column + line]; }
T* data() { return _data.data(); }
const T* data() const { return _data.data(); }
};
}
main (void)
{
matrix <double> a_data(class double, size_t 4);
double a_data(1,1)=1;
double a_data(1,2)=3;
double a_data(4,3)=5;
double a_data(4,4)=1;
for (int i=2;i<4;i++)
{
a(i,i)=1;
a(i,i+1)=3;
a(i+1,i)=5;
}
double b_data[] = { 1.0, 2.0, 3.0, 4.0 };
gsl_matrix_view m
= gsl_matrix_view_array (a_data, 4, 4);
gsl_vector_view b
= gsl_vector_view_array (b_data, 4);
gsl_vector *x = gsl_vector_alloc (4);
int s;
gsl_permutation * p = gsl_permutation_alloc (4);
gsl_linalg_LU_decomp (&m.matrix, p, &s);
gsl_linalg_LU_solve (&m.matrix, p, &b.vector, x);
printf ("x = \n");
gsl_vector_fprintf (stdout, x, "%g");
gsl_permutation_free (p);
gsl_vector_free (x);
return 0;
} |
Partager