Nektar++
Lapack.hpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: Lapack.hpp
4//
5// For more information, please see: http://www.nektar.info
6//
7// The MIT License
8//
9// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10// Department of Aeronautics, Imperial College London (UK), and Scientific
11// Computing and Imaging Institute, University of Utah (USA).
12//
13// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description: wrapper of functions around standard LAPACK routines
32//
33///////////////////////////////////////////////////////////////////////////////
34#ifndef NEKTAR_LIB_UTILITIES_LAPACK_HPP
35#define NEKTAR_LIB_UTILITIES_LAPACK_HPP
36
38
39namespace Lapack
40{
41extern "C"
42{
43 // Matrix factorisation and solves
44 void F77NAME(dsptrf)(const char &uplo, const int &n, double *ap, int *ipiv,
45 int &info);
46 void F77NAME(ssptrf)(const char &uplo, const int &n, float *ap, int *ipiv,
47 int &info);
48 void F77NAME(dsptrs)(const char &uplo, const int &n, const int &nrhs,
49 const double *ap, const int *ipiv, double *b,
50 const int &ldb, int &info);
51 void F77NAME(ssptrs)(const char &uplo, const int &n, const int &nrhs,
52 const float *ap, const int *ipiv, float *b,
53 const int &ldb, int &info);
54 void F77NAME(dsptri)(const char &uplo, const int &n, const double *ap,
55 const int *ipiv, double *work, int &info);
56 void F77NAME(ssptri)(const char &uplo, const int &n, const float *ap,
57 const int *ipiv, float *work, int &info);
58 void F77NAME(dtrtrs)(const char &uplo, const char &trans, const char &diag,
59 const int &n, const int &nrhs, const double *a,
60 const int &lda, double *b, const int &ldb, int &info);
61 void F77NAME(dtptrs)(const char &uplo, const char &trans, const char &diag,
62 const int &n, const int &nrhs, const double *a,
63 double *b, const int &ldb, int &info);
64 void F77NAME(dpptrf)(const char &uplo, const int &n, double *ap, int &info);
65 void F77NAME(dpptrs)(const char &uplo, const int &n, const int &nrhs,
66 const double *ap, double *b, const int &ldb,
67 int &info);
68 void F77NAME(dpbtrf)(const char &uplo, const int &n, const int &kd,
69 double *ab, const int &ldab, int &info);
70 void F77NAME(dpbtrs)(const char &uplo, const int &n, const int &kd,
71 const int &nrhs, const double *ab, const int &ldab,
72 double *b, const int &ldb, int &info);
73 void F77NAME(dgbtrf)(const int &m, const int &n, const int &kl,
74 const int &ku, double *a, const int &lda, int *ipiv,
75 int &info);
76 void F77NAME(dgbtrs)(const char &trans, const int &n, const int &kl,
77 const int &ku, const int &nrhs, const double *a,
78 const int &lda, const int *ipiv, double *b,
79 const int &ldb, int &info);
80 void F77NAME(dgetrf)(const int &m, const int &n, double *a, const int &lda,
81 int *ipiv, int &info);
82 void F77NAME(sgetrf)(const int &m, const int &n, float *a, const int &lda,
83 int *ipiv, int &info);
84 void F77NAME(dgetrs)(const char &trans, const int &n, const int &nrhs,
85 const double *a, const int &lda, int *ipiv, double *b,
86 const int &ldb, int &info);
87 void F77NAME(dgetri)(const int &n, double *a, const int &lda,
88 const int *ipiv, double *wk, const int &lwk,
89 int &info);
90 void F77NAME(sgetri)(const int &n, float *a, const int &lda,
91 const int *ipiv, float *wk, const int &lwk, int &info);
92 void F77NAME(dsterf)(const int &n, double *d, double *e, int &info);
93 void F77NAME(dgeev)(const char &uplo, const char &lrev, const int &n,
94 const double *a, const int &lda, double *wr, double *wi,
95 double *rev, const int &ldr, double *lev,
96 const int &ldv, double *work, const int &lwork,
97 int &info);
98 void F77NAME(sgeev)(const char &uplo, const char &lrev, const int &n,
99 const float *a, const int &lda, float *wr, float *wi,
100 float *rev, const int &ldr, float *lev, const int &ldv,
101 float *work, const int &lwork, int &info);
102 void F77NAME(dspev)(const char &jobz, const char &uplo, const int &n,
103 double *ap, double *w, double *z, const int &ldz,
104 double *work, int &info);
105 void F77NAME(dsbev)(const char &jobz, const char &uplo, const int &kl,
106 const int &ku, double *ap, const int &lda, double *w,
107 double *z, const int &ldz, double *work, int &info);
108}
109
110/// \brief factor a real packed-symmetric matrix using Bunch-Kaufman
111/// pivoting.
112static inline void DoSsptrf(const char &uplo, const int &n, double *ap,
113 int *ipiv, int &info)
114{
115 F77NAME(dsptrf)(uplo, n, ap, ipiv, info);
116}
117
118/// \brief factor a real packed-symmetric matrix using Bunch-Kaufman
119/// pivoting.
120static inline void DoSsptrf(const char &uplo, const int &n, float *ap,
121 int *ipiv, int &info)
122{
123 F77NAME(ssptrf)(uplo, n, ap, ipiv, info);
124}
125
126/// \brief factor a real packed-symmetric matrix using Bunch-Kaufman
127/// pivoting.
128static inline void Dsptrf(const char &uplo, const int &n, double *ap, int *ipiv,
129 int &info)
130{
131 F77NAME(dsptrf)(uplo, n, ap, ipiv, info);
132}
133
134/// \brief Solve a real packed-symmetric matrix problem using Bunch-Kaufman
135/// pivoting.
136static inline void DoSsptrs(const char &uplo, const int &n, const int &nrhs,
137 const double *ap, const int *ipiv, double *b,
138 const int &ldb, int &info)
139{
140 F77NAME(dsptrs)(uplo, n, nrhs, ap, ipiv, b, ldb, info);
141}
142
143/// \brief Solve a real packed-symmetric matrix problem using Bunch-Kaufman
144/// pivoting.
145static inline void DoSsptrs(const char &uplo, const int &n, const int &nrhs,
146 const float *ap, const int *ipiv, float *b,
147 const int &ldb, int &info)
148{
149 F77NAME(ssptrs)(uplo, n, nrhs, ap, ipiv, b, ldb, info);
150}
151
152/// \brief Solve a real packed-symmetric matrix problem using Bunch-Kaufman
153/// pivoting.
154static inline void Dsptrs(const char &uplo, const int &n, const int &nrhs,
155 const double *ap, const int *ipiv, double *b,
156 const int &ldb, int &info)
157{
158 F77NAME(dsptrs)(uplo, n, nrhs, ap, ipiv, b, ldb, info);
159}
160
161/// \brief Invert a real packed-symmetric matrix problem
162static inline void DoSsptri(const char &uplo, const int &n, const double *ap,
163 const int *ipiv, double *work, int &info)
164{
165 F77NAME(dsptri)(uplo, n, ap, ipiv, work, info);
166}
167
168/// \brief Invert a real packed-symmetric matrix problem
169static inline void DoSsptri(const char &uplo, const int &n, const float *ap,
170 const int *ipiv, float *work, int &info)
171{
172 F77NAME(ssptri)(uplo, n, ap, ipiv, work, info);
173}
174
175/// \brief Solve a triangular system.
176static inline void Dtrtrs(const char &uplo, const char &trans, const char &diag,
177 const int &n, const int &nrhs, const double *a,
178 const int &lda, double *b, const int &ldb, int &info)
179{
180 F77NAME(dtrtrs)(uplo, trans, diag, n, nrhs, a, lda, b, ldb, info);
181}
182
183/// \brief Solve a triangular system.
184static inline void Dtptrs(const char &uplo, const char &trans, const char &diag,
185 const int &n, const int &nrhs, const double *a,
186 double *b, const int &ldb, int &info)
187{
188 F77NAME(dtptrs)(uplo, trans, diag, n, nrhs, a, b, ldb, info);
189}
190
191/// \brief Invert a real packed-symmetric matrix problem
192static inline void Dsptri(const char &uplo, const int &n, const double *ap,
193 const int *ipiv, double *work, int &info)
194{
195 F77NAME(dsptri)(uplo, n, ap, ipiv, work, info);
196}
197
198/// \brief Cholesky factor a real positive definite packed-symmetric matrix.
199static inline void Dpptrf(const char &uplo, const int &n, double *ap, int &info)
200{
201 F77NAME(dpptrf)(uplo, n, ap, info);
202}
203
204/// \brief Solve a real positive definite symmetric matrix problem
205/// using Cholesky factorization.
206static inline void Dpptrs(const char &uplo, const int &n, const int &nrhs,
207 const double *ap, double *b, const int &ldb,
208 int &info)
209{
210 F77NAME(dpptrs)(uplo, n, nrhs, ap, b, ldb, info);
211}
212
213/// \brief Cholesky factorize a real positive definite
214/// banded-symmetric matrix
215static inline void Dpbtrf(const char &uplo, const int &n, const int &kd,
216 double *ab, const int &ldab, int &info)
217{
218 F77NAME(dpbtrf)(uplo, n, kd, ab, ldab, info);
219}
220
221/// \brief Solve a real, positive definite banded-symmetric matrix
222/// problem using Cholesky factorization.
223static inline void Dpbtrs(const char &uplo, const int &n, const int &kd,
224 const int &nrhs, const double *ab, const int &ldab,
225 double *b, const int &ldb, int &info)
226{
227 F77NAME(dpbtrs)(uplo, n, kd, nrhs, ab, ldab, b, ldb, info);
228}
229
230/// \brief General banded matrix LU factorisation
231static inline void Dgbtrf(const int &m, const int &n, const int &kl,
232 const int &ku, double *a, const int &lda, int *ipiv,
233 int &info)
234{
235 F77NAME(dgbtrf)(m, n, kl, ku, a, lda, ipiv, info);
236}
237
238/// \brief Solve general banded matrix using LU factorisation
239static inline void Dgbtrs(const char &trans, const int &n, const int &kl,
240 const int &ku, const int &nrhs, const double *a,
241 const int &lda, const int *ipiv, double *b,
242 const int &ldb, int &info)
243{
244 F77NAME(dgbtrs)(trans, n, kl, ku, nrhs, a, lda, ipiv, b, ldb, info);
245}
246
247/// \brief General matrix LU factorisation
248static inline void DoSgetrf(const int &m, const int &n, double *a,
249 const int &lda, int *ipiv, int &info)
250{
251 F77NAME(dgetrf)(m, n, a, lda, ipiv, info);
252}
253
254/// \brief General matrix LU factorisation
255static inline void DoSgetrf(const int &m, const int &n, float *a,
256 const int &lda, int *ipiv, int &info)
257{
258 F77NAME(sgetrf)(m, n, a, lda, ipiv, info);
259}
260
261/// \brief General matrix LU factorisation
262static inline void Dgetrf(const int &m, const int &n, double *a, const int &lda,
263 int *ipiv, int &info)
264{
265 F77NAME(dgetrf)(m, n, a, lda, ipiv, info);
266}
267
268/// \brief General matrix LU backsolve
269static inline void Dgetrs(const char &trans, const int &n, const int &nrhs,
270 const double *a, const int &lda, int *ipiv, double *b,
271 const int &ldb, int &info)
272{
273 F77NAME(dgetrs)(trans, n, nrhs, a, lda, ipiv, b, ldb, info);
274}
275
276/// \brief General matrix inverse
277static inline void DoSgetri(const int &n, double *a, const int &lda,
278 const int *ipiv, double *wk, const int &lwk,
279 int &info)
280{
281 F77NAME(dgetri)(n, a, lda, ipiv, wk, lwk, info);
282}
283
284/// \brief General matrix inverse
285static inline void DoSgetri(const int &n, float *a, const int &lda,
286 const int *ipiv, float *wk, const int &lwk,
287 int &info)
288{
289 F77NAME(sgetri)(n, a, lda, ipiv, wk, lwk, info);
290}
291
292/// \brief General matrix inverse
293static inline void Dgetri(const int &n, double *a, const int &lda,
294 const int *ipiv, double *wk, const int &lwk,
295 int &info)
296{
297 F77NAME(dgetri)(n, a, lda, ipiv, wk, lwk, info);
298}
299
300/// \brief Find eigenvalues of symmetric tridiagonal matrix
301static inline void Dsterf(const int &n, double *d, double *e, int &info)
302{
303 F77NAME(dsterf)(n, d, e, info);
304}
305
306/// \brief Solve general real matrix eigenproblem.
307static inline void DoSgeev(const char &uplo, const char &lrev, const int &n,
308 const double *a, const int &lda, double *wr,
309 double *wi, double *rev, const int &ldr, double *lev,
310 const int &ldv, double *work, const int &lwork,
311 int &info)
312{
314 (uplo, lrev, n, a, lda, wr, wi, rev, ldr, lev, ldv, work, lwork, info);
315}
316
317/// \brief Solve general real matrix eigenproblem.
318static inline void DoSgeev(const char &uplo, const char &lrev, const int &n,
319 const float *a, const int &lda, float *wr, float *wi,
320 float *rev, const int &ldr, float *lev,
321 const int &ldv, float *work, const int &lwork,
322 int &info)
323{
325 (uplo, lrev, n, a, lda, wr, wi, rev, ldr, lev, ldv, work, lwork, info);
326}
327
328/// \brief Solve general real matrix eigenproblem.
329static inline void Dgeev(const char &uplo, const char &lrev, const int &n,
330 const double *a, const int &lda, double *wr,
331 double *wi, double *rev, const int &ldr, double *lev,
332 const int &ldv, double *work, const int &lwork,
333 int &info)
334{
336 (uplo, lrev, n, a, lda, wr, wi, rev, ldr, lev, ldv, work, lwork, info);
337}
338
339/// \brief Solve packed-symmetric real matrix eigenproblem.
340static inline void Dspev(const char &jobz, const char &uplo, const int &n,
341 double *ap, double *w, double *z, const int &ldz,
342 double *work, int &info)
343{
344 F77NAME(dspev)(jobz, uplo, n, ap, w, z, ldz, work, info);
345}
346
347/// \brief Solve packed-banded real matrix eigenproblem.
348static inline void Dsbev(const char &jobz, const char &uplo, const int &kl,
349 const int &ku, double *ap, const int &lda, double *w,
350 double *z, const int &ldz, double *work, int &info)
351{
352 F77NAME(dsbev)(jobz, uplo, kl, ku, ap, lda, w, z, ldz, work, info);
353}
354
355} // namespace Lapack
356#endif // NEKTAR_LIB_UTILITIES_LAPACK_HPP
#define F77NAME(x)
Fortran routines need an underscore.
Definition: TransF77.hpp:46
void F77NAME() dsbev(const char &jobz, const char &uplo, const int &kl, const int &ku, double *ap, const int &lda, double *w, double *z, const int &ldz, double *work, int &info)
void F77NAME() dpbtrf(const char &uplo, const int &n, const int &kd, double *ab, const int &ldab, int &info)
static void Dspev(const char &jobz, const char &uplo, const int &n, double *ap, double *w, double *z, const int &ldz, double *work, int &info)
Solve packed-symmetric real matrix eigenproblem.
Definition: Lapack.hpp:340
void F77NAME() dpptrf(const char &uplo, const int &n, double *ap, int &info)
void F77NAME() ssptrf(const char &uplo, const int &n, float *ap, int *ipiv, int &info)
static void Dpptrs(const char &uplo, const int &n, const int &nrhs, const double *ap, double *b, const int &ldb, int &info)
Solve a real positive definite symmetric matrix problem using Cholesky factorization.
Definition: Lapack.hpp:206
static void Dgeev(const char &uplo, const char &lrev, const int &n, const double *a, const int &lda, double *wr, double *wi, double *rev, const int &ldr, double *lev, const int &ldv, double *work, const int &lwork, int &info)
Solve general real matrix eigenproblem.
Definition: Lapack.hpp:329
void F77NAME() dgetrs(const char &trans, const int &n, const int &nrhs, const double *a, const int &lda, int *ipiv, double *b, const int &ldb, int &info)
void F77NAME() dgetri(const int &n, double *a, const int &lda, const int *ipiv, double *wk, const int &lwk, int &info)
static void Dgbtrf(const int &m, const int &n, const int &kl, const int &ku, double *a, const int &lda, int *ipiv, int &info)
General banded matrix LU factorisation.
Definition: Lapack.hpp:231
void F77NAME() dgetrf(const int &m, const int &n, double *a, const int &lda, int *ipiv, int &info)
static void Dsptrs(const char &uplo, const int &n, const int &nrhs, const double *ap, const int *ipiv, double *b, const int &ldb, int &info)
Solve a real packed-symmetric matrix problem using Bunch-Kaufman pivoting.
Definition: Lapack.hpp:154
static void Dsbev(const char &jobz, const char &uplo, const int &kl, const int &ku, double *ap, const int &lda, double *w, double *z, const int &ldz, double *work, int &info)
Solve packed-banded real matrix eigenproblem.
Definition: Lapack.hpp:348
static void Dtrtrs(const char &uplo, const char &trans, const char &diag, const int &n, const int &nrhs, const double *a, const int &lda, double *b, const int &ldb, int &info)
Solve a triangular system.
Definition: Lapack.hpp:176
static void DoSsptrf(const char &uplo, const int &n, double *ap, int *ipiv, int &info)
factor a real packed-symmetric matrix using Bunch-Kaufman pivoting.
Definition: Lapack.hpp:112
void F77NAME() sgeev(const char &uplo, const char &lrev, const int &n, const float *a, const int &lda, float *wr, float *wi, float *rev, const int &ldr, float *lev, const int &ldv, float *work, const int &lwork, int &info)
static void Dgetrf(const int &m, const int &n, double *a, const int &lda, int *ipiv, int &info)
General matrix LU factorisation.
Definition: Lapack.hpp:262
void F77NAME() dgeev(const char &uplo, const char &lrev, const int &n, const double *a, const int &lda, double *wr, double *wi, double *rev, const int &ldr, double *lev, const int &ldv, double *work, const int &lwork, int &info)
static void Dpptrf(const char &uplo, const int &n, double *ap, int &info)
Cholesky factor a real positive definite packed-symmetric matrix.
Definition: Lapack.hpp:199
static void DoSsptri(const char &uplo, const int &n, const double *ap, const int *ipiv, double *work, int &info)
Invert a real packed-symmetric matrix problem.
Definition: Lapack.hpp:162
static void Dsptri(const char &uplo, const int &n, const double *ap, const int *ipiv, double *work, int &info)
Invert a real packed-symmetric matrix problem.
Definition: Lapack.hpp:192
static void Dgetrs(const char &trans, const int &n, const int &nrhs, const double *a, const int &lda, int *ipiv, double *b, const int &ldb, int &info)
General matrix LU backsolve.
Definition: Lapack.hpp:269
void F77NAME() dsterf(const int &n, double *d, double *e, int &info)
static void Dsptrf(const char &uplo, const int &n, double *ap, int *ipiv, int &info)
factor a real packed-symmetric matrix using Bunch-Kaufman pivoting.
Definition: Lapack.hpp:128
static void Dpbtrf(const char &uplo, const int &n, const int &kd, double *ab, const int &ldab, int &info)
Cholesky factorize a real positive definite banded-symmetric matrix.
Definition: Lapack.hpp:215
void F77NAME() sgetri(const int &n, float *a, const int &lda, const int *ipiv, float *wk, const int &lwk, int &info)
void F77NAME() dsptri(const char &uplo, const int &n, const double *ap, const int *ipiv, double *work, int &info)
void F77NAME() dspev(const char &jobz, const char &uplo, const int &n, double *ap, double *w, double *z, const int &ldz, double *work, int &info)
static void Dgbtrs(const char &trans, const int &n, const int &kl, const int &ku, const int &nrhs, const double *a, const int &lda, const int *ipiv, double *b, const int &ldb, int &info)
Solve general banded matrix using LU factorisation.
Definition: Lapack.hpp:239
static void DoSgetrf(const int &m, const int &n, double *a, const int &lda, int *ipiv, int &info)
General matrix LU factorisation.
Definition: Lapack.hpp:248
static void Dgetri(const int &n, double *a, const int &lda, const int *ipiv, double *wk, const int &lwk, int &info)
General matrix inverse.
Definition: Lapack.hpp:293
static void DoSgetri(const int &n, double *a, const int &lda, const int *ipiv, double *wk, const int &lwk, int &info)
General matrix inverse.
Definition: Lapack.hpp:277
void F77NAME() dgbtrs(const char &trans, const int &n, const int &kl, const int &ku, const int &nrhs, const double *a, const int &lda, const int *ipiv, double *b, const int &ldb, int &info)
static void Dsterf(const int &n, double *d, double *e, int &info)
Find eigenvalues of symmetric tridiagonal matrix.
Definition: Lapack.hpp:301
void F77NAME() dsptrf(const char &uplo, const int &n, double *ap, int *ipiv, int &info)
void F77NAME() dpptrs(const char &uplo, const int &n, const int &nrhs, const double *ap, double *b, const int &ldb, int &info)
static void DoSgeev(const char &uplo, const char &lrev, const int &n, const double *a, const int &lda, double *wr, double *wi, double *rev, const int &ldr, double *lev, const int &ldv, double *work, const int &lwork, int &info)
Solve general real matrix eigenproblem.
Definition: Lapack.hpp:307
void F77NAME() sgetrf(const int &m, const int &n, float *a, const int &lda, int *ipiv, int &info)
void F77NAME() ssptrs(const char &uplo, const int &n, const int &nrhs, const float *ap, const int *ipiv, float *b, const int &ldb, int &info)
void F77NAME() ssptri(const char &uplo, const int &n, const float *ap, const int *ipiv, float *work, int &info)
static void Dtptrs(const char &uplo, const char &trans, const char &diag, const int &n, const int &nrhs, const double *a, double *b, const int &ldb, int &info)
Solve a triangular system.
Definition: Lapack.hpp:184
static void DoSsptrs(const char &uplo, const int &n, const int &nrhs, const double *ap, const int *ipiv, double *b, const int &ldb, int &info)
Solve a real packed-symmetric matrix problem using Bunch-Kaufman pivoting.
Definition: Lapack.hpp:136
void F77NAME() dsptrs(const char &uplo, const int &n, const int &nrhs, const double *ap, const int *ipiv, double *b, const int &ldb, int &info)
void F77NAME() dtrtrs(const char &uplo, const char &trans, const char &diag, const int &n, const int &nrhs, const double *a, const int &lda, double *b, const int &ldb, int &info)
void F77NAME() dtptrs(const char &uplo, const char &trans, const char &diag, const int &n, const int &nrhs, const double *a, double *b, const int &ldb, int &info)
void F77NAME() dpbtrs(const char &uplo, const int &n, const int &kd, const int &nrhs, const double *ab, const int &ldab, double *b, const int &ldb, int &info)
void F77NAME() dgbtrf(const int &m, const int &n, const int &kl, const int &ku, double *a, const int &lda, int *ipiv, int &info)
static void Dpbtrs(const char &uplo, const int &n, const int &kd, const int &nrhs, const double *ab, const int &ldab, double *b, const int &ldb, int &info)
Solve a real, positive definite banded-symmetric matrix problem using Cholesky factorization.
Definition: Lapack.hpp:223
std::vector< double > w(NPUPPER)
std::vector< double > z(NPUPPER)
std::vector< double > d(NPUPPER *NPUPPER)