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/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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation