lapack_proto.hpp

Go to the documentation of this file.
00001 // Copyright (C) 2010 NICTA and the authors listed below
00002 // http://nicta.com.au
00003 // 
00004 // Authors:
00005 // - Conrad Sanderson (conradsand at ieee dot org)
00006 // - Edmund Highcock (edmund dot highcock at merton dot ox dot ac dot uk)
00007 // 
00008 // This file is part of the Armadillo C++ library.
00009 // It is provided without any warranty of fitness
00010 // for any purpose. You can redistribute this file
00011 // and/or modify it under the terms of the GNU
00012 // Lesser General Public License (LGPL) as published
00013 // by the Free Software Foundation, either version 3
00014 // of the License or (at your option) any later version.
00015 // (see http://www.opensource.org/licenses for more info)
00016 
00017 
00018 #ifdef ARMA_USE_LAPACK
00019 
00020 //! \namespace lapack namespace for LAPACK functions
00021 namespace lapack
00022   {
00023   //! \addtogroup LAPACK
00024   //! @{
00025   
00026   extern "C"
00027     {
00028     // LU factorisation
00029     void sgetrf_(int* m, int* n,  float* a, int* lda, int* ipiv, int* info);
00030     void dgetrf_(int* m, int* n, double* a, int* lda, int* ipiv, int* info);
00031     void cgetrf_(int* m, int* n,   void* a, int* lda, int* ipiv, int* info);
00032     void zgetrf_(int* m, int* n,   void* a, int* lda, int* ipiv, int* info);
00033     
00034     // matrix inversion
00035     void sgetri_(int* n,  float* a, int* lda, int* ipiv,  float* work, int* lwork, int* info);
00036     void dgetri_(int* n, double* a, int* lda, int* ipiv, double* work, int* lwork, int* info);
00037     void cgetri_(int* n,  void*  a, int* lda, int* ipiv,   void* work, int* lwork, int* info);
00038     void zgetri_(int* n,  void*  a, int* lda, int* ipiv,   void* work, int* lwork, int* info);
00039     
00040     // eigenvector decomposition of symmetric real matrices
00041     void ssyev_(char* jobz, char* uplo, int* n,  float* a, int* lda,  float* w,  float* work, int* lwork, int* info);
00042     void dsyev_(char* jobz, char* uplo, int* n, double* a, int* lda, double* w, double* work, int* lwork, int* info);
00043 
00044     // eigenvector decomposition of hermitian matrices (complex)
00045     void cheev_(char* jobz, char* uplo, int* n,   void* a, int* lda,  float* w,   void* work, int* lwork,  float* rwork, int* info);
00046     void zheev_(char* jobz, char* uplo, int* n,   void* a, int* lda, double* w,   void* work, int* lwork, double* rwork, int* info);
00047 
00048     // eigenvector decomposition of general real matrices
00049     void sgeev_(char* jobvl, char* jobvr, int* n,  float* a, int* lda,  float* wr,  float* wi,  float* vl, int* ldvl,  float* vr, int* ldvr,  float* work, int* lwork, int* info);
00050     void dgeev_(char* jobvl, char* jobvr, int* n, double* a, int* lda, double* wr, double* wi, double* vl, int* ldvl, double* vr, int* ldvr, double* work, int* lwork, int* info);
00051 
00052     // eigenvector decomposition of general complex matrices
00053     void cgeev_(char* jobvr, char* jobvl, int* n, void* a, int* lda, void* w, void* vl, int* ldvl, void* vr, int* ldvr, void* work, int* lwork,  float* rwork, int* info);
00054     void zgeev_(char* jobvl, char* jobvr, int* n, void* a, int* lda, void* w, void* vl, int *ldvl, void* vr, int *ldvr, void* work, int* lwork, double* rwork, int* info);
00055 
00056     // Cholesky decomposition
00057     void spotrf_(char* uplo, int* n,  float* a, int* lda, int* info);
00058     void dpotrf_(char* uplo, int* n, double* a, int* lda, int* info);
00059     void cpotrf_(char* uplo, int* n,   void* a, int* lda, int* info);
00060     void zpotrf_(char* uplo, int* n,   void* a, int* lda, int* info);
00061 
00062     // QR decomposition
00063     void sgeqrf_(int* m, int* n,  float* a, int* lda,  float* tau,  float* work, int* lwork, int* info);
00064     void dgeqrf_(int* m, int* n, double* a, int* lda, double* tau, double* work, int* lwork, int* info);
00065     void cgeqrf_(int* m, int* n,   void* a, int* lda,   void* tau,   void* work, int* lwork, int* info);
00066     void zgeqrf_(int* m, int* n,   void* a, int* lda,   void* tau,   void* work, int* lwork, int* info);
00067 
00068     // Q matrix calculation from QR decomposition (real matrices)
00069     void sorgqr_(int* m, int* n, int* k,  float* a, int* lda,  float* tau,  float* work, int* lwork, int* info);
00070     void dorgqr_(int* m, int* n, int* k, double* a, int* lda, double* tau, double* work, int* lwork, int* info);
00071     
00072     // Q matrix calculation from QR decomposition (complex matrices)
00073     void cungqr_(int* m, int* n, int* k,   void* a, int* lda,   void* tau,   void* work, int* lwork, int* info);
00074     void zungqr_(int* m, int* n, int* k,   void* a, int* lda,   void* tau,   void* work, int* lwork, int* info);
00075 
00076     // SVD (real matrices)
00077     void sgesvd_(char* jobu, char* jobvt, int* m, int* n, float*  a, int* lda, float*  s, float*  u, int* ldu, float*  vt, int* ldvt, float*  work, int* lwork, int* info);
00078     void dgesvd_(char* jobu, char* jobvt, int* m, int* n, double* a, int* lda, double* s, double* u, int* ldu, double* vt, int* ldvt, double* work, int* lwork, int* info);
00079     
00080     // SVD (complex matrices)
00081     void cgesvd_(char* jobu, char* jobvt, int* m, int* n, void*   a, int* lda, float*  s, void*   u, int* ldu, void*   vt, int* ldvt, void*   work, int* lwork, float*  rwork, int* info);
00082     void zgesvd_(char* jobu, char* jobvt, int* m, int* n, void*   a, int* lda, double* s, void*   u, int* ldu, void*   vt, int* ldvt, void*   work, int* lwork, double* rwork, int* info);
00083 
00084     // solve system of linear equations, using LU decomposition
00085     void sgesv_(int* n, int* nrhs, float*  a, int* lda, int* ipiv, float*  b, int* ldb, int* info);
00086     void dgesv_(int* n, int* nrhs, double* a, int* lda, int* ipiv, double* b, int* ldb, int* info);
00087     void cgesv_(int* n, int* nrhs, void*   a, int* lda, int* ipiv, void*   b, int* ldb, int* info);
00088     void zgesv_(int* n, int* nrhs, void*   a, int* lda, int* ipiv, void*   b, int* ldb, int* info);
00089 
00090     // solve over/underdetermined system of linear equations
00091     void sgels_(char* trans, int* m, int* n, int* nrhs, float*  a, int* lda, float*  b, int* ldb, float*  work, int* lwork, int* info);
00092     void dgels_(char* trans, int* m, int* n, int* nrhs, double* a, int* lda, double* b, int* ldb, double* work, int* lwork, int* info);
00093     void cgels_(char* trans, int* m, int* n, int* nrhs, void*   a, int *lda, void*   b, int* ldb, void*   work, int* lwork, int* info);
00094     void zgels_(char* trans, int* m, int* n, int* nrhs, void*   a, int *lda, void*   b, int* ldb, void*   work, int* lwork, int* info);
00095 
00096     // void dgeqp3_(int* m, int* n, double* a, int* lda, int* jpvt, double* tau, double* work, int* lwork, int* info);
00097     // void dormqr_(char* side, char* trans, int* m, int* n, int* k, double* a, int* lda, double* tau, double* c, int* ldc, double* work, int* lwork, int* info);
00098     // void  dposv_(char* uplo, int* n, int* nrhs, double* a, int* lda, double* b, int* ldb, int* info);
00099     // void dtrtrs_(char* uplo, char* trans, char* diag, int* n, int* nrhs, double* a, int* lda, double* b, int* ldb, int* info);
00100     // void  dgees_(char* jobvs, char* sort, int* select, int* n, double* a, int* lda, int* sdim, double* wr, double* wi, double* vs, int* ldvs, double* work, int* lwork, int* bwork, int* info);
00101     
00102     }
00103 
00104   template<typename eT>
00105   inline
00106   void
00107   getrf_(int* m, int* n, eT* a, int* lda, int* ipiv, int* info)
00108     {
00109     arma_type_check<is_supported_blas_type<eT>::value == false>::apply();
00110     
00111     if(is_float<eT>::value == true)
00112       {
00113       typedef float T;
00114       sgetrf_(m, n, (T*)a, lda, ipiv, info);
00115       }
00116     else
00117     if(is_double<eT>::value == true)
00118       {
00119       typedef double T;
00120       dgetrf_(m, n, (T*)a, lda, ipiv, info);
00121       }
00122     else
00123     if(is_supported_complex_float<eT>::value == true)
00124       {
00125       typedef std::complex<float> T;
00126       cgetrf_(m, n, (T*)a, lda, ipiv, info);
00127       }
00128     else
00129     if(is_supported_complex_double<eT>::value == true)
00130       {
00131       typedef std::complex<double> T;
00132       zgetrf_(m, n, (T*)a, lda, ipiv, info);
00133       }
00134     }
00135     
00136     
00137     
00138   template<typename eT>
00139   inline
00140   void
00141   getri_(int* n,  eT* a, int* lda, int* ipiv, eT* work, int* lwork, int* info)
00142     {
00143     arma_type_check<is_supported_blas_type<eT>::value == false>::apply();
00144     
00145     if(is_float<eT>::value == true)
00146       {
00147       typedef float T;
00148       sgetri_(n, (T*)a, lda, ipiv, (T*)work, lwork, info);
00149       }
00150     else
00151     if(is_double<eT>::value == true)
00152       {
00153       typedef double T;
00154       dgetri_(n, (T*)a, lda, ipiv, (T*)work, lwork, info);
00155       }
00156     else
00157     if(is_supported_complex_float<eT>::value == true)
00158       {
00159       typedef std::complex<float> T;
00160       cgetri_(n, (T*)a, lda, ipiv, (T*)work, lwork, info);
00161       }
00162     else
00163     if(is_supported_complex_double<eT>::value == true)
00164       {
00165       typedef std::complex<double> T;
00166       zgetri_(n, (T*)a, lda, ipiv, (T*)work, lwork, info);
00167       }
00168     }
00169   
00170   
00171   
00172   template<typename eT>
00173   inline
00174   void
00175   syev_(char* jobz, char* uplo, int* n, eT* a, int* lda, eT* w,  eT* work, int* lwork, int* info)
00176     {
00177     arma_type_check<is_supported_blas_type<eT>::value == false>::apply();
00178     
00179     if(is_float<eT>::value == true)
00180       {
00181       typedef float T;
00182       ssyev_(jobz, uplo, n, (T*)a, lda, (T*)w, (T*)work, lwork, info);
00183       }
00184     else
00185     if(is_double<eT>::value == true)
00186       {
00187       typedef double T;
00188       dsyev_(jobz, uplo, n, (T*)a, lda, (T*)w, (T*)work, lwork, info);
00189       }
00190     }
00191   
00192   
00193   
00194   template<typename eT>
00195   inline
00196   void
00197   heev_
00198     (
00199     char* jobz, char* uplo, int* n,
00200     eT* a, int* lda, typename eT::value_type* w,
00201     eT* work, int* lwork, typename eT::value_type* rwork,
00202     int* info
00203     )
00204     {
00205     arma_type_check<is_supported_blas_type<eT>::value == false>::apply();
00206     
00207     if(is_supported_complex_float<eT>::value == true)
00208       {
00209       typedef float T;
00210       typedef typename std::complex<T> cx_T;
00211       cheev_(jobz, uplo, n, (cx_T*)a, lda, (T*)w, (cx_T*)work, lwork, (T*)rwork, info);
00212       }
00213     else
00214     if(is_supported_complex_double<eT>::value == true)
00215       {
00216       typedef double T;
00217       typedef typename std::complex<T> cx_T;
00218       zheev_(jobz, uplo, n, (cx_T*)a, lda, (T*)w, (cx_T*)work, lwork, (T*)rwork, info);
00219       }
00220     }
00221   
00222    
00223   template<typename eT>
00224   inline
00225   void
00226   geev_
00227     (
00228     char* jobvl, char* jobvr, int* n, 
00229     eT* a, int* lda, eT* wr, eT* wi, eT* vl, 
00230     int* ldvl, eT* vr, int* ldvr, 
00231     eT* work, int* lwork,
00232     int* info
00233     )
00234     {
00235     arma_type_check<is_supported_blas_type<eT>::value == false>::apply();
00236 
00237     if(is_float<eT>::value == true)
00238       {
00239       typedef float T;
00240       sgeev_(jobvl, jobvr, n,  (T*)a, lda, (T*)wr, (T*)wi, (T*)vl, ldvl, (T*)vr, ldvr, (T*)work, lwork, info);
00241       }
00242     else
00243     if(is_double<eT>::value == true)
00244       {
00245       typedef double T;
00246       dgeev_(jobvl, jobvr, n,  (T*)a, lda, (T*)wr, (T*)wi, (T*)vl, ldvl, (T*)vr, ldvr, (T*)work, lwork, info);
00247       }
00248     }
00249 
00250 
00251   template<typename eT>
00252   inline
00253   void
00254   cx_geev_
00255     (
00256     char* jobvl, char* jobvr, int* n, 
00257     eT* a, int* lda, eT* w, 
00258     eT* vl, int* ldvl, 
00259     eT* vr, int* ldvr, 
00260     eT* work, int* lwork, typename eT::value_type* rwork, 
00261     int* info
00262     )
00263     {
00264     arma_type_check<is_supported_blas_type<eT>::value == false>::apply();
00265     
00266     if(is_supported_complex_float<eT>::value == true)
00267       {
00268       typedef float T;
00269       typedef typename std::complex<T> cx_T;
00270       cgeev_(jobvl, jobvr, n, (cx_T*)a, lda, (cx_T*)w, (cx_T*)vl, ldvl, (cx_T*)vr, ldvr, (cx_T*)work, lwork, (T*)rwork, info);
00271       }
00272     else
00273     if(is_supported_complex_double<eT>::value == true)
00274       {
00275       typedef double T;
00276       typedef typename std::complex<T> cx_T;
00277       zgeev_(jobvl, jobvr, n, (cx_T*)a, lda, (cx_T*)w, (cx_T*)vl, ldvl, (cx_T*)vr, ldvr, (cx_T*)work, lwork, (T*)rwork, info);
00278       }
00279     }
00280   
00281 
00282 
00283   
00284   template<typename eT>
00285   inline
00286   void
00287   potrf_(char* uplo, int* n, eT* a, int* lda, int* info)
00288     {
00289     arma_type_check<is_supported_blas_type<eT>::value == false>::apply();
00290     
00291     if(is_float<eT>::value == true)
00292       {
00293       typedef float T;
00294       spotrf_(uplo, n, (T*)a, lda, info);
00295       }
00296     else
00297     if(is_double<eT>::value == true)
00298       {
00299       typedef double T;
00300       dpotrf_(uplo, n, (T*)a, lda, info);
00301       }
00302     else
00303     if(is_supported_complex_float<eT>::value == true)
00304       {
00305       typedef std::complex<float> T;
00306       cpotrf_(uplo, n, (T*)a, lda, info);
00307       }
00308     else
00309     if(is_supported_complex_double<eT>::value == true)
00310       {
00311       typedef std::complex<double> T;
00312       zpotrf_(uplo, n, (T*)a, lda, info);
00313       }
00314     
00315     }
00316   
00317   
00318   
00319   template<typename eT>
00320   inline
00321   void
00322   geqrf_(int* m, int* n, eT* a, int* lda, eT* tau, eT* work, int* lwork, int* info)
00323     {
00324     arma_type_check<is_supported_blas_type<eT>::value == false>::apply();
00325     
00326     if(is_float<eT>::value == true)
00327       {
00328       typedef float T;
00329       sgeqrf_(m, n, (T*)a, lda, (T*)tau, (T*)work, lwork, info);
00330       }
00331     else
00332     if(is_double<eT>::value == true)
00333       {
00334       typedef double T;
00335       dgeqrf_(m, n, (T*)a, lda, (T*)tau, (T*)work, lwork, info);
00336       }
00337     else
00338     if(is_supported_complex_float<eT>::value == true)
00339       {
00340       typedef std::complex<float> T;
00341       cgeqrf_(m, n, (T*)a, lda, (T*)tau, (T*)work, lwork, info);
00342       }
00343     else
00344     if(is_supported_complex_double<eT>::value == true)
00345       {
00346       typedef std::complex<double> T;
00347       zgeqrf_(m, n, (T*)a, lda, (T*)tau, (T*)work, lwork, info);
00348       }
00349     
00350     }
00351   
00352   
00353   
00354   template<typename eT>
00355   inline
00356   void
00357   orgqr_(int* m, int* n, int* k, eT* a, int* lda, eT* tau, eT* work, int* lwork, int* info)
00358     {
00359     arma_type_check<is_supported_blas_type<eT>::value == false>::apply();
00360     
00361     if(is_float<eT>::value == true)
00362       {
00363       typedef float T;
00364       sorgqr_(m, n, k, (T*)a, lda, (T*)tau, (T*)work, lwork, info);
00365       }
00366     else
00367     if(is_double<eT>::value == true)
00368       {
00369       typedef double T;
00370       dorgqr_(m, n, k, (T*)a, lda, (T*)tau, (T*)work, lwork, info);
00371       }
00372     }
00373 
00374 
00375   
00376   template<typename eT>
00377   inline
00378   void
00379   ungqr_(int* m, int* n, int* k, eT* a, int* lda, eT* tau, eT* work, int* lwork, int* info)
00380     {
00381     arma_type_check<is_supported_blas_type<eT>::value == false>::apply();
00382     
00383     if(is_supported_complex_float<eT>::value == true)
00384       {
00385       typedef float T;
00386       cungqr_(m, n, k, (T*)a, lda, (T*)tau, (T*)work, lwork, info);
00387       }
00388     else
00389     if(is_supported_complex_double<eT>::value == true)
00390       {
00391       typedef double T;
00392       zungqr_(m, n, k, (T*)a, lda, (T*)tau, (T*)work, lwork, info);
00393       }
00394     }
00395   
00396   
00397   template<typename eT>
00398   inline
00399   void
00400   gesvd_
00401     (
00402     char* jobu, char* jobvt, int* m, int* n, eT* a, int* lda,
00403     eT* s, eT* u, int* ldu, eT* vt, int* ldvt,
00404     eT* work, int* lwork, int* info
00405     )
00406     {
00407     arma_type_check<is_supported_blas_type<eT>::value == false>::apply();
00408     
00409     if(is_float<eT>::value == true)
00410       {
00411       typedef float T;
00412       sgesvd_(jobu, jobvt, m, n, (T*)a, lda, (T*)s, (T*)u, ldu, (T*)vt, ldvt, (T*)work, lwork, info);
00413       }
00414     else
00415     if(is_double<eT>::value == true)
00416       {
00417       typedef double T;
00418       dgesvd_(jobu, jobvt, m, n, (T*)a, lda, (T*)s, (T*)u, ldu, (T*)vt, ldvt, (T*)work, lwork, info);
00419       }
00420     }
00421   
00422   
00423   
00424   template<typename T>
00425   inline
00426   void
00427   cx_gesvd_
00428     (
00429     char* jobu, char* jobvt, int* m, int* n, std::complex<T>* a, int* lda,
00430     T* s, std::complex<T>* u, int* ldu, std::complex<T>* vt, int* ldvt, 
00431     std::complex<T>* work, int* lwork, T* rwork, int* info
00432     )
00433     {
00434     arma_type_check<is_supported_blas_type<T>::value == false>::apply();
00435     arma_type_check<is_supported_blas_type< std::complex<T> >::value == false>::apply();
00436     
00437     if(is_float<T>::value == true)
00438       {
00439       typedef float bT;
00440       cgesvd_
00441         (
00442         jobu, jobvt, m, n, (std::complex<bT>*)a, lda,
00443         (bT*)s, (std::complex<bT>*)u, ldu, (std::complex<bT>*)vt, ldvt,
00444         (std::complex<bT>*)work, lwork, (bT*)rwork, info
00445         );
00446       }
00447     else
00448     if(is_double<T>::value == true)
00449       {
00450       typedef double bT;
00451       zgesvd_
00452         (
00453         jobu, jobvt, m, n, (std::complex<bT>*)a, lda,
00454         (bT*)s, (std::complex<bT>*)u, ldu, (std::complex<bT>*)vt, ldvt,
00455         (std::complex<bT>*)work, lwork, (bT*)rwork, info
00456         );
00457       }
00458     }
00459   
00460   
00461   
00462   template<typename eT>
00463   inline
00464   void
00465   gesv_(int* n, int* nrhs, eT* a, int* lda, int* ipiv, eT* b, int* ldb, int* info)
00466     {
00467     arma_type_check<is_supported_blas_type<eT>::value == false>::apply();
00468     
00469     if(is_float<eT>::value == true)
00470       {
00471       typedef float T;
00472       sgesv_(n, nrhs, (T*)a, lda, ipiv, (T*)b, ldb, info);
00473       }
00474     else
00475     if(is_double<eT>::value == true)
00476       {
00477       typedef double T;
00478       dgesv_(n, nrhs, (T*)a, lda, ipiv, (T*)b, ldb, info);
00479       }
00480     else
00481     if(is_supported_complex_float<eT>::value == true)
00482       {
00483       typedef std::complex<float> T;
00484       cgesv_(n, nrhs, (T*)a, lda, ipiv, (T*)b, ldb, info);
00485       }
00486     else
00487     if(is_supported_complex_double<eT>::value == true)
00488       {
00489       typedef std::complex<double> T;
00490       zgesv_(n, nrhs, (T*)a, lda, ipiv, (T*)b, ldb, info);
00491       }
00492     }
00493   
00494   
00495   
00496   template<typename eT>
00497   inline
00498   void
00499 //sgels_(char* trans, int* m, int* n, int* nrhs, float*  a, int* lda, float*  b, int* ldb, float*  work, int* lwork, int* info);
00500   gels_(char* trans, int* m, int* n, int* nrhs, eT* a, int* lda, eT* b, int* ldb, eT* work, int* lwork, int* info)
00501     {
00502     arma_type_check<is_supported_blas_type<eT>::value == false>::apply();
00503     
00504     if(is_float<eT>::value == true)
00505       {
00506       typedef float T;
00507       sgels_(trans, m, n, nrhs, (T*)a, lda, (T*)b, ldb, (T*)work, lwork, info);
00508       }
00509     else
00510     if(is_double<eT>::value == true)
00511       {
00512       typedef double T;
00513       dgels_(trans, m, n, nrhs, (T*)a, lda, (T*)b, ldb, (T*)work, lwork, info);
00514       }
00515     else
00516     if(is_supported_complex_float<eT>::value == true)
00517       {
00518       typedef std::complex<float> T;
00519       cgels_(trans, m, n, nrhs, (T*)a, lda, (T*)b, ldb, (T*)work, lwork, info);
00520       }
00521     else
00522     if(is_supported_complex_double<eT>::value == true)
00523       {
00524       typedef std::complex<double> T;
00525       zgels_(trans, m, n, nrhs, (T*)a, lda, (T*)b, ldb, (T*)work, lwork, info);
00526       }
00527     }
00528   
00529   
00530   
00531   //! @}
00532   }
00533 
00534 #endif
00535