[136] | 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;
|
---|
| 28 | static Matrix createMatrix(int rows, int cols) { return gsl_matrix_alloc(rows,cols); }
|
---|
| 29 | static Vector createVector(int size) { return gsl_vector_alloc(size); }
|
---|
| 30 | static void free(Matrix& m) { gsl_matrix_free(m); m=0; }
|
---|
| 31 | static void free(Vector& m) { gsl_vector_free(m); m=0; }
|
---|
| 32 | static void prod(const Matrix& m, const Vector& v, Vector& x) { gsl_blas_dgemv(CblasNoTrans,1,m,v,0,x); }
|
---|
| 33 | static void cholesky(Matrix& m) { gsl_linalg_cholesky_decomp(m); }
|
---|
| 34 | static void cholesky_solve(const Matrix& m, const Vector& b, Vector& x) { gsl_linalg_cholesky_solve(m,b,x); }
|
---|
| 35 | 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 | }
|
---|
| 45 | 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
|
---|