mxlib
c++ tools for analyzing astronomical data and other tasks by Jared R. Males. [git repo]
Loading...
Searching...
No Matches
templateLapack.cpp
Go to the documentation of this file.
1/** \file templateLapack.cpp
2 * \brief Implementation of templatized wrappers for the Lapack library
3 * \ingroup gen_math_files
4 *
5 */
6
7//***********************************************************************//
8// Copyright 2020 Jared R. Males (jaredmales@gmail.com)
9//
10// This file is part of mxlib.
11//
12// mxlib is free software: you can redistribute it and/or modify
13// it under the terms of the GNU General Public License as published by
14// the Free Software Foundation, either version 3 of the License, or
15// (at your option) any later version.
16//
17// mxlib is distributed in the hope that it will be useful,
18// but WITHOUT ANY WARRANTY; without even the implied warranty of
19// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20// GNU General Public License for more details.
21//
22// You should have received a copy of the GNU General Public License
23// along with mxlib. If not, see <http://www.gnu.org/licenses/>.
24//***********************************************************************//
25
26#include "math/templateLapack.hpp"
27
28#ifndef MXLIB_MKL
29// Generic LAPACKE headers do not declare the provider's LAED9 Fortran symbols.
30extern "C"
31{
32 void slaed9_( const MXLAPACK_INT *K,
33 const MXLAPACK_INT *KSTART,
34 const MXLAPACK_INT *KSTOP,
35 const MXLAPACK_INT *N,
36 float *D,
37 float *Q,
38 const MXLAPACK_INT *LDQ,
39 const float *RHO,
40 float *DLAMDA,
41 float *W,
42 float *S,
43 const MXLAPACK_INT *LDS,
44 MXLAPACK_INT *INFO );
45
46 void dlaed9_( const MXLAPACK_INT *K,
47 const MXLAPACK_INT *KSTART,
48 const MXLAPACK_INT *KSTOP,
49 const MXLAPACK_INT *N,
50 double *D,
51 double *Q,
52 const MXLAPACK_INT *LDQ,
53 const double *RHO,
54 double *DLAMDA,
55 double *W,
56 double *S,
57 const MXLAPACK_INT *LDS,
58 MXLAPACK_INT *INFO );
59}
60#endif
61
62namespace mx
63{
64namespace math
65{
66
67template <>
68float lamch<float>( char CMACH )
69{
70 return slamch_( &CMACH
71#ifdef LAPACK_FORTRAN_STRLEN_END
72 ,
73 1
74#endif
75 );
76}
77
78// Double specialization of lamch, a wrapper for Lapack DLAMCH
79template <>
80double lamch<double>( char CMACH )
81{
82 return dlamch_( &CMACH
83#ifdef LAPACK_FORTRAN_STRLEN_END
84 ,
85 1
86#endif
87 );
88}
89
90template <>
91MXLAPACK_INT potrf<float>( char UPLO, MXLAPACK_INT N, float *A, MXLAPACK_INT LDA, MXLAPACK_INT &INFO )
92{
93 spotrf_( &UPLO,
94 &N,
95 A,
96 &LDA,
97 &INFO
98#ifdef LAPACK_FORTRAN_STRLEN_END
99 ,
100 1
101#endif
102 );
103
104 return INFO;
105}
106
107template <>
108MXLAPACK_INT potrf<double>( char UPLO, MXLAPACK_INT N, double *A, MXLAPACK_INT LDA, MXLAPACK_INT &INFO )
109{
110 dpotrf_( &UPLO,
111 &N,
112 A,
113 &LDA,
114 &INFO
115#ifdef LAPACK_FORTRAN_STRLEN_END
116 ,
117 1
118#endif
119 );
120
121 return INFO;
122}
123
124template <>
125MXLAPACK_INT
126potrf<std::complex<float>>( char UPLO, MXLAPACK_INT N, std::complex<float> *A, MXLAPACK_INT LDA, MXLAPACK_INT &INFO )
127{
128 cpotrf_( &UPLO,
129 &N,
130 (float _Complex *)A,
131 &LDA,
132 &INFO
133#ifdef LAPACK_FORTRAN_STRLEN_END
134 ,
135 1
136#endif
137 );
138
139 return INFO;
140}
141
142template <>
143MXLAPACK_INT
144potrf<std::complex<double>>( char UPLO, MXLAPACK_INT N, std::complex<double> *A, MXLAPACK_INT LDA, MXLAPACK_INT &INFO )
145{
146 zpotrf_( &UPLO,
147 &N,
148 (double _Complex *)A,
149 &LDA,
150 &INFO
151#ifdef LAPACK_FORTRAN_STRLEN_END
152 ,
153 1
154#endif
155 );
156
157 return INFO;
158}
159
160template <>
161MXLAPACK_INT sytrd<float>( char UPLO,
162 MXLAPACK_INT N,
163 float *A,
164 MXLAPACK_INT LDA,
165 float *D,
166 float *E,
167 float *TAU,
168 float *WORK,
169 MXLAPACK_INT LWORK,
170 MXLAPACK_INT INFO )
171{
172
173 ssytrd_( &UPLO,
174 &N,
175 A,
176 &LDA,
177 D,
178 E,
179 TAU,
180 WORK,
181 &LWORK,
182 &INFO
183#ifdef LAPACK_FORTRAN_STRLEN_END
184 ,
185 1
186#endif
187 );
188
189 return INFO;
190}
191
192template <>
193MXLAPACK_INT sytrd<double>( char UPLO,
194 MXLAPACK_INT N,
195 double *A,
196 MXLAPACK_INT LDA,
197 double *D,
198 double *E,
199 double *TAU,
200 double *WORK,
201 MXLAPACK_INT LWORK,
202 MXLAPACK_INT INFO )
203{
204
205 dsytrd_( &UPLO,
206 &N,
207 A,
208 &LDA,
209 D,
210 E,
211 TAU,
212 WORK,
213 &LWORK,
214 &INFO
215#ifdef LAPACK_FORTRAN_STRLEN_END
216 ,
217 1
218#endif
219 );
220
221 return INFO;
222}
223
224template <>
225MXLAPACK_INT laed9<float>( float *D,
226 float *Q,
227 float *S,
228 MXLAPACK_INT K,
229 MXLAPACK_INT KSTART,
230 MXLAPACK_INT KSTOP,
231 MXLAPACK_INT N,
232 MXLAPACK_INT LDQ,
233 float RHO,
234 float *DLAMDA,
235 float *W,
236 MXLAPACK_INT LDS )
237{
238 MXLAPACK_INT INFO;
239
240 slaed9_( &K, &KSTART, &KSTOP, &N, D, Q, &LDQ, &RHO, DLAMDA, W, S, &LDS, &INFO );
241
242 return INFO;
243}
244
245template <>
246MXLAPACK_INT laed9<double>( double *D,
247 double *Q,
248 double *S,
249 MXLAPACK_INT K,
250 MXLAPACK_INT KSTART,
251 MXLAPACK_INT KSTOP,
252 MXLAPACK_INT N,
253 MXLAPACK_INT LDQ,
254 double RHO,
255 double *DLAMDA,
256 double *W,
257 MXLAPACK_INT LDS )
258{
259 MXLAPACK_INT INFO;
260
261 dlaed9_( &K, &KSTART, &KSTOP, &N, D, Q, &LDQ, &RHO, DLAMDA, W, S, &LDS, &INFO );
262
263 return INFO;
264}
265
266// Float specialization of syevr, a wrapper for Lapack SSYEVR
267template <>
268MXLAPACK_INT syevr<float>( char JOBZ,
269 char RANGE,
270 char UPLO,
271 MXLAPACK_INT N,
272 float *A,
273 MXLAPACK_INT LDA,
274 float VL,
275 float VU,
276 MXLAPACK_INT IL,
277 MXLAPACK_INT IU,
278 float ABSTOL,
279 MXLAPACK_INT *M,
280 float *W,
281 float *Z,
282 MXLAPACK_INT LDZ,
283 MXLAPACK_INT *ISUPPZ,
284 float *WORK,
285 MXLAPACK_INT LWORK,
286 MXLAPACK_INT *IWORK,
287 MXLAPACK_INT LIWORK )
288{
289
290 MXLAPACK_INT INFO;
291
292 ssyevr_( &JOBZ,
293 &RANGE,
294 &UPLO,
295 &N,
296 A,
297 &LDA,
298 &VL,
299 &VU,
300 &IL,
301 &IU,
302 &ABSTOL,
303 M,
304 W,
305 Z,
306 &LDZ,
307 ISUPPZ,
308 WORK,
309 &LWORK,
310 IWORK,
311 &LIWORK,
312 &INFO
313#ifdef LAPACK_FORTRAN_STRLEN_END
314 ,
315 1,
316 1,
317 1
318#endif
319 );
320
321 return INFO;
322}
323
324// Double specialization of syevr, a wrapper for Lapack DSYEVR
325template <>
326MXLAPACK_INT syevr<double>( char JOBZ,
327 char RANGE,
328 char UPLO,
329 MXLAPACK_INT N,
330 double *A,
331 MXLAPACK_INT LDA,
332 double VL,
333 double VU,
334 MXLAPACK_INT IL,
335 MXLAPACK_INT IU,
336 double ABSTOL,
337 MXLAPACK_INT *M,
338 double *W,
339 double *Z,
340 MXLAPACK_INT LDZ,
341 MXLAPACK_INT *ISUPPZ,
342 double *WORK,
343 MXLAPACK_INT LWORK,
344 MXLAPACK_INT *IWORK,
345 MXLAPACK_INT LIWORK )
346{
347
348 MXLAPACK_INT INFO;
349
350 dsyevr_( &JOBZ,
351 &RANGE,
352 &UPLO,
353 &N,
354 A,
355 &LDA,
356 &VL,
357 &VU,
358 &IL,
359 &IU,
360 &ABSTOL,
361 M,
362 W,
363 Z,
364 &LDZ,
365 ISUPPZ,
366 WORK,
367 &LWORK,
368 IWORK,
369 &LIWORK,
370 &INFO
371#ifdef LAPACK_FORTRAN_STRLEN_END
372 ,
373 1,
374 1,
375 1
376#endif
377 );
378
379 return INFO;
380}
381
382// float specialization of gesvd
383template <>
384MXLAPACK_INT gesvd<float>( char JOBU,
385 char JOBVT,
386 MXLAPACK_INT M,
387 MXLAPACK_INT N,
388 float *A,
389 MXLAPACK_INT LDA,
390 float *S,
391 float *U,
392 MXLAPACK_INT LDU,
393 float *VT,
394 MXLAPACK_INT LDVT,
395 float *WORK,
396 MXLAPACK_INT LWORK )
397{
398 MXLAPACK_INT INFO;
399
400 sgesvd_( &JOBU,
401 &JOBVT,
402 &M,
403 &N,
404 A,
405 &LDA,
406 S,
407 U,
408 &LDU,
409 VT,
410 &LDVT,
411 WORK,
412 &LWORK,
413 &INFO
414#ifdef LAPACK_FORTRAN_STRLEN_END
415 ,
416 1,
417 1
418#endif
419 );
420
421 return INFO;
422}
423
424// double specialization of gesvd
425template <>
426MXLAPACK_INT gesvd<double>( char JOBU,
427 char JOBVT,
428 MXLAPACK_INT M,
429 MXLAPACK_INT N,
430 double *A,
431 MXLAPACK_INT LDA,
432 double *S,
433 double *U,
434 MXLAPACK_INT LDU,
435 double *VT,
436 MXLAPACK_INT LDVT,
437 double *WORK,
438 MXLAPACK_INT LWORK )
439{
440 MXLAPACK_INT INFO;
441
442 dgesvd_( &JOBU,
443 &JOBVT,
444 &M,
445 &N,
446 A,
447 &LDA,
448 S,
449 U,
450 &LDU,
451 VT,
452 &LDVT,
453 WORK,
454 &LWORK,
455 &INFO
456#ifdef LAPACK_FORTRAN_STRLEN_END
457 ,
458 1,
459 1
460#endif
461 );
462
463 return INFO;
464}
465
466// float specialization of gesdd
467template <>
468MXLAPACK_INT gesdd<float>( char JOBZ,
469 MXLAPACK_INT M,
470 MXLAPACK_INT N,
471 float *A,
472 MXLAPACK_INT LDA,
473 float *S,
474 float *U,
475 MXLAPACK_INT LDU,
476 float *VT,
477 MXLAPACK_INT LDVT,
478 float *WORK,
479 MXLAPACK_INT LWORK,
480 MXLAPACK_INT *IWORK )
481{
482 MXLAPACK_INT INFO;
483
484 sgesdd_( &JOBZ,
485 &M,
486 &N,
487 A,
488 &LDA,
489 S,
490 U,
491 &LDU,
492 VT,
493 &LDVT,
494 WORK,
495 &LWORK,
496 IWORK,
497 &INFO
498#ifdef LAPACK_FORTRAN_STRLEN_END
499 ,
500 1
501#endif
502 );
503
504 return INFO;
505}
506
507// double specialization of gesdd
508template <>
509MXLAPACK_INT gesdd<double>( char JOBZ,
510 MXLAPACK_INT M,
511 MXLAPACK_INT N,
512 double *A,
513 MXLAPACK_INT LDA,
514 double *S,
515 double *U,
516 MXLAPACK_INT LDU,
517 double *VT,
518 MXLAPACK_INT LDVT,
519 double *WORK,
520 MXLAPACK_INT LWORK,
521 MXLAPACK_INT *IWORK )
522{
523
524 MXLAPACK_INT INFO;
525
526 dgesdd_( &JOBZ,
527 &M,
528 &N,
529 A,
530 &LDA,
531 S,
532 U,
533 &LDU,
534 VT,
535 &LDVT,
536 WORK,
537 &LWORK,
538 IWORK,
539 &INFO
540#ifdef LAPACK_FORTRAN_STRLEN_END
541 ,
542 1
543#endif
544 );
545
546 return INFO;
547}
548
549} // namespace math
550} // namespace mx
mx::math::sytrd
MXLAPACK_INT sytrd(char UPLO, MXLAPACK_INT N, dataT *A, MXLAPACK_INT LDA, dataT *D, dataT *E, dataT *TAU, dataT *WORK, MXLAPACK_INT LWORK, MXLAPACK_INT INFO)
Reduce a real symmetric matrix to real symmetric tridiagonal form by an orthogonal similarity transfo...
Definition templateLapack.hpp:183
mx::math::gesdd
MXLAPACK_INT gesdd(char JOBZ, MXLAPACK_INT M, MXLAPACK_INT N, dataT *A, MXLAPACK_INT LDA, dataT *S, dataT *U, MXLAPACK_INT LDU, dataT *VT, MXLAPACK_INT LDVT, dataT *WORK, MXLAPACK_INT LWORK, MXLAPACK_INT *IWORK)
Compute the singular value decomposition (SVD) of a real matrix with GESDD.
Definition templateLapack.hpp:718
mx::math::lamch
dataT lamch(char CMACH)
Determine machine parameters.
Definition templateLapack.hpp:104
mx::math::gesvd
MXLAPACK_INT gesvd(char JOBU, char JOBVT, MXLAPACK_INT M, MXLAPACK_INT N, dataT *A, MXLAPACK_INT LDA, dataT *S, dataT *U, MXLAPACK_INT LDU, dataT *VT, MXLAPACK_INT LDVT, dataT *WORK, MXLAPACK_INT LWORK)
Compute the singular value decomposition (SVD) of a real matrix.
Definition templateLapack.hpp:524
mx::math::laed9
MXLAPACK_INT laed9(dataT *D, dataT *Q, dataT *S, MXLAPACK_INT K, MXLAPACK_INT KSTART, MXLAPACK_INT KSTOP, MXLAPACK_INT N, MXLAPACK_INT LDQ, dataT RHO, dataT *DLAMDA, dataT *W, MXLAPACK_INT LDS)
Solve selected roots of a diagonal-plus-rank-one secular equation and form its eigenvectors.
Definition templateLapack.hpp:245
mx
The mxlib c++ namespace.
Definition mxlib.hpp:37
templateLapack.hpp
Declares and defines templatized wrappers for the Lapack library.
mx::math::potrf
MXLAPACK_INT potrf(char UPLO, MXLAPACK_INT N, dataT *A, MXLAPACK_INT LDA, MXLAPACK_INT &INFO)
Compute the Cholesky factorization of a real symmetric positive definite matrix A.