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