lapack_proto.hpp

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