GetFEM  5.4.2
gmm_lapack_interface.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2003-2020 Yves Renard
5 
6  This file is a part of GetFEM
7 
8  GetFEM is free software; you can redistribute it and/or modify it
9  under the terms of the GNU Lesser General Public License as published
10  by the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version along with the GCC Runtime Library
12  Exception either version 3.1 or (at your option) any later version.
13  This program is distributed in the hope that it will be useful, but
14  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  License and GCC Runtime Library Exception for more details.
17  You should have received a copy of the GNU Lesser General Public License
18  along with this program; if not, write to the Free Software Foundation,
19  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
20 
21  As a special exception, you may use this file as it is a part of a free
22  software library without restriction. Specifically, if other files
23  instantiate templates or use macros or inline functions from this file,
24  or you compile this file and link it with other files to produce an
25  executable, this file does not by itself cause the resulting executable
26  to be covered by the GNU Lesser General Public License. This exception
27  does not however invalidate any other reasons why the executable file
28  might be covered by the GNU Lesser General Public License.
29 
30 ===========================================================================*/
31 
32 /**@file gmm_lapack_interface.h
33  @author Yves Renard <Yves.Renard@insa-lyon.fr>
34  @date October 7, 2003.
35  @brief gmm interface for LAPACK
36 */
37 
38 #ifndef GMM_LAPACK_INTERFACE_H
39 #define GMM_LAPACK_INTERFACE_H
40 
41 #include "gmm_blas_interface.h"
42 #include "gmm_dense_lu.h"
43 #include "gmm_dense_qr.h"
44 
45 
46 #if defined(GMM_USES_LAPACK)
47 
48 namespace gmm {
49 
50  /* ********************************************************************** */
51  /* Operations interfaced for T = float, double, std::complex<float> */
52  /* or std::complex<double> : */
53  /* */
54  /* lu_factor(dense_matrix<T>, std::vector<long>) */
55  /* lu_solve(dense_matrix<T>, std::vector<T>, std::vector<T>) */
56  /* lu_solve(dense_matrix<T>, std::vector<long>, std::vector<T>, */
57  /* std::vector<T>) */
58  /* lu_solve_transposed(dense_matrix<T>, std::vector<long>, std::vector<T>,*/
59  /* std::vector<T>) */
60  /* lu_inverse(dense_matrix<T>) */
61  /* lu_inverse(dense_matrix<T>, std::vector<long>, dense_matrix<T>) */
62  /* */
63  /* qr_factor(dense_matrix<T>, dense_matrix<T>, dense_matrix<T>) */
64  /* */
65  /* implicit_qr_algorithm(dense_matrix<T>, std::vector<T>) */
66  /* implicit_qr_algorithm(dense_matrix<T>, std::vector<T>, */
67  /* dense_matrix<T>) */
68  /* implicit_qr_algorithm(dense_matrix<T>, std::vector<std::complex<T> >) */
69  /* implicit_qr_algorithm(dense_matrix<T>, std::vector<std::complex<T> >, */
70  /* dense_matrix<T>) */
71  /* */
72  /* geev_interface_right */
73  /* geev_interface_left */
74  /* */
75  /* schur(dense_matrix<T>, dense_matrix<T>, dense_matrix<T>) */
76  /* */
77  /* svd(dense_matrix<T>, dense_matrix<T>, dense_matrix<T>, std::vector<T>) */
78  /* svd(dense_matrix<T>, dense_matrix<T>, dense_matrix<T>, */
79  /* std::vector<std::complex<T> >) */
80  /* */
81  /* ********************************************************************** */
82 
83  /* ********************************************************************** */
84  /* LAPACK functions used. */
85  /* ********************************************************************** */
86 
87  extern "C" {
88  void sgetrf_(...); void dgetrf_(...); void cgetrf_(...); void zgetrf_(...);
89  void sgetrs_(...); void dgetrs_(...); void cgetrs_(...); void zgetrs_(...);
90  void sgetri_(...); void dgetri_(...); void cgetri_(...); void zgetri_(...);
91  void sgeqrf_(...); void dgeqrf_(...); void cgeqrf_(...); void zgeqrf_(...);
92  void sorgqr_(...); void dorgqr_(...); void cungqr_(...); void zungqr_(...);
93  void sormqr_(...); void dormqr_(...); void cunmqr_(...); void zunmqr_(...);
94  void sgees_ (...); void dgees_ (...); void cgees_ (...); void zgees_ (...);
95  void sgeev_ (...); void dgeev_ (...); void cgeev_ (...); void zgeev_ (...);
96  void sgeesx_(...); void dgeesx_(...); void cgeesx_(...); void zgeesx_(...);
97  void sgesvd_(...); void dgesvd_(...); void cgesvd_(...); void zgesvd_(...);
98  }
99 
100  /* ********************************************************************** */
101  /* LU decomposition. */
102  /* ********************************************************************** */
103 
104 # define getrf_interface(lapack_name, base_type) inline \
105  size_type lu_factor(dense_matrix<base_type > &A, lapack_ipvt &ipvt){ \
106  GMMLAPACK_TRACE("getrf_interface"); \
107  BLAS_INT m = BLAS_INT(mat_nrows(A)), n = BLAS_INT(mat_ncols(A)), lda(m); \
108  long info(-1L); \
109  if (m && n) lapack_name(&m, &n, &A(0,0), &lda, ipvt.pfirst(), &info); \
110  if ((sizeof(BLAS_INT) == 4) || \
111  ((info & 0xFFFFFFFF00000000L) && !(info & 0x00000000FFFFFFFFL))) \
112  /* For compatibility with lapack version with 32 bit integer. */ \
113  ipvt.set_to_int32(); \
114  return size_type(int(info & 0x00000000FFFFFFFFL)); \
115  }
116 
117  getrf_interface(sgetrf_, BLAS_S)
118  getrf_interface(dgetrf_, BLAS_D)
119  getrf_interface(cgetrf_, BLAS_C)
120  getrf_interface(zgetrf_, BLAS_Z)
121 
122  /* ********************************************************************* */
123  /* LU solve. */
124  /* ********************************************************************* */
125 
126 # define getrs_interface(f_name, trans1, lapack_name, base_type) inline \
127  void f_name(const dense_matrix<base_type > &A, \
128  const lapack_ipvt &ipvt, std::vector<base_type > &x, \
129  const std::vector<base_type > &b) { \
130  GMMLAPACK_TRACE("getrs_interface"); \
131  BLAS_INT n = BLAS_INT(mat_nrows(A)), info(0), nrhs(1); \
132  gmm::copy(b, x); trans1; \
133  if (n) \
134  lapack_name(&t,&n,&nrhs,&(A(0,0)),&n,ipvt.pfirst(),&x[0],&n,&info); \
135  }
136 
137 # define getrs_trans_n const char t = 'N'
138 # define getrs_trans_t const char t = 'T'
139 
140  getrs_interface(lu_solve, getrs_trans_n, sgetrs_, BLAS_S)
141  getrs_interface(lu_solve, getrs_trans_n, dgetrs_, BLAS_D)
142  getrs_interface(lu_solve, getrs_trans_n, cgetrs_, BLAS_C)
143  getrs_interface(lu_solve, getrs_trans_n, zgetrs_, BLAS_Z)
144  getrs_interface(lu_solve_transposed, getrs_trans_t, sgetrs_, BLAS_S)
145  getrs_interface(lu_solve_transposed, getrs_trans_t, dgetrs_, BLAS_D)
146  getrs_interface(lu_solve_transposed, getrs_trans_t, cgetrs_, BLAS_C)
147  getrs_interface(lu_solve_transposed, getrs_trans_t, zgetrs_, BLAS_Z)
148 
149  /* ********************************************************************* */
150  /* LU inverse. */
151  /* ********************************************************************* */
152 
153 # define getri_interface(lapack_name, base_type) inline \
154  void lu_inverse(const dense_matrix<base_type > &LU, \
155  const lapack_ipvt &ipvt, const dense_matrix<base_type > &A_) { \
156  GMMLAPACK_TRACE("getri_interface"); \
157  dense_matrix<base_type >& \
158  A = const_cast<dense_matrix<base_type > &>(A_); \
159  BLAS_INT n = BLAS_INT(mat_nrows(A)), info(0), lwork(-1); \
160  base_type work1; \
161  if (n) { \
162  gmm::copy(LU, A); \
163  lapack_name(&n, &A(0,0), &n, ipvt.pfirst(), &work1, &lwork, &info); \
164  lwork = int(gmm::real(work1)); \
165  std::vector<base_type> work(lwork); \
166  lapack_name(&n, &A(0,0), &n, ipvt.pfirst(), &work[0], &lwork,&info); \
167  } \
168  }
169 
170  getri_interface(sgetri_, BLAS_S)
171  getri_interface(dgetri_, BLAS_D)
172  getri_interface(cgetri_, BLAS_C)
173  getri_interface(zgetri_, BLAS_Z)
174 
175  /* ********************************************************************** */
176  /* QR factorization. */
177  /* ********************************************************************** */
178 
179 # define geqrf_interface(lapack_name1, base_type) inline \
180  void qr_factor(dense_matrix<base_type > &A){ \
181  GMMLAPACK_TRACE("geqrf_interface"); \
182  BLAS_INT m = BLAS_INT(mat_nrows(A)), n=BLAS_INT(mat_ncols(A)); \
183  BLAS_INT info(0), lwork(-1); \
184  base_type work1; \
185  if (m && n) { \
186  std::vector<base_type > tau(n); \
187  lapack_name1(&m, &n, &A(0,0), &m, &tau[0], &work1 , &lwork, &info); \
188  lwork = BLAS_INT(gmm::real(work1)); \
189  std::vector<base_type > work(lwork); \
190  lapack_name1(&m, &n, &A(0,0), &m, &tau[0], &work[0], &lwork, &info); \
191  GMM_ASSERT1(!info, "QR factorization failed"); \
192  } \
193  }
194 
195  geqrf_interface(sgeqrf_, BLAS_S)
196  geqrf_interface(dgeqrf_, BLAS_D)
197  // For complex values, housholder vectors are not the same as in
198  // gmm::lu_factor. Impossible to interface for the moment.
199  // geqrf_interface(cgeqrf_, BLAS_C)
200  // geqrf_interface(zgeqrf_, BLAS_Z)
201 
202 # define geqrf_interface2(lapack_name1, lapack_name2, base_type) inline \
203  void qr_factor(const dense_matrix<base_type > &A, \
204  dense_matrix<base_type > &Q, dense_matrix<base_type > &R) { \
205  GMMLAPACK_TRACE("geqrf_interface2"); \
206  BLAS_INT m = BLAS_INT(mat_nrows(A)), n=BLAS_INT(mat_ncols(A)); \
207  BLAS_INT info(0), lwork(-1); \
208  base_type work1; \
209  if (m && n) { \
210  std::copy(A.begin(), A.end(), Q.begin()); \
211  std::vector<base_type > tau(n); \
212  lapack_name1(&m, &n, &Q(0,0), &m, &tau[0], &work1 , &lwork, &info); \
213  lwork = BLAS_INT(gmm::real(work1)); \
214  std::vector<base_type > work(lwork); \
215  lapack_name1(&m, &n, &Q(0,0), &m, &tau[0], &work[0], &lwork, &info); \
216  GMM_ASSERT1(!info, "QR factorization failed"); \
217  base_type *p = &R(0,0), *q = &Q(0,0); \
218  for (BLAS_INT j = 0; j < n; ++j, q += m-n) \
219  for (BLAS_INT i = 0; i < n; ++i, ++p, ++q) \
220  *p = (j < i) ? base_type(0) : *q; \
221  lapack_name2(&m, &n, &n, &Q(0,0), &m,&tau[0],&work[0],&lwork,&info); \
222  } \
223  else gmm::clear(Q); \
224  }
225 
226  geqrf_interface2(sgeqrf_, sorgqr_, BLAS_S)
227  geqrf_interface2(dgeqrf_, dorgqr_, BLAS_D)
228  geqrf_interface2(cgeqrf_, cungqr_, BLAS_C)
229  geqrf_interface2(zgeqrf_, zungqr_, BLAS_Z)
230 
231  /* ********************************************************************** */
232  /* QR algorithm for eigenvalues search. */
233  /* ********************************************************************** */
234 
235 # define gees_interface(lapack_name, base_type) \
236  template <typename VECT> inline void implicit_qr_algorithm( \
237  const dense_matrix<base_type > &A, const VECT &eigval_, \
238  dense_matrix<base_type > &Q, \
239  double tol=gmm::default_tol(base_type()), bool compvect = true) { \
240  GMMLAPACK_TRACE("gees_interface"); \
241  typedef bool (*L_fp)(...); L_fp p = 0; \
242  BLAS_INT n=BLAS_INT(mat_nrows(A)), info(0), lwork(-1), sdim; \
243  base_type work1; \
244  if (!n) return; \
245  dense_matrix<base_type > H(n,n); gmm::copy(A, H); \
246  char jobvs = (compvect ? 'V' : 'N'), sort = 'N'; \
247  std::vector<double> rwork(n), eigv1(n), eigv2(n); \
248  lapack_name(&jobvs, &sort, p, &n, &H(0,0), &n, &sdim, &eigv1[0], \
249  &eigv2[0], &Q(0,0), &n, &work1, &lwork, &rwork[0], &info); \
250  lwork = BLAS_INT(gmm::real(work1)); \
251  std::vector<base_type > work(lwork); \
252  lapack_name(&jobvs, &sort, p, &n, &H(0,0), &n, &sdim, &eigv1[0], \
253  &eigv2[0], &Q(0,0), &n, &work[0], &lwork, &rwork[0],&info);\
254  GMM_ASSERT1(!info, "QR algorithm failed"); \
255  extract_eig(H, const_cast<VECT &>(eigval_), tol); \
256  }
257 
258 # define gees_interface2(lapack_name, base_type) \
259  template <typename VECT> inline void implicit_qr_algorithm( \
260  const dense_matrix<base_type > &A, const VECT &eigval_, \
261  dense_matrix<base_type > &Q, \
262  double tol=gmm::default_tol(base_type()), bool compvect = true) { \
263  GMMLAPACK_TRACE("gees_interface2"); \
264  typedef bool (*L_fp)(...); L_fp p = 0; \
265  BLAS_INT n=BLAS_INT(mat_nrows(A)), info(0), lwork(-1), sdim; \
266  base_type work1; \
267  if (!n) return; \
268  dense_matrix<base_type > H(n,n); gmm::copy(A, H); \
269  char jobvs = (compvect ? 'V' : 'N'), sort = 'N'; \
270  std::vector<double> rwork(n), eigvv(n*2); \
271  lapack_name(&jobvs, &sort, p, &n, &H(0,0), &n, &sdim, &eigvv[0], \
272  &Q(0,0), &n, &work1, &lwork, &rwork[0], &rwork[0], &info); \
273  lwork = BLAS_INT(gmm::real(work1)); \
274  std::vector<base_type > work(lwork); \
275  lapack_name(&jobvs, &sort, p, &n, &H(0,0), &n, &sdim, &eigvv[0], \
276  &Q(0,0), &n, &work[0], &lwork, &rwork[0], &rwork[0],&info);\
277  GMM_ASSERT1(!info, "QR algorithm failed"); \
278  extract_eig(H, const_cast<VECT &>(eigval_), tol); \
279  }
280 
281  gees_interface(sgees_, BLAS_S)
282  gees_interface(dgees_, BLAS_D)
283  gees_interface2(cgees_, BLAS_C)
284  gees_interface2(zgees_, BLAS_Z)
285 
286 
287 # define jobv_right char jobvl = 'N', jobvr = 'V';
288 # define jobv_left char jobvl = 'V', jobvr = 'N';
289 
290 # define geev_interface(lapack_name, base_type, side) \
291  template <typename VECT> inline void geev_interface_ ## side( \
292  const dense_matrix<base_type > &A, const VECT &eigval_, \
293  dense_matrix<base_type > &Q) { \
294  GMMLAPACK_TRACE("geev_interface"); \
295  BLAS_INT n = BLAS_INT(mat_nrows(A)), info(0), lwork(-1); \
296  base_type work1; \
297  if (!n) return; \
298  dense_matrix<base_type > H(n,n); gmm::copy(A, H); \
299  jobv_ ## side \
300  std::vector<base_type > eigvr(n), eigvi(n); \
301  lapack_name(&jobvl, &jobvr, &n, &H(0,0), &n, &eigvr[0], &eigvi[0], \
302  &Q(0,0), &n, &Q(0,0), &n, &work1, &lwork, &info); \
303  lwork = BLAS_INT(gmm::real(work1)); \
304  std::vector<base_type > work(lwork); \
305  lapack_name(&jobvl, &jobvr, &n, &H(0,0), &n, &eigvr[0], &eigvi[0], \
306  &Q(0,0), &n, &Q(0,0), &n, &work[0], &lwork, &info); \
307  GMM_ASSERT1(!info, "QR algorithm failed"); \
308  gmm::copy(eigvr, gmm::real_part(const_cast<VECT &>(eigval_))); \
309  gmm::copy(eigvi, gmm::imag_part(const_cast<VECT &>(eigval_))); \
310  }
311 
312 # define geev_interface2(lapack_name, base_type, side) \
313  template <typename VECT> inline void geev_interface_ ## side( \
314  const dense_matrix<base_type > &A, const VECT &eigval_, \
315  dense_matrix<base_type > &Q) { \
316  GMMLAPACK_TRACE("geev_interface"); \
317  BLAS_INT n = BLAS_INT(mat_nrows(A)), info(0), lwork(-1); \
318  base_type work1; \
319  if (!n) return; \
320  dense_matrix<base_type > H(n,n); gmm::copy(A, H); \
321  jobv_ ## side \
322  std::vector<base_type::value_type> rwork(2*n); \
323  std::vector<base_type> eigv(n); \
324  lapack_name(&jobvl, &jobvr, &n, &H(0,0), &n, &eigv[0], &Q(0,0), &n, \
325  &Q(0,0), &n, &work1, &lwork, &rwork[0], &info); \
326  lwork = BLAS_INT(gmm::real(work1)); \
327  std::vector<base_type > work(lwork); \
328  lapack_name(&jobvl, &jobvr, &n, &H(0,0), &n, &eigv[0], &Q(0,0), &n, \
329  &Q(0,0), &n, &work[0], &lwork, &rwork[0], &info); \
330  GMM_ASSERT1(!info, "QR algorithm failed"); \
331  gmm::copy(eigv, const_cast<VECT &>(eigval_)); \
332  }
333 
334  geev_interface(sgeev_, BLAS_S, right)
335  geev_interface(dgeev_, BLAS_D, right)
336  geev_interface2(cgeev_, BLAS_C, right)
337  geev_interface2(zgeev_, BLAS_Z, right)
338 
339  geev_interface(sgeev_, BLAS_S, left)
340  geev_interface(dgeev_, BLAS_D, left)
341  geev_interface2(cgeev_, BLAS_C, left)
342  geev_interface2(zgeev_, BLAS_Z, left)
343 
344 
345  /* ********************************************************************** */
346  /* SCHUR algorithm: */
347  /* A = Q*S*(Q^T), with Q orthogonal and S upper quasi-triangula */
348  /* ********************************************************************** */
349 
350 # define geesx_interface(lapack_name, base_type) inline \
351  void schur(dense_matrix<base_type> &A, \
352  dense_matrix<base_type> &S, \
353  dense_matrix<base_type> &Q) { \
354  GMMLAPACK_TRACE("geesx_interface"); \
355  BLAS_INT m = BLAS_INT(mat_nrows(A)), n = BLAS_INT(mat_ncols(A)); \
356  GMM_ASSERT1(m == n, "Schur decomposition requires square matrix"); \
357  char jobvs = 'V', sort = 'N', sense = 'N'; \
358  bool select = false; \
359  BLAS_INT lwork = 8*n, sdim = 0, liwork = 1; \
360  std::vector<base_type> work(lwork), wr(n), wi(n); \
361  std::vector<BLAS_INT> iwork(liwork); \
362  std::vector<BLAS_INT> bwork(1); \
363  resize(S, n, n); copy(A, S); \
364  resize(Q, n, n); \
365  base_type rconde(0), rcondv(0); \
366  BLAS_INT info(0); \
367  lapack_name(&jobvs, &sort, &select, &sense, &n, &S(0,0), &n, \
368  &sdim, &wr[0], &wi[0], &Q(0,0), &n, &rconde, &rcondv, \
369  &work[0], &lwork, &iwork[0], &liwork, &bwork[0], &info);\
370  GMM_ASSERT1(!info, "SCHUR algorithm failed"); \
371  }
372 
373 # define geesx_interface2(lapack_name, base_type) inline \
374  void schur(dense_matrix<base_type> &A, \
375  dense_matrix<base_type> &S, \
376  dense_matrix<base_type> &Q) { \
377  GMMLAPACK_TRACE("geesx_interface"); \
378  BLAS_INT m = BLAS_INT(mat_nrows(A)), n = BLAS_INT(mat_ncols(A)); \
379  GMM_ASSERT1(m == n, "Schur decomposition requires square matrix"); \
380  char jobvs = 'V', sort = 'N', sense = 'N'; \
381  bool select = false; \
382  BLAS_INT lwork = 8*n, sdim = 0; \
383  std::vector<base_type::value_type> rwork(lwork); \
384  std::vector<base_type> work(lwork), w(n); \
385  std::vector<BLAS_INT> bwork(1); \
386  resize(S, n, n); copy(A, S); \
387  resize(Q, n, n); \
388  base_type rconde(0), rcondv(0); \
389  BLAS_INT info(0); \
390  lapack_name(&jobvs, &sort, &select, &sense, &n, &S(0,0), &n, \
391  &sdim, &w[0], &Q(0,0), &n, &rconde, &rcondv, \
392  &work[0], &lwork, &rwork[0], &bwork[0], &info); \
393  GMM_ASSERT1(!info, "SCHUR algorithm failed"); \
394  }
395 
396  geesx_interface(sgeesx_, BLAS_S)
397  geesx_interface(dgeesx_, BLAS_D)
398  geesx_interface2(cgeesx_, BLAS_C)
399  geesx_interface2(zgeesx_, BLAS_Z)
400 
401  template <typename MAT>
402  void schur(const MAT &A_, MAT &S, MAT &Q) {
403  MAT A(A_);
404  schur(A, S, Q);
405  }
406 
407 
408  /* ********************************************************************** */
409  /* Interface to SVD. Does not correspond to a Gmm++ functionnality. */
410  /* Author : Sebastian Nowozin <sebastian.nowozin@tuebingen.mpg.de> */
411  /* ********************************************************************** */
412 
413 # define gesvd_interface(lapack_name, base_type) inline \
414  void svd(dense_matrix<base_type> &X, \
415  dense_matrix<base_type> &U, \
416  dense_matrix<base_type> &Vtransposed, \
417  std::vector<base_type> &sigma) { \
418  GMMLAPACK_TRACE("gesvd_interface"); \
419  BLAS_INT m = BLAS_INT(mat_nrows(X)), n = BLAS_INT(mat_ncols(X)); \
420  BLAS_INT mn_min = m < n ? m : n; \
421  sigma.resize(mn_min); \
422  std::vector<base_type> work(15 * mn_min); \
423  BLAS_INT lwork = BLAS_INT(work.size()); \
424  resize(U, m, m); \
425  resize(Vtransposed, n, n); \
426  char job = 'A'; \
427  BLAS_INT info(0); \
428  lapack_name(&job, &job, &m, &n, &X(0,0), &m, &sigma[0], &U(0,0), \
429  &m, &Vtransposed(0,0), &n, &work[0], &lwork, &info); \
430  }
431 
432 # define cgesvd_interface(lapack_name, base_type, base_type2) inline \
433  void svd(dense_matrix<base_type> &X, \
434  dense_matrix<base_type> &U, \
435  dense_matrix<base_type> &Vtransposed, \
436  std::vector<base_type2> &sigma) { \
437  GMMLAPACK_TRACE("gesvd_interface"); \
438  BLAS_INT m = BLAS_INT(mat_nrows(X)), n = BLAS_INT(mat_ncols(X)); \
439  BLAS_INT mn_min = m < n ? m : n; \
440  sigma.resize(mn_min); \
441  std::vector<base_type> work(15 * mn_min); \
442  std::vector<base_type2> rwork(5 * mn_min); \
443  BLAS_INT lwork = BLAS_INT(work.size()); \
444  resize(U, m, m); \
445  resize(Vtransposed, n, n); \
446  char job = 'A'; \
447  BLAS_INT info(0); \
448  lapack_name(&job, &job, &m, &n, &X(0,0), &m, &sigma[0], &U(0,0), \
449  &m, &Vtransposed(0,0), &n, &work[0], &lwork, \
450  &rwork[0], &info); \
451  }
452 
453  gesvd_interface(sgesvd_, BLAS_S)
454  gesvd_interface(dgesvd_, BLAS_D)
455  cgesvd_interface(cgesvd_, BLAS_C, BLAS_S)
456  cgesvd_interface(zgesvd_, BLAS_Z, BLAS_D)
457 
458  template <typename MAT, typename VEC>
459  void svd(const MAT &X_, MAT &U, MAT &Vtransposed, VEC &sigma) {
460  MAT X(X_);
461  svd(X, U, Vtransposed, sigma);
462  }
463 
464 
465 
466 
467 }
468 
469 #else
470 
471 namespace gmm
472 {
473 template <typename MAT>
474 void schur(const MAT &A_, MAT &S, MAT &Q)
475 {
476  GMM_ASSERT1(false, "Use of function schur(A,S,Q) requires GetFEM "
477  "to be built with Lapack");
478 }
479 
480 }// namespace gmm
481 
482 #endif // GMM_USES_LAPACK
483 
484 #endif // GMM_LAPACK_INTERFACE_H
gmm_dense_lu.h
LU factorizations and determinant computation for dense matrices.
gmm_blas_interface.h
gmm interface for fortran BLAS.
gmm_dense_qr.h
Dense QR factorization.