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