00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015 #include <shogun/lib/config.h>
00016
00017 #ifdef HAVE_LAPACK
00018 #include <shogun/mathematics/lapack.h>
00019 #include <shogun/lib/common.h>
00020 #include <shogun/base/Parallel.h>
00021 #include <pthread.h>
00022 #include <shogun/io/SGIO.h>
00023
00024 using namespace shogun;
00025
00026 #if defined(HAVE_MKL) || defined(HAVE_ACML)
00027 #define DSYEV dsyev
00028 #define DGESVD dgesvd
00029 #define DPOSV dposv
00030 #define DPOTRF dpotrf
00031 #define DPOTRI dpotri
00032 #define DGETRI dgetri
00033 #define DGETRF dgetrf
00034 #define DGEQRF dgeqrf
00035 #define DORGQR dorgqr
00036 #define DSYEVR dsyevr
00037 #define DPOTRS dpotrs
00038 #define DGETRS dgetrs
00039 #define DSYGVX dsygvx
00040 #else
00041 #define DSYEV dsyev_
00042 #define DGESVD dgesvd_
00043 #define DPOSV dposv_
00044 #define DPOTRF dpotrf_
00045 #define DPOTRI dpotri_
00046 #define DGETRI dgetri_
00047 #define DGETRF dgetrf_
00048 #define DGEQRF dgeqrf_
00049 #define DORGQR dorgqr_
00050 #define DSYEVR dsyevr_
00051 #define DGETRS dgetrs_
00052 #define DPOTRS dpotrs_
00053 #define DSYGVX dsygvx_
00054 #endif
00055
00056 #ifndef HAVE_ATLAS
00057 int clapack_dpotrf(const CBLAS_ORDER Order, const CBLAS_UPLO Uplo,
00058 const int N, double *A, const int LDA)
00059 {
00060 char uplo = 'U';
00061 int info = 0;
00062 if (Order==CblasRowMajor)
00063 {
00064 if (Uplo==CblasUpper)
00065 uplo='L';
00066 }
00067 else if (Uplo==CblasLower)
00068 {
00069 uplo='L';
00070 }
00071 #ifdef HAVE_ACML
00072 DPOTRF(uplo, N, A, LDA, &info);
00073 #else
00074 int n=N;
00075 int lda=LDA;
00076 DPOTRF(&uplo, &n, A, &lda, &info);
00077 #endif
00078 return info;
00079 }
00080 #undef DPOTRF
00081
00082 int clapack_dpotri(const CBLAS_ORDER Order, const CBLAS_UPLO Uplo,
00083 const int N, double *A, const int LDA)
00084 {
00085 char uplo = 'U';
00086 int info = 0;
00087 if (Order==CblasRowMajor)
00088 {
00089 if (Uplo==CblasUpper)
00090 uplo='L';
00091 }
00092 else if (Uplo==CblasLower)
00093 {
00094 uplo='L';
00095 }
00096 #ifdef HAVE_ACML
00097 DPOTRI(uplo, N, A, LDA, &info);
00098 #else
00099 int n=N;
00100 int lda=LDA;
00101 DPOTRI(&uplo, &n, A, &lda, &info);
00102 #endif
00103 return info;
00104 }
00105 #undef DPOTRI
00106
00107
00108
00109
00110
00111
00112 int clapack_dposv(const CBLAS_ORDER Order, const CBLAS_UPLO Uplo,
00113 const int N, const int NRHS, double *A, const int lda,
00114 double *B, const int ldb)
00115 {
00116 char uplo = 'U';
00117 int info=0;
00118 if (Order==CblasRowMajor)
00119 {
00120 if (Uplo==CblasUpper)
00121 uplo='L';
00122 }
00123 else if (Uplo==CblasLower)
00124 {
00125 uplo='L';
00126 }
00127 #ifdef HAVE_ACML
00128 DPOSV(uplo,N,NRHS,A,lda,B,ldb,&info);
00129 #else
00130 int n=N;
00131 int nrhs=NRHS;
00132 int LDA=lda;
00133 int LDB=ldb;
00134 DPOSV(&uplo, &n, &nrhs, A, &LDA, B, &LDB, &info);
00135 #endif
00136 return info;
00137 }
00138 #undef DPOSV
00139
00140 int clapack_dgetrf(const CBLAS_ORDER Order, const int M, const int N,
00141 double *A, const int lda, int *ipiv)
00142 {
00143
00144 int info=0;
00145 #ifdef HAVE_ACML
00146 DGETRF(M,N,A,lda,ipiv,&info);
00147 #else
00148 int m=M;
00149 int n=N;
00150 int LDA=lda;
00151 DGETRF(&m,&n,A,&LDA,ipiv,&info);
00152 #endif
00153 return info;
00154 }
00155 #undef DGETRF
00156
00157
00158 int clapack_dgetri(const CBLAS_ORDER Order, const int N, double *A,
00159 const int lda, int* ipiv)
00160 {
00161 int info=0;
00162 #ifdef HAVE_ACML
00163 DGETRI(N,A,lda,ipiv,&info);
00164 #else
00165 double* work;
00166 int n=N;
00167 int LDA=lda;
00168 int lwork = -1;
00169 double work1 = 0;
00170 DGETRI(&n,A,&LDA,ipiv,&work1,&lwork,&info);
00171 lwork = (int) work1;
00172 work = SG_MALLOC(double, lwork);
00173 DGETRI(&n,A,&LDA,ipiv,work,&lwork,&info);
00174 SG_FREE(work);
00175 #endif
00176 return info;
00177 }
00178 #undef DGETRI
00179
00180
00181 int clapack_dgetrs(const CBLAS_ORDER Order, const CBLAS_TRANSPOSE Transpose,
00182 const int N, const int NRHS, double *A, const int lda,
00183 int *ipiv, double *B, const int ldb)
00184 {
00185 int info = 0;
00186 char trans = 'N';
00187 if (Transpose==CblasTrans)
00188 {
00189 trans = 'T';
00190 }
00191 #ifdef HAVE_ACML
00192 DGETRS(trans,N,NRHS,A,lda,ipiv,B,ldb,info);
00193 #else
00194 int n=N;
00195 int nrhs=NRHS;
00196 int LDA=lda;
00197 int LDB=ldb;
00198 DGETRS(&trans,&n,&nrhs,A,&LDA,ipiv,B,&LDB,&info);
00199 #endif
00200 return info;
00201 }
00202 #undef DGETRS
00203
00204
00205 int clapack_dpotrs(const CBLAS_ORDER Order, const CBLAS_UPLO Uplo,
00206 const int N, const int NRHS, double *A, const int lda,
00207 double *B, const int ldb)
00208 {
00209 int info=0;
00210 char uplo = 'U';
00211 if (Uplo==CblasLower)
00212 {
00213 uplo = 'L';
00214 }
00215 #ifdef HAVE_ACML
00216 DPOTRS(uplo,N,NRHS,A,lda,B,ldb,info);
00217 #else
00218 int n=N;
00219 int nrhs=NRHS;
00220 int LDA=lda;
00221 int LDB=ldb;
00222 DPOTRS(&uplo,&n,&nrhs,A,&LDA,B,&LDB,&info);
00223 #endif
00224 return info;
00225 }
00226 #undef DPOTRS
00227 #endif //HAVE_ATLAS
00228
00229 namespace shogun
00230 {
00231
00232 void wrap_dsyev(char jobz, char uplo, int n, double *a, int lda, double *w, int *info)
00233 {
00234 #ifdef HAVE_ACML
00235 DSYEV(jobz, uplo, n, a, lda, w, info);
00236 #else
00237 int lwork=-1;
00238 double work1;
00239 DSYEV(&jobz, &uplo, &n, a, &lda, w, &work1, &lwork, info);
00240 ASSERT(*info==0);
00241 ASSERT(work1>0);
00242 lwork=(int) work1;
00243 double* work=SG_MALLOC(double, lwork);
00244 DSYEV(&jobz, &uplo, &n, a, &lda, w, work, &lwork, info);
00245 SG_FREE(work);
00246 #endif
00247 }
00248 #undef DSYEV
00249
00250 void wrap_dgesvd(char jobu, char jobvt, int m, int n, double *a, int lda, double *sing,
00251 double *u, int ldu, double *vt, int ldvt, int *info)
00252 {
00253 #ifdef HAVE_ACML
00254 DGESVD(jobu, jobvt, m, n, a, lda, sing, u, ldu, vt, ldvt, info);
00255 #else
00256 double work1 = 0;
00257 int lworka = -1;
00258 DGESVD(&jobu, &jobvt, &m, &n, a, &lda, sing, u, &ldu, vt, &ldvt, &work1, &lworka, info);
00259 ASSERT(*info==0);
00260 lworka = (int) work1;
00261 double* worka = SG_MALLOC(double, lworka);
00262 DGESVD(&jobu, &jobvt, &m, &n, a, &lda, sing, u, &ldu, vt, &ldvt, worka, &lworka, info);
00263 SG_FREE(worka);
00264 #endif
00265 }
00266 #undef DGESVD
00267
00268 void wrap_dgeqrf(int m, int n, double *a, int lda, double *tau, int *info)
00269 {
00270 #ifdef HAVE_ACML
00271 DGEQRF(m, n, a, lda, tau, info);
00272 #else
00273 int lwork = -1;
00274 double work1 = 0;
00275 DGEQRF(&m, &n, a, &lda, tau, &work1, &lwork, info);
00276 ASSERT(*info==0);
00277 lwork = (int)work1;
00278 ASSERT(lwork>0)
00279 double* work = SG_MALLOC(double, lwork);
00280 DGEQRF(&m, &n, a, &lda, tau, work, &lwork, info);
00281 SG_FREE(work);
00282 #endif
00283 }
00284 #undef DGEQRF
00285
00286 void wrap_dorgqr(int m, int n, int k, double *a, int lda, double *tau, int *info)
00287 {
00288 #ifdef HAVE_ACML
00289 DORGQR(m, n, k, a, lda, tau, info);
00290 #else
00291 int lwork = -1;
00292 double work1 = 0;
00293 DORGQR(&m, &n, &k, a, &lda, tau, &work1, &lwork, info);
00294 ASSERT(*info==0);
00295 lwork = (int)work1;
00296 ASSERT(lwork>0);
00297 double* work = SG_MALLOC(double, lwork);
00298 DORGQR(&m, &n, &k, a, &lda, tau, work, &lwork, info);
00299 SG_FREE(work);
00300 #endif
00301 }
00302 #undef DORGQR
00303
00304 void wrap_dsyevr(char jobz, char uplo, int n, double *a, int lda, int il, int iu,
00305 double *eigenvalues, double *eigenvectors, int *info)
00306 {
00307 int m;
00308 double vl,vu;
00309 double abstol = 0.0;
00310 char I = 'I';
00311 int* isuppz = SG_MALLOC(int, n);
00312 #ifdef HAVE_ACML
00313 DSYEVR(jobz,I,uplo,n,a,lda,vl,vu,il,iu,abstol,m,
00314 eigenvalues,eigenvectors,n,isuppz,info);
00315 #else
00316 int lwork = -1;
00317 int liwork = -1;
00318 double work1 = 0;
00319 int work2 = 0;
00320 DSYEVR(&jobz,&I,&uplo,&n,a,&lda,&vl,&vu,&il,&iu,&abstol,
00321 &m,eigenvalues,eigenvectors,&n,isuppz,
00322 &work1,&lwork,&work2,&liwork,info);
00323 ASSERT(*info==0);
00324 lwork = (int)work1;
00325 liwork = work2;
00326 double* work = SG_MALLOC(double, lwork);
00327 int* iwork = SG_MALLOC(int, liwork);
00328 DSYEVR(&jobz,&I,&uplo,&n,a,&lda,&vl,&vu,&il,&iu,&abstol,
00329 &m,eigenvalues,eigenvectors,&n,isuppz,
00330 work,&lwork,iwork,&liwork,info);
00331 ASSERT(*info==0);
00332 SG_FREE(work);
00333 SG_FREE(iwork);
00334 SG_FREE(isuppz);
00335 #endif
00336 }
00337 #undef DSYEVR
00338
00339 void wrap_dsygvx(int itype, char jobz, char uplo, int n, double *a, int lda, double *b,
00340 int ldb, int il, int iu, double *eigenvalues, double *eigenvectors, int *info)
00341 {
00342 int m;
00343 double abstol = 0.0;
00344 double vl,vu;
00345 int* ifail = SG_MALLOC(int, n);
00346 char I = 'I';
00347 #ifdef HAVE_ACML
00348 DSYGVX(itype,jobz,I,uplo,n,a,lda,b,ldb,vl,vu,
00349 il,iu,abstol,m,eigenvalues,
00350 eigenvectors,n,ifail,info);
00351 #else
00352 int lwork = -1;
00353 double work1 = 0;
00354 int* iwork = SG_MALLOC(int, 5*n);
00355 DSYGVX(&itype,&jobz,&I,&uplo,&n,a,&lda,b,&ldb,&vl,&vu,
00356 &il,&iu,&abstol,&m,eigenvalues,eigenvectors,
00357 &n,&work1,&lwork,iwork,ifail,info);
00358 lwork = (int)work1;
00359 double* work = SG_MALLOC(double, lwork);
00360 DSYGVX(&itype,&jobz,&I,&uplo,&n,a,&lda,b,&ldb,&vl,&vu,
00361 &il,&iu,&abstol,&m,eigenvalues,eigenvectors,
00362 &n,work,&lwork,iwork,ifail,info);
00363 SG_FREE(work);
00364 SG_FREE(iwork);
00365 SG_FREE(ifail);
00366 #endif
00367 }
00368 #undef DSYGVX
00369 }
00370 #endif //HAVE_LAPACK