GetFEM  5.4.3
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 <[email protected]>
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 #if defined(GMM_USES_LAPACK) && !defined(GMM_MATLAB_INTERFACE)
46 
47 namespace gmm {
48 
49  /* ********************************************************************** */
50  /* Operations interfaced for T = float, double, std::complex<float> */
51  /* or std::complex<double> : */
52  /* */
53  /* lu_factor(dense_matrix<T>, std::vector<long>) */
54  /* lu_solve(dense_matrix<T>, std::vector<T>, std::vector<T>) */
55  /* lu_solve(dense_matrix<T>, std::vector<long>, std::vector<T>, */
56  /* std::vector<T>) */
57  /* lu_solve_transposed(dense_matrix<T>, std::vector<long>, std::vector<T>,*/
58  /* std::vector<T>) */
59  /* lu_inverse(dense_matrix<T>) */
60  /* lu_inverse(dense_matrix<T>, std::vector<long>, dense_matrix<T>) */
61  /* */
62  /* qr_factor(dense_matrix<T>, dense_matrix<T>, dense_matrix<T>) */
63  /* */
64  /* implicit_qr_algorithm(dense_matrix<T>, std::vector<T>) */
65  /* implicit_qr_algorithm(dense_matrix<T>, std::vector<T>, */
66  /* dense_matrix<T>) */
67  /* implicit_qr_algorithm(dense_matrix<T>, std::vector<std::complex<T> >) */
68  /* implicit_qr_algorithm(dense_matrix<T>, std::vector<std::complex<T> >, */
69  /* dense_matrix<T>) */
70  /* */
71  /* geev_interface_right */
72  /* geev_interface_left */
73  /* */
74  /* schur(dense_matrix<T>, dense_matrix<T>, dense_matrix<T>) */
75  /* */
76  /* svd(dense_matrix<T>, dense_matrix<T>, dense_matrix<T>, std::vector<T>) */
77  /* svd(dense_matrix<T>, dense_matrix<T>, dense_matrix<T>, */
78  /* std::vector<std::complex<T> >) */
79  /* */
80  /* ********************************************************************** */
81 
82  /* ********************************************************************** */
83  /* LAPACK functions used. */
84  /* ********************************************************************** */
85 
86  extern "C" {
87  void sgetrf_(...); void dgetrf_(...); void cgetrf_(...); void zgetrf_(...);
88  void sgetrs_(...); void dgetrs_(...); void cgetrs_(...); void zgetrs_(...);
89  void sgetri_(...); void dgetri_(...); void cgetri_(...); void zgetri_(...);
90  void sgeqrf_(...); void dgeqrf_(...); void cgeqrf_(...); void zgeqrf_(...);
91  void sorgqr_(...); void dorgqr_(...); void cungqr_(...); void zungqr_(...);
92  void sormqr_(...); void dormqr_(...); void cunmqr_(...); void zunmqr_(...);
93  void sgees_ (...); void dgees_ (...); void cgees_ (...); void zgees_ (...);
94  void sgeev_ (...); void dgeev_ (...); void cgeev_ (...); void zgeev_ (...);
95  void sgeesx_(...); void dgeesx_(...); void cgeesx_(...); void zgeesx_(...);
96  void sgesvd_(...); void dgesvd_(...); void cgesvd_(...); void zgesvd_(...);
97  }
98 
99  /* ********************************************************************** */
100  /* LU decomposition. */
101  /* ********************************************************************** */
102 
103 # define getrf_interface(lapack_name, base_type) inline \
104  size_type lu_factor(dense_matrix<base_type> &A, lapack_ipvt &ipvt){ \
105  GMMLAPACK_TRACE("getrf_interface"); \
106  BLAS_INT m = BLAS_INT(mat_nrows(A)), n = BLAS_INT(mat_ncols(A)), lda(m); \
107  BLAS_INT info=BLAS_INT(-1); \
108  if (m && n) lapack_name(&m, &n, &A(0,0), &lda, &ipvt[0], &info); \
109  return size_type(abs(info)); \
110  }
111 
112  getrf_interface(sgetrf_, BLAS_S)
113  getrf_interface(dgetrf_, BLAS_D)
114  getrf_interface(cgetrf_, BLAS_C)
115  getrf_interface(zgetrf_, BLAS_Z)
116 
117  /* ********************************************************************* */
118  /* LU solve. */
119  /* ********************************************************************* */
120 
121 # define getrs_interface(f_name, trans1, lapack_name, base_type) inline \
122  void f_name(const dense_matrix<base_type> &A, \
123  const lapack_ipvt &ipvt, std::vector<base_type> &x, \
124  const std::vector<base_type> &b) { \
125  GMMLAPACK_TRACE("getrs_interface"); \
126  BLAS_INT n = BLAS_INT(mat_nrows(A)), info(0), nrhs(1); \
127  gmm::copy(b, x); trans1; \
128  if (n) \
129  lapack_name(&t,&n,&nrhs,&(A(0,0)),&n,&ipvt[0],&x[0],&n,&info); \
130  }
131 
132 # define getrs_trans_n const char t = 'N'
133 # define getrs_trans_t const char t = 'T'
134 
135  getrs_interface(lu_solve, getrs_trans_n, sgetrs_, BLAS_S)
136  getrs_interface(lu_solve, getrs_trans_n, dgetrs_, BLAS_D)
137  getrs_interface(lu_solve, getrs_trans_n, cgetrs_, BLAS_C)
138  getrs_interface(lu_solve, getrs_trans_n, zgetrs_, BLAS_Z)
139  getrs_interface(lu_solve_transposed, getrs_trans_t, sgetrs_, BLAS_S)
140  getrs_interface(lu_solve_transposed, getrs_trans_t, dgetrs_, BLAS_D)
141  getrs_interface(lu_solve_transposed, getrs_trans_t, cgetrs_, BLAS_C)
142  getrs_interface(lu_solve_transposed, getrs_trans_t, zgetrs_, BLAS_Z)
143 
144  /* ********************************************************************* */
145  /* LU inverse. */
146  /* ********************************************************************* */
147 
148 # define getri_interface(lapack_name, base_type) inline \
149  void lu_inverse(const dense_matrix<base_type> &LU, \
150  const lapack_ipvt &ipvt, dense_matrix<base_type> &A) { \
151  GMMLAPACK_TRACE("getri_interface"); \
152  BLAS_INT n = BLAS_INT(mat_nrows(A)), info(0), lwork(-1); \
153  base_type work1; \
154  if (n) { \
155  gmm::copy(LU, A); \
156  lapack_name(&n, &A(0,0), &n, &ipvt[0], &work1, &lwork, &info); \
157  lwork = int(gmm::real(work1)); \
158  std::vector<base_type> work(lwork); \
159  lapack_name(&n, &A(0,0), &n, &ipvt[0], &work[0], &lwork, &info); \
160  } \
161  }
162 
163  getri_interface(sgetri_, BLAS_S)
164  getri_interface(dgetri_, BLAS_D)
165  getri_interface(cgetri_, BLAS_C)
166  getri_interface(zgetri_, BLAS_Z)
167 
168  /* ********************************************************************** */
169  /* QR factorization. */
170  /* ********************************************************************** */
171 
172 # define geqrf_interface(lapack_name1, base_type) inline \
173  void qr_factor(dense_matrix<base_type> &A){ \
174  GMMLAPACK_TRACE("geqrf_interface"); \
175  BLAS_INT m = BLAS_INT(mat_nrows(A)), n=BLAS_INT(mat_ncols(A)); \
176  BLAS_INT info(0), lwork(-1); \
177  base_type work1; \
178  if (m && n) { \
179  std::vector<base_type> tau(n); \
180  lapack_name1(&m, &n, &A(0,0), &m, &tau[0], &work1 , &lwork, &info); \
181  lwork = BLAS_INT(gmm::real(work1)); \
182  std::vector<base_type> work(lwork); \
183  lapack_name1(&m, &n, &A(0,0), &m, &tau[0], &work[0], &lwork, &info); \
184  GMM_ASSERT1(!info, "QR factorization failed"); \
185  } \
186  }
187 
188  geqrf_interface(sgeqrf_, BLAS_S)
189  geqrf_interface(dgeqrf_, BLAS_D)
190  // For complex values, housholder vectors are not the same as in
191  // gmm::lu_factor. Impossible to interface for the moment.
192  // geqrf_interface(cgeqrf_, BLAS_C)
193  // geqrf_interface(zgeqrf_, BLAS_Z)
194 
195 # define geqrf_interface2(lapack_name1, lapack_name2, base_type) inline \
196  void qr_factor(const dense_matrix<base_type> &A, \
197  dense_matrix<base_type> &Q, dense_matrix<base_type> &R) { \
198  GMMLAPACK_TRACE("geqrf_interface2"); \
199  BLAS_INT m = BLAS_INT(mat_nrows(A)), n=BLAS_INT(mat_ncols(A)); \
200  BLAS_INT info(0), lwork(-1); \
201  base_type work1; \
202  if (m && n) { \
203  std::copy(A.begin(), A.end(), Q.begin()); \
204  std::vector<base_type> tau(n); \
205  lapack_name1(&m, &n, &Q(0,0), &m, &tau[0], &work1 , &lwork, &info); \
206  lwork = BLAS_INT(gmm::real(work1)); \
207  std::vector<base_type> work(lwork); \
208  lapack_name1(&m, &n, &Q(0,0), &m, &tau[0], &work[0], &lwork, &info); \
209  GMM_ASSERT1(!info, "QR factorization failed"); \
210  base_type *p = &R(0,0), *q = &Q(0,0); \
211  for (BLAS_INT j = 0; j < n; ++j, q += m-n) \
212  for (BLAS_INT i = 0; i < n; ++i, ++p, ++q) \
213  *p = (j < i) ? base_type(0) : *q; \
214  lapack_name2(&m, &n, &n, &Q(0,0), &m,&tau[0],&work[0],&lwork,&info); \
215  } \
216  else gmm::clear(Q); \
217  }
218 
219  geqrf_interface2(sgeqrf_, sorgqr_, BLAS_S)
220  geqrf_interface2(dgeqrf_, dorgqr_, BLAS_D)
221  geqrf_interface2(cgeqrf_, cungqr_, BLAS_C)
222  geqrf_interface2(zgeqrf_, zungqr_, BLAS_Z)
223 
224  /* ********************************************************************** */
225  /* QR algorithm for eigenvalues search. */
226  /* ********************************************************************** */
227 
228 # define gees_interface(lapack_name, base_type) \
229  template <typename VECT> inline void implicit_qr_algorithm( \
230  const dense_matrix<base_type> &A, VECT &eigval_, \
231  dense_matrix<base_type> &Q, \
232  double tol=gmm::default_tol(base_type()), bool compvect = true) { \
233  GMMLAPACK_TRACE("gees_interface"); \
234  typedef bool (*L_fp)(...); L_fp p = 0; \
235  BLAS_INT n=BLAS_INT(mat_nrows(A)), info(0), lwork(-1), sdim; \
236  base_type work1; \
237  if (!n) return; \
238  dense_matrix<base_type> H(n,n); gmm::copy(A, H); \
239  char jobvs = (compvect ? 'V' : 'N'), sort = 'N'; \
240  std::vector<double> rwork(n), eigv1(n), eigv2(n); \
241  lapack_name(&jobvs, &sort, p, &n, &H(0,0), &n, &sdim, &eigv1[0], \
242  &eigv2[0], &Q(0,0), &n, &work1, &lwork, &rwork[0], &info); \
243  lwork = BLAS_INT(gmm::real(work1)); \
244  std::vector<base_type> work(lwork); \
245  lapack_name(&jobvs, &sort, p, &n, &H(0,0), &n, &sdim, &eigv1[0], \
246  &eigv2[0], &Q(0,0), &n, &work[0], &lwork, &rwork[0],&info);\
247  GMM_ASSERT1(!info, "QR algorithm failed"); \
248  extract_eig(H, eigval_, tol); \
249  }
250 
251 # define gees_interface2(lapack_name, base_type) \
252  template <typename VECT> inline void implicit_qr_algorithm( \
253  const dense_matrix<base_type> &A, VECT &eigval_, \
254  dense_matrix<base_type> &Q, \
255  double tol=gmm::default_tol(base_type()), bool compvect = true) { \
256  GMMLAPACK_TRACE("gees_interface2"); \
257  typedef bool (*L_fp)(...); L_fp p = 0; \
258  BLAS_INT n=BLAS_INT(mat_nrows(A)), info(0), lwork(-1), sdim; \
259  base_type work1; \
260  if (!n) return; \
261  dense_matrix<base_type> H(n,n); gmm::copy(A, H); \
262  char jobvs = (compvect ? 'V' : 'N'), sort = 'N'; \
263  std::vector<double> rwork(n), eigvv(n*2); \
264  lapack_name(&jobvs, &sort, p, &n, &H(0,0), &n, &sdim, &eigvv[0], \
265  &Q(0,0), &n, &work1, &lwork, &rwork[0], &rwork[0], &info); \
266  lwork = BLAS_INT(gmm::real(work1)); \
267  std::vector<base_type> work(lwork); \
268  lapack_name(&jobvs, &sort, p, &n, &H(0,0), &n, &sdim, &eigvv[0], \
269  &Q(0,0), &n, &work[0], &lwork, &rwork[0], &rwork[0],&info);\
270  GMM_ASSERT1(!info, "QR algorithm failed"); \
271  extract_eig(H, eigval_, tol); \
272  }
273 
274  gees_interface(sgees_, BLAS_S)
275  gees_interface(dgees_, BLAS_D)
276  gees_interface2(cgees_, BLAS_C)
277  gees_interface2(zgees_, BLAS_Z)
278 
279 
280 # define jobv_right char jobvl = 'N', jobvr = 'V';
281 # define jobv_left char jobvl = 'V', jobvr = 'N';
282 
283 # define geev_interface(lapack_name, base_type, side) \
284  template <typename VECT> inline void geev_interface_ ## side( \
285  const dense_matrix<base_type> &A, VECT &eigval_, \
286  dense_matrix<base_type> &Q) { \
287  GMMLAPACK_TRACE("geev_interface"); \
288  BLAS_INT n = BLAS_INT(mat_nrows(A)), info(0), lwork(-1); \
289  base_type work1; \
290  if (!n) return; \
291  dense_matrix<base_type> H(n,n); gmm::copy(A, H); \
292  jobv_ ## side \
293  std::vector<base_type> eigvr(n), eigvi(n); \
294  lapack_name(&jobvl, &jobvr, &n, &H(0,0), &n, &eigvr[0], &eigvi[0], \
295  &Q(0,0), &n, &Q(0,0), &n, &work1, &lwork, &info); \
296  lwork = BLAS_INT(gmm::real(work1)); \
297  std::vector<base_type> work(lwork); \
298  lapack_name(&jobvl, &jobvr, &n, &H(0,0), &n, &eigvr[0], &eigvi[0], \
299  &Q(0,0), &n, &Q(0,0), &n, &work[0], &lwork, &info); \
300  GMM_ASSERT1(!info, "QR algorithm failed"); \
301  gmm::copy(eigvr, gmm::real_part(eigval_)); \
302  gmm::copy(eigvi, gmm::imag_part(eigval_)); \
303  }
304 
305 # define geev_interface2(lapack_name, base_type, side) \
306  template <typename VECT> inline void geev_interface_ ## side( \
307  const dense_matrix<base_type> &A, VECT &eigval_, \
308  dense_matrix<base_type> &Q) { \
309  GMMLAPACK_TRACE("geev_interface"); \
310  BLAS_INT n = BLAS_INT(mat_nrows(A)), info(0), lwork(-1); \
311  base_type work1; \
312  if (!n) return; \
313  dense_matrix<base_type> H(n,n); gmm::copy(A, H); \
314  jobv_ ## side \
315  std::vector<base_type::value_type> rwork(2*n); \
316  std::vector<base_type> eigv(n); \
317  lapack_name(&jobvl, &jobvr, &n, &H(0,0), &n, &eigv[0], &Q(0,0), &n, \
318  &Q(0,0), &n, &work1, &lwork, &rwork[0], &info); \
319  lwork = BLAS_INT(gmm::real(work1)); \
320  std::vector<base_type> work(lwork); \
321  lapack_name(&jobvl, &jobvr, &n, &H(0,0), &n, &eigv[0], &Q(0,0), &n, \
322  &Q(0,0), &n, &work[0], &lwork, &rwork[0], &info); \
323  GMM_ASSERT1(!info, "QR algorithm failed"); \
324  gmm::copy(eigv, eigval_); \
325  }
326 
327  geev_interface(sgeev_, BLAS_S, right)
328  geev_interface(dgeev_, BLAS_D, right)
329  geev_interface2(cgeev_, BLAS_C, right)
330  geev_interface2(zgeev_, BLAS_Z, right)
331 
332  geev_interface(sgeev_, BLAS_S, left)
333  geev_interface(dgeev_, BLAS_D, left)
334  geev_interface2(cgeev_, BLAS_C, left)
335  geev_interface2(zgeev_, BLAS_Z, left)
336 
337 
338  /* ********************************************************************** */
339  /* SCHUR algorithm: */
340  /* A = Q*S*(Q^T), with Q orthogonal and S upper quasi-triangula */
341  /* ********************************************************************** */
342 
343 # define geesx_interface(lapack_name, base_type) inline \
344  void schur(dense_matrix<base_type> &A, \
345  dense_matrix<base_type> &S, \
346  dense_matrix<base_type> &Q) { \
347  GMMLAPACK_TRACE("geesx_interface"); \
348  BLAS_INT m = BLAS_INT(mat_nrows(A)), n = BLAS_INT(mat_ncols(A)); \
349  GMM_ASSERT1(m == n, "Schur decomposition requires square matrix"); \
350  char jobvs = 'V', sort = 'N', sense = 'N'; \
351  bool select = false; \
352  BLAS_INT lwork = 8*n, sdim = 0, liwork = 1; \
353  std::vector<base_type> work(lwork), wr(n), wi(n); \
354  std::vector<BLAS_INT> iwork(liwork); \
355  std::vector<BLAS_INT> bwork(1); \
356  resize(S, n, n); copy(A, S); \
357  resize(Q, n, n); \
358  base_type rconde(0), rcondv(0); \
359  BLAS_INT info(0); \
360  lapack_name(&jobvs, &sort, &select, &sense, &n, &S(0,0), &n, \
361  &sdim, &wr[0], &wi[0], &Q(0,0), &n, &rconde, &rcondv, \
362  &work[0], &lwork, &iwork[0], &liwork, &bwork[0], &info);\
363  GMM_ASSERT1(!info, "SCHUR algorithm failed"); \
364  }
365 
366 # define geesx_interface2(lapack_name, base_type) inline \
367  void schur(dense_matrix<base_type> &A, \
368  dense_matrix<base_type> &S, \
369  dense_matrix<base_type> &Q) { \
370  GMMLAPACK_TRACE("geesx_interface"); \
371  BLAS_INT m = BLAS_INT(mat_nrows(A)), n = BLAS_INT(mat_ncols(A)); \
372  GMM_ASSERT1(m == n, "Schur decomposition requires square matrix"); \
373  char jobvs = 'V', sort = 'N', sense = 'N'; \
374  bool select = false; \
375  BLAS_INT lwork = 8*n, sdim = 0; \
376  std::vector<base_type::value_type> rwork(lwork); \
377  std::vector<base_type> work(lwork), w(n); \
378  std::vector<BLAS_INT> bwork(1); \
379  resize(S, n, n); copy(A, S); \
380  resize(Q, n, n); \
381  base_type rconde(0), rcondv(0); \
382  BLAS_INT info(0); \
383  lapack_name(&jobvs, &sort, &select, &sense, &n, &S(0,0), &n, \
384  &sdim, &w[0], &Q(0,0), &n, &rconde, &rcondv, \
385  &work[0], &lwork, &rwork[0], &bwork[0], &info); \
386  GMM_ASSERT1(!info, "SCHUR algorithm failed"); \
387  }
388 
389  geesx_interface(sgeesx_, BLAS_S)
390  geesx_interface(dgeesx_, BLAS_D)
391  geesx_interface2(cgeesx_, BLAS_C)
392  geesx_interface2(zgeesx_, BLAS_Z)
393 
394  template <typename MAT>
395  void schur(const MAT &A_, MAT &S, MAT &Q) {
396  MAT A(A_);
397  schur(A, S, Q);
398  }
399 
400 
401  /* ********************************************************************** */
402  /* Interface to SVD. Does not correspond to a Gmm++ functionnality. */
403  /* Author : Sebastian Nowozin <[email protected]> */
404  /* ********************************************************************** */
405 
406 # define gesvd_interface(lapack_name, base_type) inline \
407  void svd(dense_matrix<base_type> &X, \
408  dense_matrix<base_type> &U, \
409  dense_matrix<base_type> &Vtransposed, \
410  std::vector<base_type> &sigma) { \
411  GMMLAPACK_TRACE("gesvd_interface"); \
412  BLAS_INT m = BLAS_INT(mat_nrows(X)), n = BLAS_INT(mat_ncols(X)); \
413  BLAS_INT mn_min = m < n ? m : n; \
414  sigma.resize(mn_min); \
415  std::vector<base_type> work(15 * mn_min); \
416  BLAS_INT lwork = BLAS_INT(work.size()); \
417  resize(U, m, m); \
418  resize(Vtransposed, n, n); \
419  char job = 'A'; \
420  BLAS_INT info(0); \
421  lapack_name(&job, &job, &m, &n, &X(0,0), &m, &sigma[0], &U(0,0), \
422  &m, &Vtransposed(0,0), &n, &work[0], &lwork, &info); \
423  }
424 
425 # define cgesvd_interface(lapack_name, base_type, base_type2) inline \
426  void svd(dense_matrix<base_type> &X, \
427  dense_matrix<base_type> &U, \
428  dense_matrix<base_type> &Vtransposed, \
429  std::vector<base_type2> &sigma) { \
430  GMMLAPACK_TRACE("gesvd_interface"); \
431  BLAS_INT m = BLAS_INT(mat_nrows(X)), n = BLAS_INT(mat_ncols(X)); \
432  BLAS_INT mn_min = m < n ? m : n; \
433  sigma.resize(mn_min); \
434  std::vector<base_type> work(15 * mn_min); \
435  std::vector<base_type2> rwork(5 * mn_min); \
436  BLAS_INT lwork = BLAS_INT(work.size()); \
437  resize(U, m, m); \
438  resize(Vtransposed, n, n); \
439  char job = 'A'; \
440  BLAS_INT info(0); \
441  lapack_name(&job, &job, &m, &n, &X(0,0), &m, &sigma[0], &U(0,0), \
442  &m, &Vtransposed(0,0), &n, &work[0], &lwork, \
443  &rwork[0], &info); \
444  }
445 
446  gesvd_interface(sgesvd_, BLAS_S)
447  gesvd_interface(dgesvd_, BLAS_D)
448  cgesvd_interface(cgesvd_, BLAS_C, BLAS_S)
449  cgesvd_interface(zgesvd_, BLAS_Z, BLAS_D)
450 
451  template <typename MAT, typename VEC>
452  void svd(const MAT &X_, MAT &U, MAT &Vtransposed, VEC &sigma) {
453  MAT X(X_);
454  svd(X, U, Vtransposed, sigma);
455  }
456 
457 }
458 
459 #else
460 
461 namespace gmm
462 {
463 template <typename MAT>
464 void schur(const MAT &, MAT &, MAT &)
465 {
466  GMM_ASSERT1(false, "Use of function schur(A,S,Q) requires GetFEM "
467  "to be built with Lapack");
468 }
469 
470 template <typename BLAS_TYPE>
471 inline void svd(dense_matrix<BLAS_TYPE> &, dense_matrix<BLAS_TYPE> &,
472  dense_matrix<BLAS_TYPE> &, std::vector<BLAS_TYPE> &)
473 {
474  GMM_ASSERT1(false, "Use of function svd(X,U,Vtransposed,sigma) requires GetFEM "
475  "to be built with Lapack");
476 }
477 
478 }// namespace gmm
479 
480 #endif // GMM_USES_LAPACK
481 
482 #endif // GMM_LAPACK_INTERFACE_H
gmm_blas_interface.h
gmm interface for fortran BLAS.
gmm_dense_lu.h
LU factorizations and determinant computation for dense matrices.
gmm_dense_qr.h
Dense QR factorization.