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