SHOGUN
v2.0.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/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 {//A is symmetric, we switch Uplo to get result for CblasRowMajor 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 {//A is symmetric, we switch Uplo to get result for CblasRowMajor 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 /* DPOSV computes the solution to a real system of linear equations 00108 * A * X = B, 00109 * where A is an N-by-N symmetric positive definite matrix and X and B 00110 * are N-by-NRHS matrices 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 {//A is symmetric, we switch Uplo to achieve CblasColMajor 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 // no rowmajor? 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 // order not supported (yet?) 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 // order not supported (yet?) 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 // order not supported (yet?) 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