37 #ifndef GMM_LAPACK_INTERFACE_H
38 #define GMM_LAPACK_INTERFACE_H
44 #if defined(GMM_USES_LAPACK) && !defined(GMM_MATLAB_INTERFACE)
86 void sgetrf_(...);
void dgetrf_(...);
void cgetrf_(...);
void zgetrf_(...);
87 void sgetrs_(...);
void dgetrs_(...);
void cgetrs_(...);
void zgetrs_(...);
88 void sgetri_(...);
void dgetri_(...);
void cgetri_(...);
void zgetri_(...);
89 void sgeqrf_(...);
void dgeqrf_(...);
void cgeqrf_(...);
void zgeqrf_(...);
90 void sorgqr_(...);
void dorgqr_(...);
void cungqr_(...);
void zungqr_(...);
91 void sormqr_(...);
void dormqr_(...);
void cunmqr_(...);
void zunmqr_(...);
92 void sgees_ (...);
void dgees_ (...);
void cgees_ (...);
void zgees_ (...);
93 void sgeev_ (...);
void dgeev_ (...);
void cgeev_ (...);
void zgeev_ (...);
94 void sgeesx_(...);
void dgeesx_(...);
void cgeesx_(...);
void zgeesx_(...);
95 void sgesvd_(...);
void dgesvd_(...);
void cgesvd_(...);
void zgesvd_(...);
102 # define getrf_interface(lapack_name, base_type) \
104 lu_factor(dense_matrix<base_type> &A, lapack_ipvt &ipvt) { \
105 GMMLAPACK_TRACE("getrf_interface"); \
106 const BLAS_INT m=BLAS_INT(mat_nrows(A)), n=BLAS_INT(mat_ncols(A)), lda(m);\
108 if (m && n) lapack_name(&m, &n, &A(0,0), &lda, &ipvt[0], &info); \
109 return size_type(abs(info)); \
112 getrf_interface(sgetrf_, BLAS_S)
113 getrf_interface(dgetrf_, BLAS_D)
114 getrf_interface(cgetrf_, BLAS_C)
115 getrf_interface(zgetrf_, BLAS_Z)
121 # define getrs_interface(f_name, NorT, lapack_name, base_type) \
123 f_name(const dense_matrix<base_type> &A, const lapack_ipvt &ipvt, \
124 std::vector<base_type> &x, const std::vector<base_type> &b) { \
125 GMMLAPACK_TRACE("getrs_interface"); \
126 const char t=NorT; const BLAS_INT n=BLAS_INT(mat_nrows(A)), nrhs(1); \
127 BLAS_INT info(0); gmm::copy(b, x); \
128 if (n) lapack_name(&t, &n, &nrhs, &A(0,0), &n, &ipvt[0], &x[0], &n, \
132 getrs_interface(lu_solve,
'N', sgetrs_, BLAS_S)
133 getrs_interface(lu_solve,
'N', dgetrs_, BLAS_D)
134 getrs_interface(lu_solve,
'N', cgetrs_, BLAS_C)
135 getrs_interface(lu_solve,
'N', zgetrs_, BLAS_Z)
136 getrs_interface(lu_solve_transposed,
'T', sgetrs_, BLAS_S)
137 getrs_interface(lu_solve_transposed,
'T', dgetrs_, BLAS_D)
138 getrs_interface(lu_solve_transposed,
'T', cgetrs_, BLAS_C)
139 getrs_interface(lu_solve_transposed,
'T', zgetrs_, BLAS_Z)
145 # define getri_interface(lapack_name, base_type) \
146 inline void lu_inverse(const dense_matrix<base_type> &LU, \
147 const lapack_ipvt &ipvt, \
148 dense_matrix<base_type> &A) { \
149 GMMLAPACK_TRACE("getri_interface"); \
150 const BLAS_INT n=BLAS_INT(mat_nrows(A)); \
151 BLAS_INT info(0), lwork(-1); base_type work1; \
154 lapack_name(&n, &A(0,0), &n, &ipvt[0], &work1, &lwork, &info); \
155 const size_type worksize=size_type(gmm::real(work1)); \
156 std::vector<base_type> work(worksize); lwork=BLAS_INT(worksize); \
157 lapack_name(&n, &A(0,0), &n, &ipvt[0], &work[0], &lwork, &info); \
161 getri_interface(sgetri_, BLAS_S)
162 getri_interface(dgetri_, BLAS_D)
163 getri_interface(cgetri_, BLAS_C)
164 getri_interface(zgetri_, BLAS_Z)
170 # define geqrf_interface(lapack_name, base_type) \
171 inline void qr_factor(dense_matrix<base_type> &A) { \
172 GMMLAPACK_TRACE("geqrf_interface"); \
173 const size_type mm(mat_nrows(A)), nn(mat_ncols(A)); \
175 const BLAS_INT m=BLAS_INT(mm), n=BLAS_INT(nn); \
176 BLAS_INT info(0), lwork(-1); \
177 std::vector<base_type> tau(nn); base_type work1; \
178 lapack_name(&m, &n, &A(0,0), &m, &tau[0], &work1, &lwork, &info); \
179 const size_type worksize=size_type(gmm::real(work1)); \
180 std::vector<base_type> work(worksize); lwork=BLAS_INT(worksize); \
181 lapack_name(&m, &n, &A(0,0), &m, &tau[0], &work[0], &lwork, &info); \
182 GMM_ASSERT1(!info, "QR factorization failed"); \
186 geqrf_interface(sgeqrf_, BLAS_S)
187 geqrf_interface(dgeqrf_, BLAS_D)
193 # define geqrf_interface2(lapack_name1, lapack_name2, base_type) inline \
194 void qr_factor(const dense_matrix<base_type> &A, \
195 dense_matrix<base_type> &Q, dense_matrix<base_type> &R) { \
196 GMMLAPACK_TRACE("geqrf_interface2"); \
197 const size_type mm(mat_nrows(A)), nn(mat_ncols(A)); \
199 const BLAS_INT m=BLAS_INT(mm), n=BLAS_INT(nn); \
200 BLAS_INT info(0), lwork(-1); \
201 std::copy(A.begin(), A.end(), Q.begin()); \
202 std::vector<base_type> tau(nn); base_type work1; \
203 lapack_name1(&m, &n, &Q(0,0), &m, &tau[0], &work1, &lwork, &info); \
204 const size_type worksize=size_type(gmm::real(work1)); \
205 std::vector<base_type> work(worksize); lwork=BLAS_INT(worksize); \
206 lapack_name1(&m, &n, &Q(0,0), &m, &tau[0], &work[0], &lwork, &info); \
207 GMM_ASSERT1(!info, "QR factorization failed"); \
208 base_type *p = &R(0,0), *q = &Q(0,0); \
209 for (BLAS_INT j = 0; j < n; ++j, q += m-n) \
210 for (BLAS_INT i = 0; i < n; ++i, ++p, ++q) \
211 *p = (j < i) ? base_type(0) : *q; \
212 lapack_name2(&m, &n, &n, &Q(0,0), &m,&tau[0],&work[0],&lwork,&info); \
214 else gmm::clear(Q); \
217 geqrf_interface2(sgeqrf_, sorgqr_, BLAS_S)
218 geqrf_interface2(dgeqrf_, dorgqr_, BLAS_D)
219 geqrf_interface2(cgeqrf_, cungqr_, BLAS_C)
220 geqrf_interface2(zgeqrf_, zungqr_, BLAS_Z)
226 # define gees_interface(lapack_name, base_type) \
227 template <typename VECT> inline void implicit_qr_algorithm( \
228 const dense_matrix<base_type> &A, VECT &eigval_, \
229 dense_matrix<base_type> &Q, \
230 double tol=gmm::default_tol(base_type()), bool compvect = true) { \
231 GMMLAPACK_TRACE("gees_interface"); \
232 typedef bool (*L_fp)(...); L_fp p = 0; \
233 const size_type nn(mat_nrows(A)); \
235 dense_matrix<base_type> H(nn,nn); gmm::copy(A, H); \
236 const char jobvs=(compvect ? 'V' : 'N'), sort='N'; \
237 const BLAS_INT n=BLAS_INT(nn); \
238 std::vector<double> rwork(nn), eigv1(nn), eigv2(nn); \
239 BLAS_INT info(0), lwork(-1), sdim; base_type work1; \
240 lapack_name(&jobvs, &sort, p, &n, &H(0,0), &n, &sdim, &eigv1[0], \
241 &eigv2[0], &Q(0,0), &n, &work1, &lwork, &rwork[0], &info); \
242 const size_type worksize=size_type(gmm::real(work1)); \
243 std::vector<base_type> work(worksize); lwork=BLAS_INT(worksize); \
244 lapack_name(&jobvs, &sort, p, &n, &H(0,0), &n, &sdim, &eigv1[0], \
245 &eigv2[0], &Q(0,0), &n, &work[0], &lwork, &rwork[0],&info);\
246 GMM_ASSERT1(!info, "QR algorithm failed"); \
247 extract_eig(H, eigval_, tol); \
250 # define gees_interface_cplx(lapack_name, base_type) \
251 template <typename VECT> inline void implicit_qr_algorithm( \
252 const dense_matrix<base_type> &A, VECT &eigval_, \
253 dense_matrix<base_type> &Q, \
254 double tol=gmm::default_tol(base_type()), bool compvect = true) { \
255 GMMLAPACK_TRACE("gees_interface2"); \
256 typedef bool (*L_fp)(...); L_fp p = 0; \
257 const size_type nn(mat_nrows(A)); \
259 dense_matrix<base_type> H(nn,nn); gmm::copy(A, H); \
260 const char jobvs=(compvect ? 'V' : 'N'), sort='N'; \
261 const BLAS_INT n=BLAS_INT(nn); \
262 std::vector<double> rwork(nn), eigvv(nn*2); \
263 BLAS_INT info(0), lwork(-1), sdim; base_type work1; \
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 const size_type worksize=size_type(gmm::real(work1)); \
267 std::vector<base_type> work(worksize); lwork=BLAS_INT(worksize); \
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); \
274 gees_interface(sgees_, BLAS_S)
275 gees_interface(dgees_, BLAS_D)
276 gees_interface_cplx(cgees_, BLAS_C)
277 gees_interface_cplx(zgees_, BLAS_Z)
280 # define geev_interface(f_name, lapack_name, lapack_name_cplx, \
281 base_type, cplx_type, _jobvl, _jobvr) \
282 template <typename VECT> inline void \
283 f_name(const dense_matrix<base_type> &A, VECT &eigval_, \
284 dense_matrix<base_type> &Q) { \
285 GMMLAPACK_TRACE("geev_interface"); \
286 const size_type nn(mat_nrows(A)); \
288 dense_matrix<base_type> H(nn,nn); gmm::copy(A, H); \
289 const char jobvl=_jobvl, jobvr=_jobvl; \
290 const BLAS_INT n=BLAS_INT(nn); BLAS_INT info(0), lwork(-1); \
291 std::vector<base_type> eigvr(nn), eigvi(nn); \
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 const size_type worksize=size_type(gmm::real(work1)); \
296 std::vector<base_type> work(worksize); lwork=BLAS_INT(worksize); \
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(eigval_)); \
301 gmm::copy(eigvi, gmm::imag_part(eigval_)); \
303 template <typename VECT> inline void \
304 f_name(const dense_matrix<cplx_type> &A, VECT &eigval_, \
305 dense_matrix<cplx_type> &Q) { \
306 GMMLAPACK_TRACE("geev_interface"); \
307 const size_type nn(mat_nrows(A)); \
309 dense_matrix<cplx_type> H(nn,nn); gmm::copy(A, H); \
310 const char jobvl=_jobvl, jobvr=_jobvl; \
311 const BLAS_INT n=BLAS_INT(nn); BLAS_INT info(0), lwork(-1); \
312 std::vector<cplx_type> eigv(nn); std::vector<base_type> rwork(2*nn); \
314 lapack_name(&jobvl, &jobvr, &n, &H(0,0), &n, &eigv[0], &Q(0,0), &n, \
315 &Q(0,0), &n, &work1, &lwork, &rwork[0], &info); \
316 const size_type worksize=size_type(gmm::real(work1)); \
317 std::vector<cplx_type> work(worksize); lwork=BLAS_INT(worksize); \
318 lapack_name_cplx(&jobvl, &jobvr, &n, &H(0,0), &n, &eigv[0], &Q(0,0), \
319 &n, &Q(0,0), &n, &work[0], &lwork, &rwork[0], &info);\
320 GMM_ASSERT1(!info, "QR algorithm failed"); \
321 gmm::copy(eigv, eigval_); \
324 geev_interface(geev_interface_right, sgeev_, cgeev_, BLAS_S, BLAS_C,
'N',
'V')
325 geev_interface(geev_interface_right, dgeev_, zgeev_, BLAS_D, BLAS_Z, 'N', 'V')
326 geev_interface(geev_interface_left, sgeev_, cgeev_, BLAS_S, BLAS_C, 'V', 'N')
327 geev_interface(geev_interface_left, dgeev_, zgeev_, BLAS_D, BLAS_Z, 'V', 'N')
334 # define geesx_interface(lapack_name, lapack_name_cplx, \
335 base_type, cplx_type) \
336 inline void schur(const dense_matrix<base_type> &A, \
337 dense_matrix<base_type> &S, \
338 dense_matrix<base_type> &Q) { \
339 GMMLAPACK_TRACE("geesx_interface"); \
340 const size_type mm(mat_nrows(A)), nn(mat_ncols(A)), worksize(8*nn); \
341 GMM_ASSERT1(mm==nn, "Schur decomposition requires square matrix"); \
342 resize(S, nn, nn); copy(A, S); resize(Q, nn, nn); \
343 const char jobvs='V', sort='N', sense='N'; const bool select=false; \
344 const BLAS_INT n=BLAS_INT(nn), lwork=BLAS_INT(worksize), liwork(1); \
345 BLAS_INT info(0), sdim(0); base_type rconde(0), rcondv(0); \
346 std::vector<base_type> work(worksize), wr(nn), wi(nn); \
347 std::vector<BLAS_INT> iwork(1), bwork(1); \
348 lapack_name(&jobvs, &sort, &select, &sense, &n, &S(0,0), &n, \
349 &sdim, &wr[0], &wi[0], &Q(0,0), &n, &rconde, &rcondv, \
350 &work[0], &lwork, &iwork[0], &liwork, &bwork[0], &info);\
351 GMM_ASSERT1(!info, "SCHUR algorithm failed"); \
353 inline void schur(const dense_matrix<cplx_type> &A, \
354 dense_matrix<cplx_type> &S, \
355 dense_matrix<cplx_type> &Q) { \
356 GMMLAPACK_TRACE("geesx_interface"); \
357 const size_type mm(mat_nrows(A)), nn(mat_ncols(A)), worksize(8*nn); \
358 GMM_ASSERT1(mm==nn, "Schur decomposition requires square matrix"); \
359 resize(S, nn, nn); copy(A, S); resize(Q, nn, nn); \
360 const char jobvs='V', sort='N', sense='N'; const bool select=false; \
361 const BLAS_INT n=BLAS_INT(nn), lwork=BLAS_INT(worksize); \
362 BLAS_INT info(0), sdim(0); base_type rconde(0), rcondv(0); \
363 std::vector<cplx_type> work(worksize), w(nn); \
364 std::vector<base_type> rwork(worksize); \
365 std::vector<BLAS_INT> bwork(1); \
366 lapack_name_cplx(&jobvs, &sort, &select, &sense, &n, &S(0,0), &n, \
367 &sdim, &w[0], &Q(0,0), &n, &rconde, &rcondv, \
368 &work[0], &lwork, &rwork[0], &bwork[0], &info); \
369 GMM_ASSERT1(!info, "SCHUR algorithm failed"); \
372 geesx_interface(sgeesx_, cgeesx_, BLAS_S, BLAS_C)
373 geesx_interface(dgeesx_, zgeesx_, BLAS_D, BLAS_Z)
381 # define gesvd_interface(lapack_name, lapack_name_cplx, \
382 base_type, cplx_type) \
383 inline void svd(const dense_matrix<base_type> &X, \
384 dense_matrix<base_type> &U, \
385 dense_matrix<base_type> &Vtransposed, \
386 std::vector<base_type> &sigma) { \
387 GMMLAPACK_TRACE("gesvd_interface"); \
388 const size_type mm(mat_nrows(X)), nn(mat_ncols(X)), \
389 mn_min(mm < nn ? mm : nn), worksize(15*mn_min); \
390 sigma.resize(mn_min); resize(U, mm, mm); resize(Vtransposed, nn, nn); \
391 const char job='A'; const BLAS_INT m=BLAS_INT(mm), n=BLAS_INT(nn), \
392 lwork=BLAS_INT(worksize); \
393 std::vector<base_type> work(worksize); BLAS_INT info(0); \
394 lapack_name(&job, &job, &m, &n, &X(0,0), &m, &sigma[0], &U(0,0), \
395 &m, &Vtransposed(0,0), &n, &work[0], &lwork, &info); \
397 inline void svd(const dense_matrix<cplx_type> &X, \
398 dense_matrix<cplx_type> &U, \
399 dense_matrix<cplx_type> &Vtransposed, \
400 std::vector<base_type> &sigma) { \
401 GMMLAPACK_TRACE("gesvd_interface"); \
402 const size_type mm(mat_nrows(X)), nn(mat_ncols(X)), \
403 mn_min(mm < nn ? mm : nn), worksize(15*mn_min); \
404 sigma.resize(mn_min); resize(U, mm, mm); resize(Vtransposed, nn, nn); \
405 const char job='A'; const BLAS_INT m=BLAS_INT(mm), n=BLAS_INT(nn), \
406 lwork=BLAS_INT(worksize); \
407 std::vector<cplx_type> work(worksize); \
408 std::vector<base_type> rwork(5*mn_min); BLAS_INT info(0); \
409 lapack_name_cplx(&job, &job, &m, &n, &X(0,0), &m, &sigma[0], &U(0,0), \
410 &m, &Vtransposed(0,0), &n, &work[0], &lwork, \
414 gesvd_interface(sgesvd_, cgesvd_, BLAS_S, BLAS_C)
415 gesvd_interface(dgesvd_, zgesvd_, BLAS_D, BLAS_Z)
423 template <
typename MAT>
424 void schur(
const MAT &, MAT &, MAT &)
426 GMM_ASSERT1(
false,
"Use of function schur(A,S,Q) requires GetFEM "
427 "to be built with Lapack");
430 template <
typename BLAS_TYPE>
431 inline void svd(dense_matrix<BLAS_TYPE> &, dense_matrix<BLAS_TYPE> &,
432 dense_matrix<BLAS_TYPE> &, std::vector<BLAS_TYPE> &)
434 GMM_ASSERT1(
false,
"Use of function svd(X,U,Vtransposed,sigma) requires GetFEM "
435 "to be built with Lapack");
gmm interface for fortran BLAS.
LU factorizations and determinant computation for dense matrices.