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.
00003 //
00004 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.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 #include <gsl/gsl_poly.h>
00037 
00038 namespace Eigen {
00039 
00040 template<typename Scalar, bool IsComplex = NumTraits<Scalar>::IsComplex> struct GslTraits
00041 {
00042   typedef gsl_matrix* Matrix;
00043   typedef gsl_vector* Vector;
00044   static Matrix createMatrix(int rows, int cols) { return gsl_matrix_alloc(rows,cols); }
00045   static Vector createVector(int size) { return gsl_vector_alloc(size); }
00046   static void free(Matrix& m) { gsl_matrix_free(m); m=0; }
00047   static void free(Vector& m) { gsl_vector_free(m); m=0; }
00048   static void prod(const Matrix& m, const Vector& v, Vector& x) { gsl_blas_dgemv(CblasNoTrans,1,m,v,0,x); }
00049   static void cholesky(Matrix& m) { gsl_linalg_cholesky_decomp(m); }
00050   static void cholesky_solve(const Matrix& m, const Vector& b, Vector& x) { gsl_linalg_cholesky_solve(m,b,x); }
00051   static void eigen_symm(const Matrix& m, Vector& eval, Matrix& evec)
00052   {
00053     gsl_eigen_symmv_workspace * w = gsl_eigen_symmv_alloc(m->size1);
00054     Matrix a = createMatrix(m->size1, m->size2);
00055     gsl_matrix_memcpy(a, m);
00056     gsl_eigen_symmv(a,eval,evec,w);
00057     gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_VAL_ASC);
00058     gsl_eigen_symmv_free(w);
00059     free(a);
00060   }
00061   static void eigen_symm_gen(const Matrix& m, const Matrix& _b, Vector& eval, Matrix& evec)
00062   {
00063     gsl_eigen_gensymmv_workspace * w = gsl_eigen_gensymmv_alloc(m->size1);
00064     Matrix a = createMatrix(m->size1, m->size2);
00065     Matrix b = createMatrix(_b->size1, _b->size2);
00066     gsl_matrix_memcpy(a, m);
00067     gsl_matrix_memcpy(b, _b);
00068     gsl_eigen_gensymmv(a,b,eval,evec,w);
00069     gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_VAL_ASC);
00070     gsl_eigen_gensymmv_free(w);
00071     free(a);
00072   }
00073 
00074   template<class EIGEN_VECTOR, class EIGEN_ROOTS>
00075   static void eigen_poly_solve(const EIGEN_VECTOR& poly, EIGEN_ROOTS& roots )
00076   {
00077     const int deg = poly.size()-1;
00078     double *z = new double[2*deg];
00079     double *a = new double[poly.size()];
00080     for( int i=0; i<poly.size(); ++i ){ a[i] = poly[i]; }
00081 
00082     gsl_poly_complex_workspace * w = gsl_poly_complex_workspace_alloc (poly.size());
00083     gsl_poly_complex_solve(a, poly.size(), w, z);
00084     gsl_poly_complex_workspace_free (w);
00085     for( int i=0; i<deg; ++i )
00086     {
00087       roots[i].real() = z[2*i];
00088       roots[i].imag() = z[2*i+1];
00089     }
00090 
00091     delete[] a;
00092     delete[] z;
00093   }
00094 };
00095 
00096 template<typename Scalar> struct GslTraits<Scalar,true>
00097 {
00098   typedef gsl_matrix_complex* Matrix;
00099   typedef gsl_vector_complex* Vector;
00100   static Matrix createMatrix(int rows, int cols) { return gsl_matrix_complex_alloc(rows,cols); }
00101   static Vector createVector(int size) { return gsl_vector_complex_alloc(size); }
00102   static void free(Matrix& m) { gsl_matrix_complex_free(m); m=0; }
00103   static void free(Vector& m) { gsl_vector_complex_free(m); m=0; }
00104   static void cholesky(Matrix& m) { gsl_linalg_complex_cholesky_decomp(m); }
00105   static void cholesky_solve(const Matrix& m, const Vector& b, Vector& x) { gsl_linalg_complex_cholesky_solve(m,b,x); }
00106   static void prod(const Matrix& m, const Vector& v, Vector& x)
00107   { gsl_blas_zgemv(CblasNoTrans,gsl_complex_rect(1,0),m,v,gsl_complex_rect(0,0),x); }
00108   static void eigen_symm(const Matrix& m, gsl_vector* &eval, Matrix& evec)
00109   {
00110     gsl_eigen_hermv_workspace * w = gsl_eigen_hermv_alloc(m->size1);
00111     Matrix a = createMatrix(m->size1, m->size2);
00112     gsl_matrix_complex_memcpy(a, m);
00113     gsl_eigen_hermv(a,eval,evec,w);
00114     gsl_eigen_hermv_sort(eval, evec, GSL_EIGEN_SORT_VAL_ASC);
00115     gsl_eigen_hermv_free(w);
00116     free(a);
00117   }
00118   static void eigen_symm_gen(const Matrix& m, const Matrix& _b, gsl_vector* &eval, Matrix& evec)
00119   {
00120     gsl_eigen_genhermv_workspace * w = gsl_eigen_genhermv_alloc(m->size1);
00121     Matrix a = createMatrix(m->size1, m->size2);
00122     Matrix b = createMatrix(_b->size1, _b->size2);
00123     gsl_matrix_complex_memcpy(a, m);
00124     gsl_matrix_complex_memcpy(b, _b);
00125     gsl_eigen_genhermv(a,b,eval,evec,w);
00126     gsl_eigen_hermv_sort(eval, evec, GSL_EIGEN_SORT_VAL_ASC);
00127     gsl_eigen_genhermv_free(w);
00128     free(a);
00129   }
00130 };
00131 
00132 template<typename MatrixType>
00133 void convert(const MatrixType& m, gsl_matrix* &res)
00134 {
00135 //   if (res)
00136 //     gsl_matrix_free(res);
00137   res = gsl_matrix_alloc(m.rows(), m.cols());
00138   for (int i=0 ; i<m.rows() ; ++i)
00139     for (int j=0 ; j<m.cols(); ++j)
00140       gsl_matrix_set(res, i, j, m(i,j));
00141 }
00142 
00143 template<typename MatrixType>
00144 void convert(const gsl_matrix* m, MatrixType& res)
00145 {
00146   res.resize(int(m->size1), int(m->size2));
00147   for (int i=0 ; i<res.rows() ; ++i)
00148     for (int j=0 ; j<res.cols(); ++j)
00149       res(i,j) = gsl_matrix_get(m,i,j);
00150 }
00151 
00152 template<typename VectorType>
00153 void convert(const VectorType& m, gsl_vector* &res)
00154 {
00155   if (res) gsl_vector_free(res);
00156   res = gsl_vector_alloc(m.size());
00157   for (int i=0 ; i<m.size() ; ++i)
00158       gsl_vector_set(res, i, m[i]);
00159 }
00160 
00161 template<typename VectorType>
00162 void convert(const gsl_vector* m, VectorType& res)
00163 {
00164   res.resize (m->size);
00165   for (int i=0 ; i<res.rows() ; ++i)
00166     res[i] = gsl_vector_get(m, i);
00167 }
00168 
00169 template<typename MatrixType>
00170 void convert(const MatrixType& m, gsl_matrix_complex* &res)
00171 {
00172   res = gsl_matrix_complex_alloc(m.rows(), m.cols());
00173   for (int i=0 ; i<m.rows() ; ++i)
00174     for (int j=0 ; j<m.cols(); ++j)
00175     {
00176       gsl_matrix_complex_set(res, i, j,
00177         gsl_complex_rect(m(i,j).real(), m(i,j).imag()));
00178     }
00179 }
00180 
00181 template<typename MatrixType>
00182 void convert(const gsl_matrix_complex* m, MatrixType& res)
00183 {
00184   res.resize(int(m->size1), int(m->size2));
00185   for (int i=0 ; i<res.rows() ; ++i)
00186     for (int j=0 ; j<res.cols(); ++j)
00187       res(i,j) = typename MatrixType::Scalar(
00188         GSL_REAL(gsl_matrix_complex_get(m,i,j)),
00189         GSL_IMAG(gsl_matrix_complex_get(m,i,j)));
00190 }
00191 
00192 template<typename VectorType>
00193 void convert(const VectorType& m, gsl_vector_complex* &res)
00194 {
00195   res = gsl_vector_complex_alloc(m.size());
00196   for (int i=0 ; i<m.size() ; ++i)
00197       gsl_vector_complex_set(res, i, gsl_complex_rect(m[i].real(), m[i].imag()));
00198 }
00199 
00200 template<typename VectorType>
00201 void convert(const gsl_vector_complex* m, VectorType& res)
00202 {
00203   res.resize(m->size);
00204   for (int i=0 ; i<res.rows() ; ++i)
00205     res[i] = typename VectorType::Scalar(
00206         GSL_REAL(gsl_vector_complex_get(m, i)),
00207         GSL_IMAG(gsl_vector_complex_get(m, i)));
00208 }
00209 
00210 }
00211 
00212 #endif // EIGEN_GSL_HELPER


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