00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #ifdef ARMA_USE_LAPACK
00019
00020
00021 namespace lapack
00022 {
00023
00024
00025
00026 extern "C"
00027 {
00028
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
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
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
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
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
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
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
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
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
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
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
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
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
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
00097
00098
00099
00100
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
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