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