SHOGUN  v3.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
lapack.cpp
Go to the documentation of this file.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation