SHOGUN  v2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
lapack.cpp
Go to the documentation of this file.
1 /*
2  * This program is free software; you can redistribute it and/or modify
3  * it under the terms of the GNU General Public License as published by
4  * the Free Software Foundation; either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * Written (W) 1999-2009 Soeren Sonnenburg
8  * Written (W) 1999-2008 Gunnar Raetsch
9  * Written (W) 2006-2007 Mikio L. Braun
10  * Written (W) 2008 Jochen Garcke
11  * Written (W) 2011 Sergey Lisitsyn
12  * Copyright (C) 1999-2009 Fraunhofer Institute FIRST and Max-Planck-Society
13  */
14 
15 #include <shogun/lib/config.h>
16 
17 #ifdef HAVE_LAPACK
18 #include <shogun/mathematics/lapack.h>
19 #include <shogun/lib/common.h>
20 #include <shogun/base/Parallel.h>
21 #include <pthread.h>
22 #include <shogun/io/SGIO.h>
23 
24 using namespace shogun;
25 
26 #if defined(HAVE_MKL) || defined(HAVE_ACML)
27 #define DSYEV dsyev
28 #define DGESVD dgesvd
29 #define DPOSV dposv
30 #define DPOTRF dpotrf
31 #define DPOTRI dpotri
32 #define DGETRI dgetri
33 #define DGETRF dgetrf
34 #define DGEQRF dgeqrf
35 #define DORGQR dorgqr
36 #define DSYEVR dsyevr
37 #define DPOTRS dpotrs
38 #define DGETRS dgetrs
39 #define DSYGVX dsygvx
40 #else
41 #define DSYEV dsyev_
42 #define DGESVD dgesvd_
43 #define DPOSV dposv_
44 #define DPOTRF dpotrf_
45 #define DPOTRI dpotri_
46 #define DGETRI dgetri_
47 #define DGETRF dgetrf_
48 #define DGEQRF dgeqrf_
49 #define DORGQR dorgqr_
50 #define DSYEVR dsyevr_
51 #define DGETRS dgetrs_
52 #define DPOTRS dpotrs_
53 #define DSYGVX dsygvx_
54 #endif
55 
56 #ifndef HAVE_ATLAS
57 int clapack_dpotrf(const CBLAS_ORDER Order, const CBLAS_UPLO Uplo,
58  const int N, double *A, const int LDA)
59 {
60  char uplo = 'U';
61  int info = 0;
62  if (Order==CblasRowMajor)
63  {//A is symmetric, we switch Uplo to get result for CblasRowMajor
64  if (Uplo==CblasUpper)
65  uplo='L';
66  }
67  else if (Uplo==CblasLower)
68  {
69  uplo='L';
70  }
71 #ifdef HAVE_ACML
72  DPOTRF(uplo, N, A, LDA, &info);
73 #else
74  int n=N;
75  int lda=LDA;
76  DPOTRF(&uplo, &n, A, &lda, &info);
77 #endif
78  return info;
79 }
80 #undef DPOTRF
81 
82 int clapack_dpotri(const CBLAS_ORDER Order, const CBLAS_UPLO Uplo,
83  const int N, double *A, const int LDA)
84 {
85  char uplo = 'U';
86  int info = 0;
87  if (Order==CblasRowMajor)
88  {//A is symmetric, we switch Uplo to get result for CblasRowMajor
89  if (Uplo==CblasUpper)
90  uplo='L';
91  }
92  else if (Uplo==CblasLower)
93  {
94  uplo='L';
95  }
96 #ifdef HAVE_ACML
97  DPOTRI(uplo, N, A, LDA, &info);
98 #else
99  int n=N;
100  int lda=LDA;
101  DPOTRI(&uplo, &n, A, &lda, &info);
102 #endif
103  return info;
104 }
105 #undef DPOTRI
106 
107 /* DPOSV computes the solution to a real system of linear equations
108  * A * X = B,
109  * where A is an N-by-N symmetric positive definite matrix and X and B
110  * are N-by-NRHS matrices
111  */
112 int clapack_dposv(const CBLAS_ORDER Order, const CBLAS_UPLO Uplo,
113  const int N, const int NRHS, double *A, const int lda,
114  double *B, const int ldb)
115 {
116  char uplo = 'U';
117  int info=0;
118  if (Order==CblasRowMajor)
119  {//A is symmetric, we switch Uplo to achieve CblasColMajor
120  if (Uplo==CblasUpper)
121  uplo='L';
122  }
123  else if (Uplo==CblasLower)
124  {
125  uplo='L';
126  }
127 #ifdef HAVE_ACML
128  DPOSV(uplo,N,NRHS,A,lda,B,ldb,&info);
129 #else
130  int n=N;
131  int nrhs=NRHS;
132  int LDA=lda;
133  int LDB=ldb;
134  DPOSV(&uplo, &n, &nrhs, A, &LDA, B, &LDB, &info);
135 #endif
136  return info;
137 }
138 #undef DPOSV
139 
140 int clapack_dgetrf(const CBLAS_ORDER Order, const int M, const int N,
141  double *A, const int lda, int *ipiv)
142 {
143  // no rowmajor?
144  int info=0;
145 #ifdef HAVE_ACML
146  DGETRF(M,N,A,lda,ipiv,&info);
147 #else
148  int m=M;
149  int n=N;
150  int LDA=lda;
151  DGETRF(&m,&n,A,&LDA,ipiv,&info);
152 #endif
153  return info;
154 }
155 #undef DGETRF
156 
157 // order not supported (yet?)
158 int clapack_dgetri(const CBLAS_ORDER Order, const int N, double *A,
159  const int lda, int* ipiv)
160 {
161  int info=0;
162 #ifdef HAVE_ACML
163  DGETRI(N,A,lda,ipiv,&info);
164 #else
165  double* work;
166  int n=N;
167  int LDA=lda;
168  int lwork = -1;
169  double work1 = 0;
170  DGETRI(&n,A,&LDA,ipiv,&work1,&lwork,&info);
171  lwork = (int) work1;
172  work = SG_MALLOC(double, lwork);
173  DGETRI(&n,A,&LDA,ipiv,work,&lwork,&info);
174  SG_FREE(work);
175 #endif
176  return info;
177 }
178 #undef DGETRI
179 
180 // order not supported (yet?)
181 int clapack_dgetrs(const CBLAS_ORDER Order, const CBLAS_TRANSPOSE Transpose,
182  const int N, const int NRHS, double *A, const int lda,
183  int *ipiv, double *B, const int ldb)
184 {
185  int info = 0;
186  char trans = 'N';
187  if (Transpose==CblasTrans)
188  {
189  trans = 'T';
190  }
191 #ifdef HAVE_ACML
192  DGETRS(trans,N,NRHS,A,lda,ipiv,B,ldb,info);
193 #else
194  int n=N;
195  int nrhs=NRHS;
196  int LDA=lda;
197  int LDB=ldb;
198  DGETRS(&trans,&n,&nrhs,A,&LDA,ipiv,B,&LDB,&info);
199 #endif
200  return info;
201 }
202 #undef DGETRS
203 
204 // order not supported (yet?)
205 int clapack_dpotrs(const CBLAS_ORDER Order, const CBLAS_UPLO Uplo,
206  const int N, const int NRHS, double *A, const int lda,
207  double *B, const int ldb)
208 {
209  int info=0;
210  char uplo = 'U';
211  if (Uplo==CblasLower)
212  {
213  uplo = 'L';
214  }
215 #ifdef HAVE_ACML
216  DPOTRS(uplo,N,NRHS,A,lda,B,ldb,info);
217 #else
218  int n=N;
219  int nrhs=NRHS;
220  int LDA=lda;
221  int LDB=ldb;
222  DPOTRS(&uplo,&n,&nrhs,A,&LDA,B,&LDB,&info);
223 #endif
224  return info;
225 }
226 #undef DPOTRS
227 #endif //HAVE_ATLAS
228 
229 namespace shogun
230 {
231 
232 void wrap_dsyev(char jobz, char uplo, int n, double *a, int lda, double *w, int *info)
233 {
234 #ifdef HAVE_ACML
235  DSYEV(jobz, uplo, n, a, lda, w, info);
236 #else
237  int lwork=-1;
238  double work1;
239  DSYEV(&jobz, &uplo, &n, a, &lda, w, &work1, &lwork, info);
240  ASSERT(*info==0);
241  ASSERT(work1>0);
242  lwork=(int) work1;
243  double* work=SG_MALLOC(double, lwork);
244  DSYEV(&jobz, &uplo, &n, a, &lda, w, work, &lwork, info);
245  SG_FREE(work);
246 #endif
247 }
248 #undef DSYEV
249 
250 void wrap_dgesvd(char jobu, char jobvt, int m, int n, double *a, int lda, double *sing,
251  double *u, int ldu, double *vt, int ldvt, int *info)
252 {
253 #ifdef HAVE_ACML
254  DGESVD(jobu, jobvt, m, n, a, lda, sing, u, ldu, vt, ldvt, info);
255 #else
256  double work1 = 0;
257  int lworka = -1;
258  DGESVD(&jobu, &jobvt, &m, &n, a, &lda, sing, u, &ldu, vt, &ldvt, &work1, &lworka, info);
259  ASSERT(*info==0);
260  lworka = (int) work1;
261  double* worka = SG_MALLOC(double, lworka);
262  DGESVD(&jobu, &jobvt, &m, &n, a, &lda, sing, u, &ldu, vt, &ldvt, worka, &lworka, info);
263  SG_FREE(worka);
264 #endif
265 }
266 #undef DGESVD
267 
268 void wrap_dgeqrf(int m, int n, double *a, int lda, double *tau, int *info)
269 {
270 #ifdef HAVE_ACML
271  DGEQRF(m, n, a, lda, tau, info);
272 #else
273  int lwork = -1;
274  double work1 = 0;
275  DGEQRF(&m, &n, a, &lda, tau, &work1, &lwork, info);
276  ASSERT(*info==0);
277  lwork = (int)work1;
278  ASSERT(lwork>0)
279  double* work = SG_MALLOC(double, lwork);
280  DGEQRF(&m, &n, a, &lda, tau, work, &lwork, info);
281  SG_FREE(work);
282 #endif
283 }
284 #undef DGEQRF
285 
286 void wrap_dorgqr(int m, int n, int k, double *a, int lda, double *tau, int *info)
287 {
288 #ifdef HAVE_ACML
289  DORGQR(m, n, k, a, lda, tau, info);
290 #else
291  int lwork = -1;
292  double work1 = 0;
293  DORGQR(&m, &n, &k, a, &lda, tau, &work1, &lwork, info);
294  ASSERT(*info==0);
295  lwork = (int)work1;
296  ASSERT(lwork>0);
297  double* work = SG_MALLOC(double, lwork);
298  DORGQR(&m, &n, &k, a, &lda, tau, work, &lwork, info);
299  SG_FREE(work);
300 #endif
301 }
302 #undef DORGQR
303 
304 void wrap_dsyevr(char jobz, char uplo, int n, double *a, int lda, int il, int iu,
305  double *eigenvalues, double *eigenvectors, int *info)
306 {
307  int m;
308  double vl,vu;
309  double abstol = 0.0;
310  char I = 'I';
311  int* isuppz = SG_MALLOC(int, n);
312 #ifdef HAVE_ACML
313  DSYEVR(jobz,I,uplo,n,a,lda,vl,vu,il,iu,abstol,m,
314  eigenvalues,eigenvectors,n,isuppz,info);
315 #else
316  int lwork = -1;
317  int liwork = -1;
318  double work1 = 0;
319  int work2 = 0;
320  DSYEVR(&jobz,&I,&uplo,&n,a,&lda,&vl,&vu,&il,&iu,&abstol,
321  &m,eigenvalues,eigenvectors,&n,isuppz,
322  &work1,&lwork,&work2,&liwork,info);
323  ASSERT(*info==0);
324  lwork = (int)work1;
325  liwork = work2;
326  double* work = SG_MALLOC(double, lwork);
327  int* iwork = SG_MALLOC(int, liwork);
328  DSYEVR(&jobz,&I,&uplo,&n,a,&lda,&vl,&vu,&il,&iu,&abstol,
329  &m,eigenvalues,eigenvectors,&n,isuppz,
330  work,&lwork,iwork,&liwork,info);
331  ASSERT(*info==0);
332  SG_FREE(work);
333  SG_FREE(iwork);
334  SG_FREE(isuppz);
335 #endif
336 }
337 #undef DSYEVR
338 
339 void wrap_dsygvx(int itype, char jobz, char uplo, int n, double *a, int lda, double *b,
340  int ldb, int il, int iu, double *eigenvalues, double *eigenvectors, int *info)
341 {
342  int m;
343  double abstol = 0.0;
344  double vl,vu;
345  int* ifail = SG_MALLOC(int, n);
346  char I = 'I';
347 #ifdef HAVE_ACML
348  DSYGVX(itype,jobz,I,uplo,n,a,lda,b,ldb,vl,vu,
349  il,iu,abstol,m,eigenvalues,
350  eigenvectors,n,ifail,info);
351 #else
352  int lwork = -1;
353  double work1 = 0;
354  int* iwork = SG_MALLOC(int, 5*n);
355  DSYGVX(&itype,&jobz,&I,&uplo,&n,a,&lda,b,&ldb,&vl,&vu,
356  &il,&iu,&abstol,&m,eigenvalues,eigenvectors,
357  &n,&work1,&lwork,iwork,ifail,info);
358  lwork = (int)work1;
359  double* work = SG_MALLOC(double, lwork);
360  DSYGVX(&itype,&jobz,&I,&uplo,&n,a,&lda,b,&ldb,&vl,&vu,
361  &il,&iu,&abstol,&m,eigenvalues,eigenvectors,
362  &n,work,&lwork,iwork,ifail,info);
363  SG_FREE(work);
364  SG_FREE(iwork);
365  SG_FREE(ifail);
366 #endif
367 }
368 #undef DSYGVX
369 }
370 #endif //HAVE_LAPACK
SHOGUN Machine Learning Toolbox - Documentation