00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
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 {
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
00059
00060
00061
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 {
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
00093
00094
00095
00096
00097
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