GetFEM++  5.3
gmm_lapack_interface.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2003-2017 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  long m = long(mat_nrows(A)), n = long(mat_ncols(A)), lda(m), info(-1L); \
108  if (m && n) lapack_name(&m, &n, &A(0,0), &lda, ipvt.pfirst(), &info); \
109  if ((info & 0xFFFFFFFF00000000L) && !(info & 0x00000000FFFFFFFFL)) \
110  /* For compatibility with lapack version with 32 bit integer. */ \
111  ipvt.set_to_int32(); \
112  return size_type(int(info & 0x00000000FFFFFFFFL)); \
113  }
114 
115  getrf_interface(sgetrf_, BLAS_S)
116  getrf_interface(dgetrf_, BLAS_D)
117  getrf_interface(cgetrf_, BLAS_C)
118  getrf_interface(zgetrf_, BLAS_Z)
119 
120  /* ********************************************************************* */
121  /* LU solve. */
122  /* ********************************************************************* */
123 
124 # define getrs_interface(f_name, trans1, lapack_name, base_type) inline \
125  void f_name(const dense_matrix<base_type > &A, \
126  const lapack_ipvt &ipvt, std::vector<base_type > &x, \
127  const std::vector<base_type > &b) { \
128  GMMLAPACK_TRACE("getrs_interface"); \
129  long n = long(mat_nrows(A)), info(0), nrhs(1); \
130  gmm::copy(b, x); trans1; \
131  if (n) \
132  lapack_name(&t,&n,&nrhs,&(A(0,0)),&n,ipvt.pfirst(),&x[0],&n,&info); \
133  }
134 
135 # define getrs_trans_n const char t = 'N'
136 # define getrs_trans_t const char t = 'T'
137 
138  getrs_interface(lu_solve, getrs_trans_n, sgetrs_, BLAS_S)
139  getrs_interface(lu_solve, getrs_trans_n, dgetrs_, BLAS_D)
140  getrs_interface(lu_solve, getrs_trans_n, cgetrs_, BLAS_C)
141  getrs_interface(lu_solve, getrs_trans_n, zgetrs_, BLAS_Z)
142  getrs_interface(lu_solve_transposed, getrs_trans_t, sgetrs_, BLAS_S)
143  getrs_interface(lu_solve_transposed, getrs_trans_t, dgetrs_, BLAS_D)
144  getrs_interface(lu_solve_transposed, getrs_trans_t, cgetrs_, BLAS_C)
145  getrs_interface(lu_solve_transposed, getrs_trans_t, zgetrs_, BLAS_Z)
146 
147  /* ********************************************************************* */
148  /* LU inverse. */
149  /* ********************************************************************* */
150 
151 # define getri_interface(lapack_name, base_type) inline \
152  void lu_inverse(const dense_matrix<base_type > &LU, \
153  const lapack_ipvt &ipvt, const dense_matrix<base_type > &A_) { \
154  GMMLAPACK_TRACE("getri_interface"); \
155  dense_matrix<base_type >& \
156  A = const_cast<dense_matrix<base_type > &>(A_); \
157  long n = int(mat_nrows(A)), info(0), lwork(-1); base_type work1; \
158  if (n) { \
159  gmm::copy(LU, A); \
160  lapack_name(&n, &A(0,0), &n, ipvt.pfirst(), &work1, &lwork, &info); \
161  lwork = int(gmm::real(work1)); \
162  std::vector<base_type > work(lwork); \
163  lapack_name(&n, &A(0,0), &n, ipvt.pfirst(), &work[0], &lwork,&info); \
164  } \
165  }
166 
167  getri_interface(sgetri_, BLAS_S)
168  getri_interface(dgetri_, BLAS_D)
169  getri_interface(cgetri_, BLAS_C)
170  getri_interface(zgetri_, BLAS_Z)
171 
172  /* ********************************************************************** */
173  /* QR factorization. */
174  /* ********************************************************************** */
175 
176 # define geqrf_interface(lapack_name1, base_type) inline \
177  void qr_factor(dense_matrix<base_type > &A){ \
178  GMMLAPACK_TRACE("geqrf_interface"); \
179  long m = long(mat_nrows(A)), n=long(mat_ncols(A)), info(0), lwork(-1); \
180  base_type work1; \
181  if (m && n) { \
182  std::vector<base_type > tau(n); \
183  lapack_name1(&m, &n, &A(0,0), &m, &tau[0], &work1 , &lwork, &info); \
184  lwork = long(gmm::real(work1)); \
185  std::vector<base_type > work(lwork); \
186  lapack_name1(&m, &n, &A(0,0), &m, &tau[0], &work[0], &lwork, &info); \
187  GMM_ASSERT1(!info, "QR factorization failed"); \
188  } \
189  }
190 
191  geqrf_interface(sgeqrf_, BLAS_S)
192  geqrf_interface(dgeqrf_, BLAS_D)
193  // For complex values, housholder vectors are not the same as in
194  // gmm::lu_factor. Impossible to interface for the moment.
195  // geqrf_interface(cgeqrf_, BLAS_C)
196  // geqrf_interface(zgeqrf_, BLAS_Z)
197 
198 # define geqrf_interface2(lapack_name1, lapack_name2, base_type) inline \
199  void qr_factor(const dense_matrix<base_type > &A, \
200  dense_matrix<base_type > &Q, dense_matrix<base_type > &R) { \
201  GMMLAPACK_TRACE("geqrf_interface2"); \
202  long m = long(mat_nrows(A)), n=long(mat_ncols(A)), info(0), lwork(-1); \
203  base_type work1; \
204  if (m && n) { \
205  std::copy(A.begin(), A.end(), Q.begin()); \
206  std::vector<base_type > tau(n); \
207  lapack_name1(&m, &n, &Q(0,0), &m, &tau[0], &work1 , &lwork, &info); \
208  lwork = long(gmm::real(work1)); \
209  std::vector<base_type > work(lwork); \
210  lapack_name1(&m, &n, &Q(0,0), &m, &tau[0], &work[0], &lwork, &info); \
211  GMM_ASSERT1(!info, "QR factorization failed"); \
212  base_type *p = &R(0,0), *q = &Q(0,0); \
213  for (long j = 0; j < n; ++j, q += m-n) \
214  for (long i = 0; i < n; ++i, ++p, ++q) \
215  *p = (j < i) ? base_type(0) : *q; \
216  lapack_name2(&m, &n, &n, &Q(0,0), &m,&tau[0],&work[0],&lwork,&info); \
217  } \
218  else gmm::clear(Q); \
219  }
220 
221  geqrf_interface2(sgeqrf_, sorgqr_, BLAS_S)
222  geqrf_interface2(dgeqrf_, dorgqr_, BLAS_D)
223  geqrf_interface2(cgeqrf_, cungqr_, BLAS_C)
224  geqrf_interface2(zgeqrf_, zungqr_, BLAS_Z)
225 
226  /* ********************************************************************** */
227  /* QR algorithm for eigenvalues search. */
228  /* ********************************************************************** */
229 
230 # define gees_interface(lapack_name, base_type) \
231  template <typename VECT> inline void implicit_qr_algorithm( \
232  const dense_matrix<base_type > &A, const VECT &eigval_, \
233  dense_matrix<base_type > &Q, \
234  double tol=gmm::default_tol(base_type()), bool compvect = true) { \
235  GMMLAPACK_TRACE("gees_interface"); \
236  typedef bool (*L_fp)(...); L_fp p = 0; \
237  long n=long(mat_nrows(A)), info(0), lwork(-1), sdim; base_type work1; \
238  if (!n) return; \
239  dense_matrix<base_type > H(n,n); gmm::copy(A, H); \
240  char jobvs = (compvect ? 'V' : 'N'), sort = 'N'; \
241  std::vector<double> rwork(n), eigv1(n), eigv2(n); \
242  lapack_name(&jobvs, &sort, p, &n, &H(0,0), &n, &sdim, &eigv1[0], \
243  &eigv2[0], &Q(0,0), &n, &work1, &lwork, &rwork[0], &info); \
244  lwork = long(gmm::real(work1)); \
245  std::vector<base_type > work(lwork); \
246  lapack_name(&jobvs, &sort, p, &n, &H(0,0), &n, &sdim, &eigv1[0], \
247  &eigv2[0], &Q(0,0), &n, &work[0], &lwork, &rwork[0],&info);\
248  GMM_ASSERT1(!info, "QR algorithm failed"); \
249  extract_eig(H, const_cast<VECT &>(eigval_), tol); \
250  }
251 
252 # define gees_interface2(lapack_name, base_type) \
253  template <typename VECT> inline void implicit_qr_algorithm( \
254  const dense_matrix<base_type > &A, const VECT &eigval_, \
255  dense_matrix<base_type > &Q, \
256  double tol=gmm::default_tol(base_type()), bool compvect = true) { \
257  GMMLAPACK_TRACE("gees_interface2"); \
258  typedef bool (*L_fp)(...); L_fp p = 0; \
259  long n=long(mat_nrows(A)), info(0), lwork(-1), sdim; 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 = long(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, const_cast<VECT &>(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, const VECT &eigval_, \
286  dense_matrix<base_type > &Q) { \
287  GMMLAPACK_TRACE("geev_interface"); \
288  long n = long(mat_nrows(A)), info(0), lwork(-1); base_type work1; \
289  if (!n) return; \
290  dense_matrix<base_type > H(n,n); gmm::copy(A, H); \
291  jobv_ ## side \
292  std::vector<base_type > eigvr(n), eigvi(n); \
293  lapack_name(&jobvl, &jobvr, &n, &H(0,0), &n, &eigvr[0], &eigvi[0], \
294  &Q(0,0), &n, &Q(0,0), &n, &work1, &lwork, &info); \
295  lwork = long(gmm::real(work1)); \
296  std::vector<base_type > work(lwork); \
297  lapack_name(&jobvl, &jobvr, &n, &H(0,0), &n, &eigvr[0], &eigvi[0], \
298  &Q(0,0), &n, &Q(0,0), &n, &work[0], &lwork, &info); \
299  GMM_ASSERT1(!info, "QR algorithm failed"); \
300  gmm::copy(eigvr, gmm::real_part(const_cast<VECT &>(eigval_))); \
301  gmm::copy(eigvi, gmm::imag_part(const_cast<VECT &>(eigval_))); \
302  }
303 
304 # define geev_interface2(lapack_name, base_type, side) \
305  template <typename VECT> inline void geev_interface_ ## side( \
306  const dense_matrix<base_type > &A, const VECT &eigval_, \
307  dense_matrix<base_type > &Q) { \
308  GMMLAPACK_TRACE("geev_interface"); \
309  long n = long(mat_nrows(A)), info(0), lwork(-1); base_type work1; \
310  if (!n) return; \
311  dense_matrix<base_type > H(n,n); gmm::copy(A, H); \
312  jobv_ ## side \
313  std::vector<base_type::value_type> rwork(2*n); \
314  std::vector<base_type> eigv(n); \
315  lapack_name(&jobvl, &jobvr, &n, &H(0,0), &n, &eigv[0], &Q(0,0), &n, \
316  &Q(0,0), &n, &work1, &lwork, &rwork[0], &info); \
317  lwork = long(gmm::real(work1)); \
318  std::vector<base_type > work(lwork); \
319  lapack_name(&jobvl, &jobvr, &n, &H(0,0), &n, &eigv[0], &Q(0,0), &n, \
320  &Q(0,0), &n, &work[0], &lwork, &rwork[0], &info); \
321  GMM_ASSERT1(!info, "QR algorithm failed"); \
322  gmm::copy(eigv, const_cast<VECT &>(eigval_)); \
323  }
324 
325  geev_interface(sgeev_, BLAS_S, right)
326  geev_interface(dgeev_, BLAS_D, right)
327  geev_interface2(cgeev_, BLAS_C, right)
328  geev_interface2(zgeev_, BLAS_Z, right)
329 
330  geev_interface(sgeev_, BLAS_S, left)
331  geev_interface(dgeev_, BLAS_D, left)
332  geev_interface2(cgeev_, BLAS_C, left)
333  geev_interface2(zgeev_, BLAS_Z, left)
334 
335 
336  /* ********************************************************************** */
337  /* SCHUR algorithm: */
338  /* A = Q*S*(Q^T), with Q orthogonal and S upper quasi-triangula */
339  /* ********************************************************************** */
340 
341 # define geesx_interface(lapack_name, base_type) inline \
342  void schur(dense_matrix<base_type> &A, \
343  dense_matrix<base_type> &S, \
344  dense_matrix<base_type> &Q) { \
345  GMMLAPACK_TRACE("geesx_interface"); \
346  long m = long(mat_nrows(A)), n = long(mat_ncols(A)); \
347  GMM_ASSERT1(m == n, "Schur decomposition requires square matrix"); \
348  char jobvs = 'V', sort = 'N', sense = 'N'; \
349  bool select = false; \
350  long lwork = 8*n, sdim = 0, liwork = 1; \
351  std::vector<base_type> work(lwork), wr(n), wi(n); \
352  std::vector<long> iwork(liwork); \
353  std::vector<long> bwork(1); \
354  resize(S, n, n); copy(A, S); \
355  resize(Q, n, n); \
356  base_type rconde(0), rcondv(0); \
357  long info(0); \
358  lapack_name(&jobvs, &sort, &select, &sense, &n, &S(0,0), &n, \
359  &sdim, &wr[0], &wi[0], &Q(0,0), &n, &rconde, &rcondv, \
360  &work[0], &lwork, &iwork[0], &liwork, &bwork[0], &info);\
361  GMM_ASSERT1(!info, "SCHUR algorithm failed"); \
362  }
363 
364 # define geesx_interface2(lapack_name, base_type) inline \
365  void schur(dense_matrix<base_type> &A, \
366  dense_matrix<base_type> &S, \
367  dense_matrix<base_type> &Q) { \
368  GMMLAPACK_TRACE("geesx_interface"); \
369  long m = long(mat_nrows(A)), n = long(mat_ncols(A)); \
370  GMM_ASSERT1(m == n, "Schur decomposition requires square matrix"); \
371  char jobvs = 'V', sort = 'N', sense = 'N'; \
372  bool select = false; \
373  long lwork = 8*n, sdim = 0; \
374  std::vector<base_type::value_type> rwork(lwork); \
375  std::vector<base_type> work(lwork), w(n); \
376  std::vector<long> bwork(1); \
377  resize(S, n, n); copy(A, S); \
378  resize(Q, n, n); \
379  base_type rconde(0), rcondv(0); \
380  long info(0); \
381  lapack_name(&jobvs, &sort, &select, &sense, &n, &S(0,0), &n, \
382  &sdim, &w[0], &Q(0,0), &n, &rconde, &rcondv, \
383  &work[0], &lwork, &rwork[0], &bwork[0], &info); \
384  GMM_ASSERT1(!info, "SCHUR algorithm failed"); \
385  }
386 
387  geesx_interface(sgeesx_, BLAS_S)
388  geesx_interface(dgeesx_, BLAS_D)
389  geesx_interface2(cgeesx_, BLAS_C)
390  geesx_interface2(zgeesx_, BLAS_Z)
391 
392  template <typename MAT>
393  void schur(const MAT &A_, MAT &S, MAT &Q) {
394  MAT A(A_);
395  schur(A, S, Q);
396  }
397 
398 
399  /* ********************************************************************** */
400  /* Interface to SVD. Does not correspond to a Gmm++ functionnality. */
401  /* Author : Sebastian Nowozin <sebastian.nowozin@tuebingen.mpg.de> */
402  /* ********************************************************************** */
403 
404 # define gesvd_interface(lapack_name, base_type) inline \
405  void svd(dense_matrix<base_type> &X, \
406  dense_matrix<base_type> &U, \
407  dense_matrix<base_type> &Vtransposed, \
408  std::vector<base_type> &sigma) { \
409  GMMLAPACK_TRACE("gesvd_interface"); \
410  long m = long(mat_nrows(X)), n = long(mat_ncols(X)); \
411  long mn_min = m < n ? m : n; \
412  sigma.resize(mn_min); \
413  std::vector<base_type> work(15 * mn_min); \
414  long lwork = long(work.size()); \
415  resize(U, m, m); \
416  resize(Vtransposed, n, n); \
417  char job = 'A'; \
418  long info(0); \
419  lapack_name(&job, &job, &m, &n, &X(0,0), &m, &sigma[0], &U(0,0), \
420  &m, &Vtransposed(0,0), &n, &work[0], &lwork, &info); \
421  }
422 
423 # define cgesvd_interface(lapack_name, base_type, base_type2) inline \
424  void svd(dense_matrix<base_type> &X, \
425  dense_matrix<base_type> &U, \
426  dense_matrix<base_type> &Vtransposed, \
427  std::vector<base_type2> &sigma) { \
428  GMMLAPACK_TRACE("gesvd_interface"); \
429  long m = long(mat_nrows(X)), n = long(mat_ncols(X)); \
430  long mn_min = m < n ? m : n; \
431  sigma.resize(mn_min); \
432  std::vector<base_type> work(15 * mn_min); \
433  std::vector<base_type2> rwork(5 * mn_min); \
434  long lwork = long(work.size()); \
435  resize(U, m, m); \
436  resize(Vtransposed, n, n); \
437  char job = 'A'; \
438  long info(0); \
439  lapack_name(&job, &job, &m, &n, &X(0,0), &m, &sigma[0], &U(0,0), \
440  &m, &Vtransposed(0,0), &n, &work[0], &lwork, \
441  &rwork[0], &info); \
442  }
443 
444  gesvd_interface(sgesvd_, BLAS_S)
445  gesvd_interface(dgesvd_, BLAS_D)
446  cgesvd_interface(cgesvd_, BLAS_C, BLAS_S)
447  cgesvd_interface(zgesvd_, BLAS_Z, BLAS_D)
448 
449  template <typename MAT, typename VEC>
450  void svd(const MAT &X_, MAT &U, MAT &Vtransposed, VEC &sigma) {
451  MAT X(X_);
452  svd(X, U, Vtransposed, sigma);
453  }
454 
455 
456 
457 
458 }
459 
460 #else
461 
462 namespace gmm
463 {
464 template <typename MAT>
465 void schur(const MAT &A_, MAT &S, MAT &Q)
466 {
467  GMM_ASSERT1(false, "Use of function schur(A,S,Q) requires GetFEM++ "
468  "to be built with Lapack");
469 }
470 
471 }// namespace gmm
472 
473 #endif // GMM_USES_LAPACK
474 
475 #endif // GMM_LAPACK_INTERFACE_H
gmm_dense_qr.h
Dense QR factorization.
gmm_dense_lu.h
LU factorizations and determinant computation for dense matrices.
gmm_blas_interface.h
gmm interface for fortran BLAS.