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