1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra. Eigen itself is part of the KDE project.
3 //
4 // Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
5 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 
10 #ifndef EIGEN_GSL_HELPER
11 #define EIGEN_GSL_HELPER
12 
13 #include <Eigen/Core>
14 
15 #include <gsl/gsl_blas.h>
16 #include <gsl/gsl_multifit.h>
17 #include <gsl/gsl_eigen.h>
18 #include <gsl/gsl_linalg.h>
19 #include <gsl/gsl_complex.h>
20 #include <gsl/gsl_complex_math.h>
21 
22 namespace Eigen {
23 
24 template<typename Scalar, bool IsComplex = NumTraits<Scalar>::IsComplex> struct GslTraits
25 {
26   typedef gsl_matrix* Matrix;
27   typedef gsl_vector* Vector;
createMatrixGslTraits28   static Matrix createMatrix(int rows, int cols) { return gsl_matrix_alloc(rows,cols); }
createVectorGslTraits29   static Vector createVector(int size) { return gsl_vector_alloc(size); }
freeGslTraits30   static void free(Matrix& m) { gsl_matrix_free(m); m=0; }
freeGslTraits31   static void free(Vector& m) { gsl_vector_free(m); m=0; }
prodGslTraits32   static void prod(const Matrix& m, const Vector& v, Vector& x) { gsl_blas_dgemv(CblasNoTrans,1,m,v,0,x); }
choleskyGslTraits33   static void cholesky(Matrix& m) { gsl_linalg_cholesky_decomp(m); }
cholesky_solveGslTraits34   static void cholesky_solve(const Matrix& m, const Vector& b, Vector& x) { gsl_linalg_cholesky_solve(m,b,x); }
eigen_symmGslTraits35   static void eigen_symm(const Matrix& m, Vector& eval, Matrix& evec)
36   {
37     gsl_eigen_symmv_workspace * w = gsl_eigen_symmv_alloc(m->size1);
38     Matrix a = createMatrix(m->size1, m->size2);
39     gsl_matrix_memcpy(a, m);
40     gsl_eigen_symmv(a,eval,evec,w);
41     gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_VAL_ASC);
42     gsl_eigen_symmv_free(w);
43     free(a);
44   }
eigen_symm_genGslTraits45   static void eigen_symm_gen(const Matrix& m, const Matrix& _b, Vector& eval, Matrix& evec)
46   {
47     gsl_eigen_gensymmv_workspace * w = gsl_eigen_gensymmv_alloc(m->size1);
48     Matrix a = createMatrix(m->size1, m->size2);
49     Matrix b = createMatrix(_b->size1, _b->size2);
50     gsl_matrix_memcpy(a, m);
51     gsl_matrix_memcpy(b, _b);
52     gsl_eigen_gensymmv(a,b,eval,evec,w);
53     gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_VAL_ASC);
54     gsl_eigen_gensymmv_free(w);
55     free(a);
56   }
57 };
58 
59 template<typename Scalar> struct GslTraits<Scalar,true>
60 {
61   typedef gsl_matrix_complex* Matrix;
62   typedef gsl_vector_complex* Vector;
63   static Matrix createMatrix(int rows, int cols) { return gsl_matrix_complex_alloc(rows,cols); }
64   static Vector createVector(int size) { return gsl_vector_complex_alloc(size); }
65   static void free(Matrix& m) { gsl_matrix_complex_free(m); m=0; }
66   static void free(Vector& m) { gsl_vector_complex_free(m); m=0; }
67   static void cholesky(Matrix& m) { gsl_linalg_complex_cholesky_decomp(m); }
68   static void cholesky_solve(const Matrix& m, const Vector& b, Vector& x) { gsl_linalg_complex_cholesky_solve(m,b,x); }
69   static void prod(const Matrix& m, const Vector& v, Vector& x)
70   { gsl_blas_zgemv(CblasNoTrans,gsl_complex_rect(1,0),m,v,gsl_complex_rect(0,0),x); }
71   static void eigen_symm(const Matrix& m, gsl_vector* &eval, Matrix& evec)
72   {
73     gsl_eigen_hermv_workspace * w = gsl_eigen_hermv_alloc(m->size1);
74     Matrix a = createMatrix(m->size1, m->size2);
75     gsl_matrix_complex_memcpy(a, m);
76     gsl_eigen_hermv(a,eval,evec,w);
77     gsl_eigen_hermv_sort(eval, evec, GSL_EIGEN_SORT_VAL_ASC);
78     gsl_eigen_hermv_free(w);
79     free(a);
80   }
81   static void eigen_symm_gen(const Matrix& m, const Matrix& _b, gsl_vector* &eval, Matrix& evec)
82   {
83     gsl_eigen_genhermv_workspace * w = gsl_eigen_genhermv_alloc(m->size1);
84     Matrix a = createMatrix(m->size1, m->size2);
85     Matrix b = createMatrix(_b->size1, _b->size2);
86     gsl_matrix_complex_memcpy(a, m);
87     gsl_matrix_complex_memcpy(b, _b);
88     gsl_eigen_genhermv(a,b,eval,evec,w);
89     gsl_eigen_hermv_sort(eval, evec, GSL_EIGEN_SORT_VAL_ASC);
90     gsl_eigen_genhermv_free(w);
91     free(a);
92   }
93 };
94 
95 template<typename MatrixType>
96 void convert(const MatrixType& m, gsl_matrix* &res)
97 {
98 //   if (res)
99 //     gsl_matrix_free(res);
100   res = gsl_matrix_alloc(m.rows(), m.cols());
101   for (int i=0 ; i<m.rows() ; ++i)
102     for (int j=0 ; j<m.cols(); ++j)
103       gsl_matrix_set(res, i, j, m(i,j));
104 }
105 
106 template<typename MatrixType>
107 void convert(const gsl_matrix* m, MatrixType& res)
108 {
109   res.resize(int(m->size1), int(m->size2));
110   for (int i=0 ; i<res.rows() ; ++i)
111     for (int j=0 ; j<res.cols(); ++j)
112       res(i,j) = gsl_matrix_get(m,i,j);
113 }
114 
115 template<typename VectorType>
116 void convert(const VectorType& m, gsl_vector* &res)
117 {
118   if (res) gsl_vector_free(res);
119   res = gsl_vector_alloc(m.size());
120   for (int i=0 ; i<m.size() ; ++i)
121       gsl_vector_set(res, i, m[i]);
122 }
123 
124 template<typename VectorType>
125 void convert(const gsl_vector* m, VectorType& res)
126 {
127   res.resize (m->size);
128   for (int i=0 ; i<res.rows() ; ++i)
129     res[i] = gsl_vector_get(m, i);
130 }
131 
132 template<typename MatrixType>
133 void convert(const MatrixType& m, gsl_matrix_complex* &res)
134 {
135   res = gsl_matrix_complex_alloc(m.rows(), m.cols());
136   for (int i=0 ; i<m.rows() ; ++i)
137     for (int j=0 ; j<m.cols(); ++j)
138     {
139       gsl_matrix_complex_set(res, i, j,
140         gsl_complex_rect(m(i,j).real(), m(i,j).imag()));
141     }
142 }
143 
144 template<typename MatrixType>
145 void convert(const gsl_matrix_complex* m, MatrixType& res)
146 {
147   res.resize(int(m->size1), int(m->size2));
148   for (int i=0 ; i<res.rows() ; ++i)
149     for (int j=0 ; j<res.cols(); ++j)
150       res(i,j) = typename MatrixType::Scalar(
151         GSL_REAL(gsl_matrix_complex_get(m,i,j)),
152         GSL_IMAG(gsl_matrix_complex_get(m,i,j)));
153 }
154 
155 template<typename VectorType>
156 void convert(const VectorType& m, gsl_vector_complex* &res)
157 {
158   res = gsl_vector_complex_alloc(m.size());
159   for (int i=0 ; i<m.size() ; ++i)
160       gsl_vector_complex_set(res, i, gsl_complex_rect(m[i].real(), m[i].imag()));
161 }
162 
163 template<typename VectorType>
164 void convert(const gsl_vector_complex* m, VectorType& res)
165 {
166   res.resize(m->size);
167   for (int i=0 ; i<res.rows() ; ++i)
168     res[i] = typename VectorType::Scalar(
169         GSL_REAL(gsl_vector_complex_get(m, i)),
170         GSL_IMAG(gsl_vector_complex_get(m, i)));
171 }
172 
173 }
174 
175 #endif // EIGEN_GSL_HELPER
176