SHOGUN
v3.2.0
|
00001 /* 00002 * This program is free software; you can redistribute it and/or modify 00003 * it under the terms of the GNU General Public License as published by 00004 * the Free Software Foundation; either version 3 of the License, or 00005 * (at your option) any later version. 00006 * 00007 * Written (W) 1999-2009 Soeren Sonnenburg 00008 * Written (W) 1999-2008 Gunnar Raetsch 00009 * Written (W) 2006-2007 Mikio L. Braun 00010 * Written (W) 2008 Jochen Garcke 00011 * Written (W) 2011 Sergey Lisitsyn 00012 * Copyright (C) 1999-2009 Fraunhofer Institute FIRST and Max-Planck-Society 00013 */ 00014 00015 #include <shogun/lib/config.h> 00016 00017 #ifdef HAVE_LAPACK 00018 #include <shogun/lib/common.h> 00019 #include <shogun/mathematics/lapack.h> 00020 #include <shogun/base/Parallel.h> 00021 #include <shogun/io/SGIO.h> 00022 00023 #include <pthread.h> 00024 00025 using namespace shogun; 00026 00027 #if defined(HAVE_MKL) || defined(HAVE_ACML) 00028 #define DSYEV dsyev 00029 #define DGESVD dgesvd 00030 #define DPOSV dposv 00031 #define DPOTRF dpotrf 00032 #define DPOTRI dpotri 00033 #define DGETRI dgetri 00034 #define DGETRF dgetrf 00035 #define DGEQRF dgeqrf 00036 #define DORGQR dorgqr 00037 #define DSYEVR dsyevr 00038 #define DPOTRS dpotrs 00039 #define DGETRS dgetrs 00040 #define DSYGVX dsygvx 00041 #define DSTEMR dstemr 00042 #else 00043 #define DSYEV dsyev_ 00044 #define DGESVD dgesvd_ 00045 #define DPOSV dposv_ 00046 #define DPOTRF dpotrf_ 00047 #define DPOTRI dpotri_ 00048 #define DGETRI dgetri_ 00049 #define DGETRF dgetrf_ 00050 #define DGEQRF dgeqrf_ 00051 #define DORGQR dorgqr_ 00052 #define DSYEVR dsyevr_ 00053 #define DGETRS dgetrs_ 00054 #define DPOTRS dpotrs_ 00055 #define DSYGVX dsygvx_ 00056 #define DSTEMR dstemr_ 00057 #endif 00058 00059 #ifndef HAVE_ATLAS 00060 int clapack_dpotrf(const CBLAS_ORDER Order, const CBLAS_UPLO Uplo, 00061 const int N, double *A, const int LDA) 00062 { 00063 char uplo = 'U'; 00064 int info = 0; 00065 if (Order==CblasRowMajor) 00066 {//A is symmetric, we switch Uplo to get result for CblasRowMajor 00067 if (Uplo==CblasUpper) 00068 uplo='L'; 00069 } 00070 else if (Uplo==CblasLower) 00071 { 00072 uplo='L'; 00073 } 00074 #ifdef HAVE_ACML 00075 DPOTRF(uplo, N, A, LDA, &info); 00076 #else 00077 int n=N; 00078 int lda=LDA; 00079 DPOTRF(&uplo, &n, A, &lda, &info); 00080 #endif 00081 return info; 00082 } 00083 #undef DPOTRF 00084 00085 int clapack_dpotri(const CBLAS_ORDER Order, const CBLAS_UPLO Uplo, 00086 const int N, double *A, const int LDA) 00087 { 00088 char uplo = 'U'; 00089 int info = 0; 00090 if (Order==CblasRowMajor) 00091 {//A is symmetric, we switch Uplo to get result for CblasRowMajor 00092 if (Uplo==CblasUpper) 00093 uplo='L'; 00094 } 00095 else if (Uplo==CblasLower) 00096 { 00097 uplo='L'; 00098 } 00099 #ifdef HAVE_ACML 00100 DPOTRI(uplo, N, A, LDA, &info); 00101 #else 00102 int n=N; 00103 int lda=LDA; 00104 DPOTRI(&uplo, &n, A, &lda, &info); 00105 #endif 00106 return info; 00107 } 00108 #undef DPOTRI 00109 00110 /* DPOSV computes the solution to a real system of linear equations 00111 * A * X = B, 00112 * where A is an N-by-N symmetric positive definite matrix and X and B 00113 * are N-by-NRHS matrices 00114 */ 00115 int clapack_dposv(const CBLAS_ORDER Order, const CBLAS_UPLO Uplo, 00116 const int N, const int NRHS, double *A, const int lda, 00117 double *B, const int ldb) 00118 { 00119 char uplo = 'U'; 00120 int info=0; 00121 if (Order==CblasRowMajor) 00122 {//A is symmetric, we switch Uplo to achieve CblasColMajor 00123 if (Uplo==CblasUpper) 00124 uplo='L'; 00125 } 00126 else if (Uplo==CblasLower) 00127 { 00128 uplo='L'; 00129 } 00130 #ifdef HAVE_ACML 00131 DPOSV(uplo,N,NRHS,A,lda,B,ldb,&info); 00132 #else 00133 int n=N; 00134 int nrhs=NRHS; 00135 int LDA=lda; 00136 int LDB=ldb; 00137 DPOSV(&uplo, &n, &nrhs, A, &LDA, B, &LDB, &info); 00138 #endif 00139 return info; 00140 } 00141 #undef DPOSV 00142 00143 int clapack_dgetrf(const CBLAS_ORDER Order, const int M, const int N, 00144 double *A, const int lda, int *ipiv) 00145 { 00146 // no rowmajor? 00147 int info=0; 00148 #ifdef HAVE_ACML 00149 DGETRF(M,N,A,lda,ipiv,&info); 00150 #else 00151 int m=M; 00152 int n=N; 00153 int LDA=lda; 00154 DGETRF(&m,&n,A,&LDA,ipiv,&info); 00155 #endif 00156 return info; 00157 } 00158 #undef DGETRF 00159 00160 // order not supported (yet?) 00161 int clapack_dgetri(const CBLAS_ORDER Order, const int N, double *A, 00162 const int lda, int* ipiv) 00163 { 00164 int info=0; 00165 #ifdef HAVE_ACML 00166 DGETRI(N,A,lda,ipiv,&info); 00167 #else 00168 double* work; 00169 int n=N; 00170 int LDA=lda; 00171 int lwork = -1; 00172 double work1 = 0; 00173 DGETRI(&n,A,&LDA,ipiv,&work1,&lwork,&info); 00174 lwork = (int) work1; 00175 work = SG_MALLOC(double, lwork); 00176 DGETRI(&n,A,&LDA,ipiv,work,&lwork,&info); 00177 SG_FREE(work); 00178 #endif 00179 return info; 00180 } 00181 #undef DGETRI 00182 00183 // order not supported (yet?) 00184 int clapack_dgetrs(const CBLAS_ORDER Order, const CBLAS_TRANSPOSE Transpose, 00185 const int N, const int NRHS, double *A, const int lda, 00186 int *ipiv, double *B, const int ldb) 00187 { 00188 int info = 0; 00189 char trans = 'N'; 00190 if (Transpose==CblasTrans) 00191 { 00192 trans = 'T'; 00193 } 00194 #ifdef HAVE_ACML 00195 DGETRS(trans,N,NRHS,A,lda,ipiv,B,ldb,info); 00196 #else 00197 int n=N; 00198 int nrhs=NRHS; 00199 int LDA=lda; 00200 int LDB=ldb; 00201 DGETRS(&trans,&n,&nrhs,A,&LDA,ipiv,B,&LDB,&info); 00202 #endif 00203 return info; 00204 } 00205 #undef DGETRS 00206 00207 // order not supported (yet?) 00208 int clapack_dpotrs(const CBLAS_ORDER Order, const CBLAS_UPLO Uplo, 00209 const int N, const int NRHS, double *A, const int lda, 00210 double *B, const int ldb) 00211 { 00212 int info=0; 00213 char uplo = 'U'; 00214 if (Uplo==CblasLower) 00215 { 00216 uplo = 'L'; 00217 } 00218 #ifdef HAVE_ACML 00219 DPOTRS(uplo,N,NRHS,A,lda,B,ldb,info); 00220 #else 00221 int n=N; 00222 int nrhs=NRHS; 00223 int LDA=lda; 00224 int LDB=ldb; 00225 DPOTRS(&uplo,&n,&nrhs,A,&LDA,B,&LDB,&info); 00226 #endif 00227 return info; 00228 } 00229 #undef DPOTRS 00230 #endif //HAVE_ATLAS 00231 00232 namespace shogun 00233 { 00234 00235 void wrap_dsyev(char jobz, char uplo, int n, double *a, int lda, double *w, int *info) 00236 { 00237 #ifdef HAVE_ACML 00238 DSYEV(jobz, uplo, n, a, lda, w, info); 00239 #else 00240 int lwork=-1; 00241 double work1; 00242 DSYEV(&jobz, &uplo, &n, a, &lda, w, &work1, &lwork, info); 00243 ASSERT(*info==0) 00244 ASSERT(work1>0) 00245 lwork=(int) work1; 00246 double* work=SG_MALLOC(double, lwork); 00247 DSYEV(&jobz, &uplo, &n, a, &lda, w, work, &lwork, info); 00248 SG_FREE(work); 00249 #endif 00250 } 00251 #undef DSYEV 00252 00253 void wrap_dgesvd(char jobu, char jobvt, int m, int n, double *a, int lda, double *sing, 00254 double *u, int ldu, double *vt, int ldvt, int *info) 00255 { 00256 #ifdef HAVE_ACML 00257 DGESVD(jobu, jobvt, m, n, a, lda, sing, u, ldu, vt, ldvt, info); 00258 #else 00259 double work1 = 0; 00260 int lworka = -1; 00261 DGESVD(&jobu, &jobvt, &m, &n, a, &lda, sing, u, &ldu, vt, &ldvt, &work1, &lworka, info); 00262 ASSERT(*info==0) 00263 lworka = (int) work1; 00264 double* worka = SG_MALLOC(double, lworka); 00265 DGESVD(&jobu, &jobvt, &m, &n, a, &lda, sing, u, &ldu, vt, &ldvt, worka, &lworka, info); 00266 SG_FREE(worka); 00267 #endif 00268 } 00269 #undef DGESVD 00270 00271 void wrap_dgeqrf(int m, int n, double *a, int lda, double *tau, int *info) 00272 { 00273 #ifdef HAVE_ACML 00274 DGEQRF(m, n, a, lda, tau, info); 00275 #else 00276 int lwork = -1; 00277 double work1 = 0; 00278 DGEQRF(&m, &n, a, &lda, tau, &work1, &lwork, info); 00279 ASSERT(*info==0) 00280 lwork = (int)work1; 00281 ASSERT(lwork>0) 00282 double* work = SG_MALLOC(double, lwork); 00283 DGEQRF(&m, &n, a, &lda, tau, work, &lwork, info); 00284 SG_FREE(work); 00285 #endif 00286 } 00287 #undef DGEQRF 00288 00289 void wrap_dorgqr(int m, int n, int k, double *a, int lda, double *tau, int *info) 00290 { 00291 #ifdef HAVE_ACML 00292 DORGQR(m, n, k, a, lda, tau, info); 00293 #else 00294 int lwork = -1; 00295 double work1 = 0; 00296 DORGQR(&m, &n, &k, a, &lda, tau, &work1, &lwork, info); 00297 ASSERT(*info==0) 00298 lwork = (int)work1; 00299 ASSERT(lwork>0) 00300 double* work = SG_MALLOC(double, lwork); 00301 DORGQR(&m, &n, &k, a, &lda, tau, work, &lwork, info); 00302 SG_FREE(work); 00303 #endif 00304 } 00305 #undef DORGQR 00306 00307 void wrap_dsyevr(char jobz, char uplo, int n, double *a, int lda, int il, int iu, 00308 double *eigenvalues, double *eigenvectors, int *info) 00309 { 00310 int m; 00311 double vl,vu; 00312 double abstol = 0.0; 00313 char I = 'I'; 00314 int* isuppz = SG_MALLOC(int, n); 00315 #ifdef HAVE_ACML 00316 DSYEVR(jobz,I,uplo,n,a,lda,vl,vu,il,iu,abstol,m, 00317 eigenvalues,eigenvectors,n,isuppz,info); 00318 #else 00319 int lwork = -1; 00320 int liwork = -1; 00321 double work1 = 0; 00322 int work2 = 0; 00323 DSYEVR(&jobz,&I,&uplo,&n,a,&lda,&vl,&vu,&il,&iu,&abstol, 00324 &m,eigenvalues,eigenvectors,&n,isuppz, 00325 &work1,&lwork,&work2,&liwork,info); 00326 ASSERT(*info==0) 00327 lwork = (int)work1; 00328 liwork = work2; 00329 double* work = SG_MALLOC(double, lwork); 00330 int* iwork = SG_MALLOC(int, liwork); 00331 DSYEVR(&jobz,&I,&uplo,&n,a,&lda,&vl,&vu,&il,&iu,&abstol, 00332 &m,eigenvalues,eigenvectors,&n,isuppz, 00333 work,&lwork,iwork,&liwork,info); 00334 ASSERT(*info==0) 00335 SG_FREE(work); 00336 SG_FREE(iwork); 00337 SG_FREE(isuppz); 00338 #endif 00339 } 00340 #undef DSYEVR 00341 00342 void wrap_dsygvx(int itype, char jobz, char uplo, int n, double *a, int lda, double *b, 00343 int ldb, int il, int iu, double *eigenvalues, double *eigenvectors, int *info) 00344 { 00345 int m; 00346 double abstol = 0.0; 00347 double vl,vu; 00348 int* ifail = SG_MALLOC(int, n); 00349 char I = 'I'; 00350 #ifdef HAVE_ACML 00351 DSYGVX(itype,jobz,I,uplo,n,a,lda,b,ldb,vl,vu, 00352 il,iu,abstol,m,eigenvalues, 00353 eigenvectors,n,ifail,info); 00354 #else 00355 int lwork = -1; 00356 double work1 = 0; 00357 int* iwork = SG_MALLOC(int, 5*n); 00358 DSYGVX(&itype,&jobz,&I,&uplo,&n,a,&lda,b,&ldb,&vl,&vu, 00359 &il,&iu,&abstol,&m,eigenvalues,eigenvectors, 00360 &n,&work1,&lwork,iwork,ifail,info); 00361 lwork = (int)work1; 00362 double* work = SG_MALLOC(double, lwork); 00363 DSYGVX(&itype,&jobz,&I,&uplo,&n,a,&lda,b,&ldb,&vl,&vu, 00364 &il,&iu,&abstol,&m,eigenvalues,eigenvectors, 00365 &n,work,&lwork,iwork,ifail,info); 00366 SG_FREE(work); 00367 SG_FREE(iwork); 00368 SG_FREE(ifail); 00369 #endif 00370 } 00371 #undef DSYGVX 00372 00373 void wrap_dstemr(char jobz, char range, int n, double* diag, double *subdiag, 00374 double vl, double vu, int il, int iu, int* m, double* w, double* z__, 00375 int ldz, int nzc, int *isuppz, int tryrac, int *info) 00376 { 00377 #ifdef HAVE_ACML 00378 SG_SNOTIMPLEMENTED 00379 #else 00380 int lwork=-1; 00381 int liwork=-1; 00382 double work1=0; 00383 int iwork1=0; 00384 DSTEMR(&jobz, &range, &n, diag, subdiag, &vl, &vu, 00385 &il, &iu, m, w, z__, &ldz, &nzc, isuppz, &tryrac, 00386 &work1, &lwork, &iwork1, &liwork, info); 00387 lwork=(int)work1; 00388 liwork=iwork1; 00389 double* work=SG_MALLOC(double, lwork); 00390 int* iwork=SG_MALLOC(int, liwork); 00391 DSTEMR(&jobz, &range, &n, diag, subdiag, &vl, &vu, 00392 &il, &iu, m, w, z__, &ldz, &nzc, isuppz, &tryrac, 00393 work, &lwork, iwork, &liwork, info); 00394 00395 SG_FREE(work); 00396 SG_FREE(iwork); 00397 #endif 00398 } 00399 #undef DSTEMR 00400 00401 } 00402 #endif //HAVE_LAPACK