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  * Copyright (C) 1999-2009 Fraunhofer Institute FIRST and Max-Planck-Society
00012  */
00013 
00014 #include "lib/config.h"
00015 
00016 #ifdef HAVE_LAPACK
00017 #include "lib/lapack.h"
00018 #include "lib/common.h"
00019 #include "lib/io.h"
00020 
00021 #if defined(HAVE_MKL) || defined(HAVE_ACML) 
00022 #define DSYEV dsyev
00023 #define DGESVD dgesvd
00024 #define DPOSV dposv
00025 #define DPOTRF dpotrf
00026 #else
00027 #define DSYEV dsyev_
00028 #define DGESVD dgesvd_
00029 #define DPOSV dposv_
00030 #define DPOTRF dpotrf_
00031 #endif
00032 
00033 #ifndef HAVE_ATLAS
00034 int clapack_dpotrf(const CBLAS_ORDER Order, const CBLAS_UPLO Uplo,
00035         const int N, double *A, const int LDA)
00036 {
00037     char uplo = 'U';
00038     int info = 0;
00039     if (Order==CblasRowMajor)
00040     {//A is symmetric, we switch Uplo to get result for CblasRowMajor
00041         if (Uplo==CblasUpper)
00042             uplo='L';
00043     }
00044     else
00045         if (Uplo==CblasLower)
00046             uplo='L';
00047 #ifdef HAVE_ACML
00048     DPOTRF(uplo, N, A, LDA, &info);
00049 #else
00050     int n=N;
00051     int lda=LDA;
00052     DPOTRF(&uplo, &n, A, &lda, &info);
00053 #endif
00054     return info;
00055 }
00056 #undef DPOTRF
00057 
00058 /* DPOSV computes the solution to a real system of linear equations
00059  * A * X = B,
00060  * where A is an N-by-N symmetric positive definite matrix and X and B
00061  * are N-by-NRHS matrices
00062  */
00063 int clapack_dposv(const CBLAS_ORDER Order, const CBLAS_UPLO Uplo,
00064         const int N, const int NRHS, double *A, const int lda,
00065         double *B, const int ldb)
00066 {
00067     char uplo = 'U';
00068     int info=0;
00069     if (Order==CblasRowMajor)
00070     {//A is symmetric, we switch Uplo to achieve CblasColMajor
00071         if (Uplo==CblasUpper)
00072             uplo='L';
00073     }
00074     else
00075         if (Uplo==CblasLower)
00076             uplo='L';
00077 #ifdef HAVE_ACML
00078     DPOSV(uplo,N,NRHS,A,lda,B,ldb,&info);
00079 #else
00080     int n=N;
00081     int nrhs=NRHS;
00082     int LDA=lda;
00083     int LDB=ldb;
00084     DPOSV(&uplo, &n, &nrhs, A, &LDA, B, &LDB, &info);
00085 #endif
00086     return info;
00087 }
00088 #undef DPOSV
00089 #endif //HAVE_ATLAS
00090 
00091 /*
00092  * Wrapper files for LAPACK if there isn't a clapack interface
00093  * 
00094  */
00095 
00096 /*  DSYEV computes all eigenvalues and, optionally, eigenvectors of a
00097  *  real symmetric matrix A.
00098  */
00099 void wrap_dsyev(char jobz, char uplo, int n, double *a, int lda, double *w, int *info)
00100 {
00101 #ifdef HAVE_ACML
00102     DSYEV(jobz, uplo, n, a, lda, w, info);
00103 #else
00104     int lwork=-1;
00105     double work1;
00106     DSYEV(&jobz, &uplo, &n, a, &lda, w, &work1, &lwork, info);
00107     ASSERT(*info==0);
00108     ASSERT(work1>0);
00109     lwork=(int) work1;
00110     double* work=new double[lwork];
00111     DSYEV(&jobz, &uplo, &n, a, &lda, w, work, &lwork, info);
00112     delete[] work;
00113 #endif
00114 }
00115 #undef DSYEV
00116 
00117 void wrap_dgesvd(char jobu, char jobvt, int m, int n, double *a, int lda, double *sing, 
00118         double *u, int ldu, double *vt, int ldvt, int *info)
00119 {
00120 #ifdef HAVE_ACML
00121     DGESVD(jobu, jobvt, m, n, a, lda, sing, u, ldu, vt, ldvt, info);
00122 #else
00123     int lwork=-1;
00124     double work1;
00125     DGESVD(&jobu, &jobvt, &m, &n, a, &lda, sing, u, &ldu, vt, &ldvt, &work1, &lwork, info);
00126     ASSERT(*info==0);
00127     ASSERT(work1>0);
00128     lwork=(int) work1;
00129     double* work=new double[lwork];
00130     DGESVD(&jobu, &jobvt, &m, &n, a, &lda, sing, u, &ldu, vt, &ldvt, work, &lwork, info);
00131     delete[] work;
00132 #endif
00133 }
00134 #undef DGESVD
00135 #endif //HAVE_LAPACK

SHOGUN Machine Learning Toolbox - Documentation