00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014 #include "lib/config.h"
00015
00016 #ifdef HAVE_LAPACK
00017 #include "lib/lapack.h"
00018 #include "lib/common.h"
00019 #include "lib/io.h"
00020
00021 using namespace shogun;
00022
00023 #if defined(HAVE_MKL) || defined(HAVE_ACML)
00024 #define DSYEV dsyev
00025 #define DGESVD dgesvd
00026 #define DPOSV dposv
00027 #define DPOTRF dpotrf
00028 #else
00029 #define DSYEV dsyev_
00030 #define DGESVD dgesvd_
00031 #define DPOSV dposv_
00032 #define DPOTRF dpotrf_
00033 #endif
00034
00035 #ifndef HAVE_ATLAS
00036 int clapack_dpotrf(const CBLAS_ORDER Order, const CBLAS_UPLO Uplo,
00037 const int N, double *A, const int LDA)
00038 {
00039 char uplo = 'U';
00040 int info = 0;
00041 if (Order==CblasRowMajor)
00042 {
00043 if (Uplo==CblasUpper)
00044 uplo='L';
00045 }
00046 else
00047 if (Uplo==CblasLower)
00048 uplo='L';
00049 #ifdef HAVE_ACML
00050 DPOTRF(uplo, N, A, LDA, &info);
00051 #else
00052 int n=N;
00053 int lda=LDA;
00054 DPOTRF(&uplo, &n, A, &lda, &info);
00055 #endif
00056 return info;
00057 }
00058 #undef DPOTRF
00059
00060
00061
00062
00063
00064
00065 int clapack_dposv(const CBLAS_ORDER Order, const CBLAS_UPLO Uplo,
00066 const int N, const int NRHS, double *A, const int lda,
00067 double *B, const int ldb)
00068 {
00069 char uplo = 'U';
00070 int info=0;
00071 if (Order==CblasRowMajor)
00072 {
00073 if (Uplo==CblasUpper)
00074 uplo='L';
00075 }
00076 else
00077 if (Uplo==CblasLower)
00078 uplo='L';
00079 #ifdef HAVE_ACML
00080 DPOSV(uplo,N,NRHS,A,lda,B,ldb,&info);
00081 #else
00082 int n=N;
00083 int nrhs=NRHS;
00084 int LDA=lda;
00085 int LDB=ldb;
00086 DPOSV(&uplo, &n, &nrhs, A, &LDA, B, &LDB, &info);
00087 #endif
00088 return info;
00089 }
00090 #undef DPOSV
00091 #endif //HAVE_ATLAS
00092
00093
00094
00095
00096
00097
00098 namespace shogun
00099 {
00100
00101
00102
00103 void wrap_dsyev(char jobz, char uplo, int n, double *a, int lda, double *w, int *info)
00104 {
00105 #ifdef HAVE_ACML
00106 DSYEV(jobz, uplo, n, a, lda, w, info);
00107 #else
00108 int lwork=-1;
00109 double work1;
00110 DSYEV(&jobz, &uplo, &n, a, &lda, w, &work1, &lwork, info);
00111 ASSERT(*info==0);
00112 ASSERT(work1>0);
00113 lwork=(int) work1;
00114 double* work=new double[lwork];
00115 DSYEV(&jobz, &uplo, &n, a, &lda, w, work, &lwork, info);
00116 delete[] work;
00117 #endif
00118 }
00119 #undef DSYEV
00120
00121 void wrap_dgesvd(char jobu, char jobvt, int m, int n, double *a, int lda, double *sing,
00122 double *u, int ldu, double *vt, int ldvt, int *info)
00123 {
00124 #ifdef HAVE_ACML
00125 DGESVD(jobu, jobvt, m, n, a, lda, sing, u, ldu, vt, ldvt, info);
00126 #else
00127 int lwork=-1;
00128 double work1;
00129 DGESVD(&jobu, &jobvt, &m, &n, a, &lda, sing, u, &ldu, vt, &ldvt, &work1, &lwork, info);
00130 ASSERT(*info==0);
00131 ASSERT(work1>0);
00132 lwork=(int) work1;
00133 double* work=new double[lwork];
00134 DGESVD(&jobu, &jobvt, &m, &n, a, &lda, sing, u, &ldu, vt, &ldvt, work, &lwork, info);
00135 delete[] work;
00136 #endif
00137 }
00138 }
00139 #undef DGESVD
00140 #endif //HAVE_LAPACK