gsl_helper.h
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra. Eigen itself is part of the KDE project.
00003 //
00004 // Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
00005 //
00006 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00024 
00025 #ifndef EIGEN_GSL_HELPER
00026 #define EIGEN_GSL_HELPER
00027 
00028 #include <Eigen/Core>
00029 
00030 #include <gsl/gsl_blas.h>
00031 #include <gsl/gsl_multifit.h>
00032 #include <gsl/gsl_eigen.h>
00033 #include <gsl/gsl_linalg.h>
00034 #include <gsl/gsl_complex.h>
00035 #include <gsl/gsl_complex_math.h>
00036 
00037 namespace Eigen {
00038 
00039 template<typename Scalar, bool IsComplex = NumTraits<Scalar>::IsComplex> struct GslTraits
00040 {
00041   typedef gsl_matrix* Matrix;
00042   typedef gsl_vector* Vector;
00043   static Matrix createMatrix(int rows, int cols) { return gsl_matrix_alloc(rows,cols); }
00044   static Vector createVector(int size) { return gsl_vector_alloc(size); }
00045   static void free(Matrix& m) { gsl_matrix_free(m); m=0; }
00046   static void free(Vector& m) { gsl_vector_free(m); m=0; }
00047   static void prod(const Matrix& m, const Vector& v, Vector& x) { gsl_blas_dgemv(CblasNoTrans,1,m,v,0,x); }
00048   static void cholesky(Matrix& m) { gsl_linalg_cholesky_decomp(m); }
00049   static void cholesky_solve(const Matrix& m, const Vector& b, Vector& x) { gsl_linalg_cholesky_solve(m,b,x); }
00050   static void eigen_symm(const Matrix& m, Vector& eval, Matrix& evec)
00051   {
00052     gsl_eigen_symmv_workspace * w = gsl_eigen_symmv_alloc(m->size1);
00053     Matrix a = createMatrix(m->size1, m->size2);
00054     gsl_matrix_memcpy(a, m);
00055     gsl_eigen_symmv(a,eval,evec,w);
00056     gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_VAL_ASC);
00057     gsl_eigen_symmv_free(w);
00058     free(a);
00059   }
00060   static void eigen_symm_gen(const Matrix& m, const Matrix& _b, Vector& eval, Matrix& evec)
00061   {
00062     gsl_eigen_gensymmv_workspace * w = gsl_eigen_gensymmv_alloc(m->size1);
00063     Matrix a = createMatrix(m->size1, m->size2);
00064     Matrix b = createMatrix(_b->size1, _b->size2);
00065     gsl_matrix_memcpy(a, m);
00066     gsl_matrix_memcpy(b, _b);
00067     gsl_eigen_gensymmv(a,b,eval,evec,w);
00068     gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_VAL_ASC);
00069     gsl_eigen_gensymmv_free(w);
00070     free(a);
00071   }
00072 };
00073 
00074 template<typename Scalar> struct GslTraits<Scalar,true>
00075 {
00076   typedef gsl_matrix_complex* Matrix;
00077   typedef gsl_vector_complex* Vector;
00078   static Matrix createMatrix(int rows, int cols) { return gsl_matrix_complex_alloc(rows,cols); }
00079   static Vector createVector(int size) { return gsl_vector_complex_alloc(size); }
00080   static void free(Matrix& m) { gsl_matrix_complex_free(m); m=0; }
00081   static void free(Vector& m) { gsl_vector_complex_free(m); m=0; }
00082   static void cholesky(Matrix& m) { gsl_linalg_complex_cholesky_decomp(m); }
00083   static void cholesky_solve(const Matrix& m, const Vector& b, Vector& x) { gsl_linalg_complex_cholesky_solve(m,b,x); }
00084   static void prod(const Matrix& m, const Vector& v, Vector& x)
00085   { gsl_blas_zgemv(CblasNoTrans,gsl_complex_rect(1,0),m,v,gsl_complex_rect(0,0),x); }
00086   static void eigen_symm(const Matrix& m, gsl_vector* &eval, Matrix& evec)
00087   {
00088     gsl_eigen_hermv_workspace * w = gsl_eigen_hermv_alloc(m->size1);
00089     Matrix a = createMatrix(m->size1, m->size2);
00090     gsl_matrix_complex_memcpy(a, m);
00091     gsl_eigen_hermv(a,eval,evec,w);
00092     gsl_eigen_hermv_sort(eval, evec, GSL_EIGEN_SORT_VAL_ASC);
00093     gsl_eigen_hermv_free(w);
00094     free(a);
00095   }
00096   static void eigen_symm_gen(const Matrix& m, const Matrix& _b, gsl_vector* &eval, Matrix& evec)
00097   {
00098     gsl_eigen_genhermv_workspace * w = gsl_eigen_genhermv_alloc(m->size1);
00099     Matrix a = createMatrix(m->size1, m->size2);
00100     Matrix b = createMatrix(_b->size1, _b->size2);
00101     gsl_matrix_complex_memcpy(a, m);
00102     gsl_matrix_complex_memcpy(b, _b);
00103     gsl_eigen_genhermv(a,b,eval,evec,w);
00104     gsl_eigen_hermv_sort(eval, evec, GSL_EIGEN_SORT_VAL_ASC);
00105     gsl_eigen_genhermv_free(w);
00106     free(a);
00107   }
00108 };
00109 
00110 template<typename MatrixType>
00111 void convert(const MatrixType& m, gsl_matrix* &res)
00112 {
00113 //   if (res)
00114 //     gsl_matrix_free(res);
00115   res = gsl_matrix_alloc(m.rows(), m.cols());
00116   for (int i=0 ; i<m.rows() ; ++i)
00117     for (int j=0 ; j<m.cols(); ++j)
00118       gsl_matrix_set(res, i, j, m(i,j));
00119 }
00120 
00121 template<typename MatrixType>
00122 void convert(const gsl_matrix* m, MatrixType& res)
00123 {
00124   res.resize(int(m->size1), int(m->size2));
00125   for (int i=0 ; i<res.rows() ; ++i)
00126     for (int j=0 ; j<res.cols(); ++j)
00127       res(i,j) = gsl_matrix_get(m,i,j);
00128 }
00129 
00130 template<typename VectorType>
00131 void convert(const VectorType& m, gsl_vector* &res)
00132 {
00133   if (res) gsl_vector_free(res);
00134   res = gsl_vector_alloc(m.size());
00135   for (int i=0 ; i<m.size() ; ++i)
00136       gsl_vector_set(res, i, m[i]);
00137 }
00138 
00139 template<typename VectorType>
00140 void convert(const gsl_vector* m, VectorType& res)
00141 {
00142   res.resize (m->size);
00143   for (int i=0 ; i<res.rows() ; ++i)
00144     res[i] = gsl_vector_get(m, i);
00145 }
00146 
00147 template<typename MatrixType>
00148 void convert(const MatrixType& m, gsl_matrix_complex* &res)
00149 {
00150   res = gsl_matrix_complex_alloc(m.rows(), m.cols());
00151   for (int i=0 ; i<m.rows() ; ++i)
00152     for (int j=0 ; j<m.cols(); ++j)
00153     {
00154       gsl_matrix_complex_set(res, i, j,
00155         gsl_complex_rect(m(i,j).real(), m(i,j).imag()));
00156     }
00157 }
00158 
00159 template<typename MatrixType>
00160 void convert(const gsl_matrix_complex* m, MatrixType& res)
00161 {
00162   res.resize(int(m->size1), int(m->size2));
00163   for (int i=0 ; i<res.rows() ; ++i)
00164     for (int j=0 ; j<res.cols(); ++j)
00165       res(i,j) = typename MatrixType::Scalar(
00166         GSL_REAL(gsl_matrix_complex_get(m,i,j)),
00167         GSL_IMAG(gsl_matrix_complex_get(m,i,j)));
00168 }
00169 
00170 template<typename VectorType>
00171 void convert(const VectorType& m, gsl_vector_complex* &res)
00172 {
00173   res = gsl_vector_complex_alloc(m.size());
00174   for (int i=0 ; i<m.size() ; ++i)
00175       gsl_vector_complex_set(res, i, gsl_complex_rect(m[i].real(), m[i].imag()));
00176 }
00177 
00178 template<typename VectorType>
00179 void convert(const gsl_vector_complex* m, VectorType& res)
00180 {
00181   res.resize(m->size);
00182   for (int i=0 ; i<res.rows() ; ++i)
00183     res[i] = typename VectorType::Scalar(
00184         GSL_REAL(gsl_vector_complex_get(m, i)),
00185         GSL_IMAG(gsl_vector_complex_get(m, i)));
00186 }
00187 
00188 }
00189 
00190 #endif // EIGEN_GSL_HELPER


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:32:45