Loading [MathJax]/extensions/tex2jax.js
mxlib
c++ tools for analyzing astronomical data and other tasks by Jared R. Males. [git repo]
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Modules Pages
templateLapack.hpp
Go to the documentation of this file.
1/** \file templateLapack.hpp
2 * \brief Declares and defines templatized wrappers for the Lapack library
3 * \ingroup gen_math_files
4 * \author Jared R. Males (jaredmales@gmail.com)
5 *
6 */
7
8//***********************************************************************//
9// Copyright 2015, 2016, 2017 Jared R. Males (jaredmales@gmail.com)
10//
11// This file is part of mxlib.
12//
13// mxlib is free software: you can redistribute it and/or modify
14// it under the terms of the GNU General Public License as published by
15// the Free Software Foundation, either version 3 of the License, or
16// (at your option) any later version.
17//
18// mxlib is distributed in the hope that it will be useful,
19// but WITHOUT ANY WARRANTY; without even the implied warranty of
20// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21// GNU General Public License for more details.
22//
23// You should have received a copy of the GNU General Public License
24// along with mxlib. If not, see <http://www.gnu.org/licenses/>.
25//***********************************************************************//
26
27#ifndef math_templateLapack_hpp
28#define math_templateLapack_hpp
29
30#include <complex>
31
32// If using MKL:
33extern "C"
34{
35#if defined( MXLIB_MKL )
36
37#define MKL_Complex8 float _Complex
38#define MKL_Complex16 double _Complex
39
40#include <mkl.h>
41
42#endif
43}
44
45// MKL can use 64-bit integer types, standard BLAS and LAPACK use int
46// MKL declares some pointer args const. Compiling with ATLAS doesn't seem to care, but just to be safe...
47#ifdef MXLIB_MKL
48typedef MKL_INT MXLAPACK_INT;
49#define MKL_CONST_PTR const
50#else
51#ifndef MXLAPACK_INT
52typedef int MXLAPACK_INT;
53#endif
54
55#define MKL_CONST_PTR
56#endif
57
58// lapack.h probably has a typedef for lapack_int that we need to synchronize
59#ifndef MXNODEF_LAPACK_INT
60#ifndef lapack_int
61#define lapack_int MXLAPACK_INT
62#endif
63#endif
64
65// Now if not using MKL, we bring in lapack
66extern "C"
67{
68#ifndef MXLIB_MKL
69
70#include <lapacke.h>
71
72#endif
73}
74
75namespace mx
76{
77namespace math
78{
79
80/// Determine machine parameters.
81/** Wrapper for Lapack xLAMCH
82 *
83 * See more details at http://www.netlib.org/lapack/lapack-3.1.1/html/slamch.f.html.
84 *
85 * * \tparam dataT is the data type of the arrays, and determines which underlying Lapack routine is called.
86 *
87 * \param[in] CMACH Specifies the value to be returned by SLAMCH: <pre>
88 * = 'E' or 'e', returns eps, the relative machine precision
89 * = 'S' or 's , returns sfmin, the safe minimum, such that 1/sfmin does not overflow
90 * = 'B' or 'b', returns base, the base of the machine
91 * = 'P' or 'p', returns eps*base
92 * = 'N' or 'n', returns t, the number of (base) digits in the mantissa
93 * = 'R' or 'r', returns rnd, 1.0 when rounding occurs in addition, 0.0 otherwise
94 * = 'M' or 'm', returns emin, the minimum exponent before (gradual) underflow
95 * = 'U' or 'u', returns rmin, the underflow threshold - base^(emin-1)
96 * = 'L' or 'l', returns emax, the largest exponent before overflow
97 * = 'O' or 'o', returns rmax, the overflow threshold - (base^emax)*(1-eps)
98 * </pre>
99 *
100 * \returns the value of the specified machine parameters for the specified precision
101 *
102 * \ingroup template_lapack
103 */
104template <typename dataT>
105dataT lamch( char CMACH )
106{
107 return -1;
108}
109
110/*
111extern "C"
112{
113 extern float slamch_ (MKL_CONST_PTR char *CMACHp);
114 extern double dlamch_ (MKL_CONST_PTR char *CMACHp);
115}*/
116
117// Float specialization of lamch, a wrapper for Lapack SLAMCH
118template <>
119float lamch<float>( char CMACH );
120
121// Double specialization of lamch, a wrapper for Lapack DLAMCH
122template <>
123double lamch<double>( char CMACH );
124
125/// Compute the Cholesky factorization of a real symmetric positive definite matrix A.
126/**
127 * The factorization has the form
128 * A = U**T * U, if UPLO = 'U', or
129 * A = L * L**T, if UPLO = 'L',
130 * where U is an upper triangular matrix and L is lower triangular.
131 */
132template <typename dataT>
133MXLAPACK_INT potrf(
134 char UPLO, ///< [in] 'U' if upper triangle of A is stored, 'L' if lower triangle of A is stored.
135 MXLAPACK_INT N, ///< [in] The order of the matrix A, >= 0.
136 dataT *A, ///< [in.out] Symmetric matrix of dimension (LDA,N), stored as specified in UPLO. Note that the opposite
137 ///< half is not referenced.
138 MXLAPACK_INT LDA, ///< [in] The leading dimension of A.
139 MXLAPACK_INT &INFO ///< [out] 0 on success, < 0 -INFO means the i-th argument had an illegal value, >0 the leading
140 ///< minor of order INFO is not positive definite, and the factorization could not be completed.
141);
142
143template <>
144MXLAPACK_INT potrf<float>( char UPLO, MXLAPACK_INT N, float *A, MXLAPACK_INT LDA, MXLAPACK_INT &INFO );
145
146template <>
147MXLAPACK_INT potrf<double>( char UPLO, MXLAPACK_INT N, double *A, MXLAPACK_INT LDA, MXLAPACK_INT &INFO );
148
149template <>
150MXLAPACK_INT
151potrf<std::complex<float>>( char UPLO, MXLAPACK_INT N, std::complex<float> *A, MXLAPACK_INT LDA, MXLAPACK_INT &INFO );
152
153template <>
154MXLAPACK_INT
155potrf<std::complex<double>>( char UPLO, MXLAPACK_INT N, std::complex<double> *A, MXLAPACK_INT LDA, MXLAPACK_INT &INFO );
156
157/// Reduce a real symmetric matrix to real symmetric tridiagonal form by an orthogonal similarity transformation
158/** xSYTRD reduces a real symmetric matrix A to real symmetric
159 * tridiagonal form T by an orthogonal similarity transformation:
160 *
161 * \f$ Q^T * A * Q = T. \f$
162 *
163 * For more see: http://www.netlib.org/lapack/lapack-3.1.1/html/ssytrd.f.html
164 *
165 * \tparam dataT is the data type
166 *
167 * \param UPLO 'U': Upper triangle of A is stored, 'L': Lower triangle of A is stored.
168 * \param N The order of the matrix A. N >= 0.
169 * \param A array, dimension (LDA,N)
170 * \param LDA The leading dimension of the array A. LDA >= max(1,N).
171 * \param D (output) array, dimension (N), the diagonal elements of the tridiagonal matrix T:
172 * \param E (output) array, dimension (N-1), the off-diagonal elements of the tridiagonal matrix T:
173 * \param TAU (output) array, dimension (N-1), the scalar factors of the elementary reflectors (see Further Details).
174 * \param WORK (workspace/output) array, dimension (MAX(1,LWORK)) On exit, if INFO = 0, WORK(1) returns the optimal
175 * LWORK. \param LWORK (input) The dimension of the array WORK. LWORK >= 1. For optimum performance LWORK >= N*NB,
176 * where NB is the optimal blocksize. \param INFO (output) 0: successful exit < 0: if INFO = -i, the i-th argument
177 * had an illegal value
178 *
179 * \returns the value of INFO from the LAPACK routine
180 *
181 * \ingroup template_lapack
182 */
183template <typename dataT>
184MXLAPACK_INT sytrd( char UPLO,
185 MXLAPACK_INT N,
186 dataT *A,
187 MXLAPACK_INT LDA,
188 dataT *D,
189 dataT *E,
190 dataT *TAU,
191 dataT *WORK,
192 MXLAPACK_INT LWORK,
193 MXLAPACK_INT INFO )
194{
195 return -1;
196}
197
198/*
199//Declarations of the actual LAPACK functions
200extern "C"
201{
202 void ssytrd_( MKL_CONST_PTR char * UPLO, MKL_CONST_PTR MXLAPACK_INT * N, float * A, MKL_CONST_PTR MXLAPACK_INT * LDA,
203float *D, float *E, float* TAU, float *WORK, MKL_CONST_PTR MXLAPACK_INT *LWORK, MXLAPACK_INT *INFO ); void dsytrd_(
204MKL_CONST_PTR char * UPLO, MKL_CONST_PTR MXLAPACK_INT * N, double * A, MKL_CONST_PTR MXLAPACK_INT * LDA, double *D,
205double *E, double* TAU, double *WORK, MKL_CONST_PTR MXLAPACK_INT *LWORK, MXLAPACK_INT *INFO );
206}*/
207
208template <>
209MXLAPACK_INT sytrd<float>( char UPLO,
210 MXLAPACK_INT N,
211 float *A,
212 MXLAPACK_INT LDA,
213 float *D,
214 float *E,
215 float *TAU,
216 float *WORK,
217 MXLAPACK_INT LWORK,
218 MXLAPACK_INT INFO );
219
220template <>
221MXLAPACK_INT sytrd<double>( char UPLO,
222 MXLAPACK_INT N,
223 double *A,
224 MXLAPACK_INT LDA,
225 double *D,
226 double *E,
227 double *TAU,
228 double *WORK,
229 MXLAPACK_INT LWORK,
230 MXLAPACK_INT INFO );
231
232/// Compute selected eigenvalues and, optionally, eigenvectors of a real symmetric matrix
233/** xSYEVR computes selected eigenvalues and, optionally, eigenvectors
234 * of a real symmetric matrix A. Eigenvalues and eigenvectors can be
235 * selected by specifying either a range of values or a range of
236 * indices for the desired eigenvalues.
237 *
238 * See more details: http://www.netlib.org/lapack/lapack-3.1.1/html/ssyevr.f.html.
239 *
240 * \ingroup template_lapack
241 */
242template <typename dataT>
243MXLAPACK_INT syevr( char JOBZ,
244 char RANGE,
245 char UPLO,
246 MXLAPACK_INT N,
247 dataT *A,
248 MXLAPACK_INT LDA,
249 dataT VL,
250 dataT VU,
251 MXLAPACK_INT IL,
252 MXLAPACK_INT IU,
253 dataT ABSTOL,
254 MXLAPACK_INT *M,
255 dataT *W,
256 dataT *Z,
257 MXLAPACK_INT LDZ,
258 MXLAPACK_INT *ISUPPZ,
259 dataT *WORK,
260 MXLAPACK_INT LWORK,
261 MXLAPACK_INT *IWORK,
262 MXLAPACK_INT LIWORK )
263{
264 return -1;
265}
266
267/*
268extern "C"
269{
270 void ssyevr_ (MKL_CONST_PTR char *JOBZp, MKL_CONST_PTR char *RANGEp, MKL_CONST_PTR char *UPLOp, MKL_CONST_PTR
271MXLAPACK_INT *Np, float *A, MKL_CONST_PTR MXLAPACK_INT *LDAp, MKL_CONST_PTR float *VLp, MKL_CONST_PTR float *VUp,
272 MKL_CONST_PTR MXLAPACK_INT *ILp, MKL_CONST_PTR MXLAPACK_INT *IUp, MKL_CONST_PTR float *ABSTOLp,
273MXLAPACK_INT *Mp, float *W, float *Z, MKL_CONST_PTR MXLAPACK_INT *LDZp, MXLAPACK_INT *ISUPPZ, float *WORK, MKL_CONST_PTR
274MXLAPACK_INT *LWORKp, MXLAPACK_INT *IWORK, MKL_CONST_PTR MXLAPACK_INT *LIWORKp, MXLAPACK_INT *INFOp);
275
276
277 void dsyevr_ (MKL_CONST_PTR char *JOBZp, MKL_CONST_PTR char *RANGEp, MKL_CONST_PTR char *UPLOp, MKL_CONST_PTR
278MXLAPACK_INT *Np, double *A, MKL_CONST_PTR MXLAPACK_INT *LDAp, MKL_CONST_PTR double *VLp, MKL_CONST_PTR double *VUp,
279 MKL_CONST_PTR MXLAPACK_INT *ILp, MKL_CONST_PTR MXLAPACK_INT *IUp, MKL_CONST_PTR double *ABSTOLp,
280MXLAPACK_INT *Mp, double *W, double *Z, MKL_CONST_PTR MXLAPACK_INT *LDZp, MXLAPACK_INT *ISUPPZ, double *WORK,
281MKL_CONST_PTR MXLAPACK_INT *LWORKp, MXLAPACK_INT *IWORK, MKL_CONST_PTR MXLAPACK_INT *LIWORKp, MXLAPACK_INT *INFOp);
282}
283*/
284
285// Float specialization of syevr, a wrapper for Lapack SSYEVR
286template <>
287MXLAPACK_INT syevr<float>( char JOBZ,
288 char RANGE,
289 char UPLO,
290 MXLAPACK_INT N,
291 float *A,
292 MXLAPACK_INT LDA,
293 float VL,
294 float VU,
295 MXLAPACK_INT IL,
296 MXLAPACK_INT IU,
297 float ABSTOL,
298 MXLAPACK_INT *M,
299 float *W,
300 float *Z,
301 MXLAPACK_INT LDZ,
302 MXLAPACK_INT *ISUPPZ,
303 float *WORK,
304 MXLAPACK_INT LWORK,
305 MXLAPACK_INT *IWORK,
306 MXLAPACK_INT LIWORK );
307
308// Double specialization of syevr, a wrapper for Lapack DSYEVR
309template <>
310MXLAPACK_INT syevr<double>( char JOBZ,
311 char RANGE,
312 char UPLO,
313 MXLAPACK_INT N,
314 double *A,
315 MXLAPACK_INT LDA,
316 double VL,
317 double VU,
318 MXLAPACK_INT IL,
319 MXLAPACK_INT IU,
320 double ABSTOL,
321 MXLAPACK_INT *M,
322 double *W,
323 double *Z,
324 MXLAPACK_INT LDZ,
325 MXLAPACK_INT *ISUPPZ,
326 double *WORK,
327 MXLAPACK_INT LWORK,
328 MXLAPACK_INT *IWORK,
329 MXLAPACK_INT LIWORK );
330
331/// Compute the singular value decomposition (SVD) of a real matrix
332/** xGESVD computes the singular value decomposition (SVD) of a real
333 * M-by-N matrix A, optionally computing the left and/or right singular
334 * vectors. The SVD is written
335 *
336 * \f$ A = U * \Sigma * V^T \f$
337 *
338 * where \f$ \Sigma \f$ is an M-by-N matrix which is zero except for its
339 * min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and
340 * V is an N-by-N orthogonal matrix. The diagonal elements of \f$ \Sigma \f$
341 * are the singular values of A; they are real and non-negative, and
342 * are returned in descending order. The first min(m,n) columns of
343 * U and V are the left and right singular vectors of A.
344 *
345 * Note that the routine returns \f$ V^T \f$, not \f$ V \f$.
346 *
347 * See more details: http://www.netlib.org/lapack/explore-html/d8/d49/sgesvd_8f.html.
348 * This documentation is taken from there.
349 *
350 * \tparam dataT is the data type of the arrays, and determines which underlying Lapack routine is called.
351 *
352 * \param[in] JOBU
353 * \parblock
354 (char) Specifies options for computing all or part of the matrix U: <br />
355 = 'A': all M columns of U are returned in array U <br />
356 = 'S': the first min(m,n) columns of U (the left singular vectors)
357 are returned in the array U <br />
358 = 'O': the first min(m,n) columns of U (the left singular <br />
359 vectors) are overwritten on the array A <br />
360 = 'N': no columns of U (no left singular vectors) are
361 computed.
362
363 * \param[in] JOBVT
364 * \parblock
365 (char) Specifies options for computing all or part of the matrix V**T:
366 = 'A': all N rows of V**T are returned in the array VT <br />
367 = 'S': the first min(m,n) rows of V**T (the right singular
368 vectors) are returned in the array VT <br />
369 = 'O': the first min(m,n) rows of V**T (the right singular
370 vectors) are overwritten on the array A<br />
371 = 'N': no rows of V**T (no right singular vectors) are
372 computed.
373
374 JOBVT and JOBU cannot both be 'O'.
375 *
376 * \param[in] M
377 * \parblock
378 (MXLAPACK_INT)
379 The number of rows of the input matrix A. M >= 0.
380 *
381 * \param[in] N
382 * \parblock
383 (MXLAPACK_INT)
384 The number of columns of the input matrix A. N >= 0.
385 *
386 * \param[in,out] A
387 * \parblock
388 (dataT *, dimension (LDA,N))
389 On entry: the M-by-N matrix A.<br />
390 On exit:<br />
391 if JOBU = 'O', A is overwritten with the first min(m,n)
392 columns of U (the left singular vectors,
393 stored columnwise)<br />
394 if JOBVT = 'O', A is overwritten with the first min(m,n)
395 rows of V**T (the right singular vectors,
396 stored rowwise)<br />
397 if JOBU != 'O' and JOBVT .ne. 'O', the contents of A
398 are destroyed.
399 *
400 * \param[in] LDA
401 * \parblock
402 (MXLAPACK_INT)
403 The leading dimension of the array A. LDA >= max(1,M).
404 *
405 * \param[out] S
406 * \parblock
407 (dataT *, dimension (min(M,N))
408 The singular values of A, sorted so that S(i) >= S(i+1).
409 *
410 * \param[out] U
411 * \parblock
412 (dataT *, dimension (LDU,UCOL))
413 (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'.<br />
414 If JOBU = 'A', U contains the M-by-M orthogonal matrix U<br />
415 if JOBU = 'S', U contains the first min(m,n) columns of U
416 (the left singular vectors, stored columnwise)<br />
417 if JOBU = 'N' or 'O', U is not referenced.
418 *
419 * \param[in] LDU
420 * \parblock
421 (MXLAPACK_INT)
422 The leading dimension of the array U. <br />
423 LDU >= 1<br />
424 if JOBU == 'S' or 'A', LDU >= M.
425 *
426 * \param[out] VT
427 * \parblock
428 (dataT *, dimension (LDVT,N))
429 If JOBVT = 'A', VT contains the N-by-N orthogonal matrix
430 V**T <br />
431 if JOBVT = 'S', VT contains the first min(m,n) rows of
432 V**T (the right singular vectors, stored rowwise) <br />
433 if JOBVT = 'N' or 'O', VT is not referenced.
434 *
435 * \param[in] LDVT
436 * \parblock
437 (MXLAPACK_INT)
438 The leading dimension of the array VT. <br />
439 LDVT >= 1<br />
440 if JOBVT = 'A', LDVT >= N<br />
441 if JOBVT = 'S', LDVT >= min(M,N).
442 *
443 * \param[out] WORK
444 \parblock
445 (dataT *, dimension (MAX(1,LWORK)) )<br />
446 On exit, if INFO = 0, WORK[1] returns the optimal LWORK<br />
447 if INFO > 0, WORK(2:MIN(M,N)) contains the unconverged
448 superdiagonal elements of an upper bidiagonal matrix B
449 whose diagonal is in S (not necessarily sorted). B
450 satisfies A = U * B * VT, so it has the same singular values
451 as A, and singular vectors related by U and VT.
452 *
453 * \param[in] LWORK
454 \parblock
455 (MXLAPACK_INT)
456 The dimension of the array WORK.<br />
457 LWORK >= MAX(1,5*MIN(M,N)) for the paths (see comments inside code):
458 - PATH 1 (M much larger than N, JOBU='N')
459 - PATH 1t (N much larger than M, JOBVT='N')
460 LWORK >= MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(M,N)) for the other paths
461 For good performance, LWORK should generally be larger.
462
463 If LWORK = -1, then a workspace query is assumed; the routine
464 only calculates the optimal size of the WORK array, returns
465 this value as the first entry of the WORK array, and no error
466 message related to LWORK is issued by XERBLA.
467 * \endparblock
468 * \returns
469 * \parblock
470 =0: successful exit. <br />
471 <0: if INFO = -i, the i-th argument had an illegal value. <br />
472 >0: if SBDSQR did not converge, INFO specifies how many
473 superdiagonals of an MXLAPACK_INTermediate bidiagonal form B
474 did not converge to zero. See the description of WORK
475 above for details.
476 \endparblock
477 *
478 * \ingroup template_lapack
479 */
480template <typename dataT>
481MXLAPACK_INT gesvd( char JOBU,
482 char JOBVT,
483 MXLAPACK_INT M,
484 MXLAPACK_INT N,
485 dataT *A,
486 MXLAPACK_INT LDA,
487 dataT *S,
488 dataT *U,
489 MXLAPACK_INT LDU,
490 dataT *VT,
491 MXLAPACK_INT LDVT,
492 dataT *WORK,
493 MXLAPACK_INT LWORK )
494{
495 return -1;
496}
497
498/*
499//Declarations of the lapack calls
500extern "C"
501{
502 void sgesvd_( MKL_CONST_PTR char *JOBUp, MKL_CONST_PTR char *JOBVTp, MKL_CONST_PTR MXLAPACK_INT *Mp, MKL_CONST_PTR
503MXLAPACK_INT *Np, float * A, MKL_CONST_PTR MXLAPACK_INT *LDAp, float * S, float *U, MKL_CONST_PTR MXLAPACK_INT *LDUp,
504 float * VT, MKL_CONST_PTR MXLAPACK_INT *LDVTp, float * WORK, MKL_CONST_PTR MXLAPACK_INT *LWORKp,
505MXLAPACK_INT *INFOp);
506
507 void dgesvd_( MKL_CONST_PTR char *JOBUp, MKL_CONST_PTR char *JOBVTp, MKL_CONST_PTR MXLAPACK_INT *Mp, MKL_CONST_PTR
508MXLAPACK_INT *Np, double * A, MKL_CONST_PTR MXLAPACK_INT *LDAp, double * S, double *U, MKL_CONST_PTR MXLAPACK_INT *LDUp,
509 double * VT, MKL_CONST_PTR MXLAPACK_INT *LDVTp, double * WORK, MKL_CONST_PTR MXLAPACK_INT *LWORKp,
510MXLAPACK_INT *INFOp);
511}*/
512
513// float specialization of gesvd
514template <>
515MXLAPACK_INT gesvd<float>( char JOBU,
516 char JOBVT,
517 MXLAPACK_INT M,
518 MXLAPACK_INT N,
519 float *A,
520 MXLAPACK_INT LDA,
521 float *S,
522 float *U,
523 MXLAPACK_INT LDU,
524 float *VT,
525 MXLAPACK_INT LDVT,
526 float *WORK,
527 MXLAPACK_INT LWORK );
528
529// double specialization of gesvd
530template <>
531MXLAPACK_INT gesvd<double>( char JOBU,
532 char JOBVT,
533 MXLAPACK_INT M,
534 MXLAPACK_INT N,
535 double *A,
536 MXLAPACK_INT LDA,
537 double *S,
538 double *U,
539 MXLAPACK_INT LDU,
540 double *VT,
541 MXLAPACK_INT LDVT,
542 double *WORK,
543 MXLAPACK_INT LWORK );
544
545/// Compute the singular value decomposition (SVD) of a real matrix with GESDD
546/**
547 This documentation copied shamelessly from the LAPACK source at <a
548href="http://www.netlib.org/lapack/explore-html/d4/dca/group__real_g_esing.html#gac2cd4f1079370ac908186d77efcd5ea8">netlib</a>.
549
550 SGESDD computes the singular value decomposition (SVD) of a real
551 M-by-N matrix A, optionally computing the left and right singular
552 vectors. If singular vectors are desired, it uses a
553 divide-and-conquer algorithm.
554
555 The SVD is written
556
557\f[ A = U \Sigma V^T \f]
558
559 where \f$ \Sigma \f$ is an M-by-N matrix which is zero except for its
560 min(m,n) diagonal elements, \f$ U \f$ is an M-by-M orthogonal matrix, and
561 \f$ V is an N-by-N orthogonal matrix. The diagonal elements of \f$ \Sigma \f$
562 are the singular values of A; they are real and non-negative, and
563 are returned in descending order. The first min(m,n) columns of
564 U and V are the left and right singular vectors of A.
565
566 Note that the routine returns \f$ V^T \f$, not \f$ V \f$.
567
568 The divide and conquer algorithm makes very mild assumptions about
569 floating point arithmetic. It will work on machines with a guard
570 digit in add/subtract, or on those binary machines without guard
571 digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
572 Cray-2. It could conceivably fail on hexadecimal or decimal machines
573 without guard digits, but the authors know of none.
574 \endparblock
575
576 \param[in] JOBZ
577 \parblock
578 JOBZ is CHARACTER*1
579 Specifies options for computing all or part of the matrix U:
580 = 'A': all M columns of U and all N rows of V**T are
581 returned in the arrays U and VT;
582 = 'S': the first min(M,N) columns of U and the first
583 min(M,N) rows of V**T are returned in the arrays U
584 and VT;
585 = 'O': If M >= N, the first N columns of U are overwritten
586 on the array A and all rows of V**T are returned in
587 the array VT;
588 otherwise, all columns of U are returned in the
589 array U and the first M rows of V**T are overwritten
590 in the array A;
591 = 'N': no columns of U or rows of V**T are computed.
592 \param[in] M
593 \parblock
594 M is INTEGER
595 The number of rows of the input matrix A. M >= 0.
596 \param[in] N
597 \parblock
598 N is INTEGER
599 The number of columns of the input matrix A. N >= 0.
600 \param[in,out] A
601 \parblock
602 A is REAL array, dimension (LDA,N)
603 On entry, the M-by-N matrix A.
604 On exit,
605 if JOBZ = 'O', A is overwritten with the first N columns
606 of U (the left singular vectors, stored
607 columnwise) if M >= N;
608 A is overwritten with the first M rows
609 of V**T (the right singular vectors, stored
610 rowwise) otherwise.
611 if JOBZ .ne. 'O', the contents of A are destroyed.
612 \param[in] LDA
613 \parblock
614 LDA is INTEGER
615 The leading dimension of the array A. LDA >= max(1,M).
616
617
618 \param[out] S
619 \parblock
620 S is REAL array, dimension (min(M,N))
621 The singular values of A, sorted so that S(i) >= S(i+1).
622
623 \param[out] U
624 \parblock
625 U is REAL array, dimension (LDU,UCOL)
626 UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N;
627 UCOL = min(M,N) if JOBZ = 'S'.
628 If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M
629 orthogonal matrix U;
630 if JOBZ = 'S', U contains the first min(M,N) columns of U
631 (the left singular vectors, stored columnwise);
632 if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced.
633
634 \param[in] LDU
635 \parblock
636 LDU is INTEGER
637 The leading dimension of the array U. LDU >= 1; if
638 JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M.
639
640 \param[out] VT
641 \parblock
642 VT is REAL array, dimension (LDVT,N)
643 If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the
644 N-by-N orthogonal matrix V**T;
645 if JOBZ = 'S', VT contains the first min(M,N) rows of
646 V**T (the right singular vectors, stored rowwise);
647 if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced.
648
649 \param[in] LDVT
650 \parblock
651 LDVT is INTEGER
652 The leading dimension of the array VT. LDVT >= 1; if
653 JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N;
654 if JOBZ = 'S', LDVT >= min(M,N).
655
656 \param[out] WORK
657 \parblock
658 WORK is REAL array, dimension (MAX(1,LWORK))
659 On exit, if INFO = 0, WORK(1) returns the optimal LWORK;
660
661 \param[in] LWORK
662 \parblock
663 LWORK is INTEGER
664 The dimension of the array WORK. LWORK >= 1.
665 If JOBZ = 'N',
666 LWORK >= 3*min(M,N) + max(max(M,N),6*min(M,N)).
667 If JOBZ = 'O',
668 LWORK >= 3*min(M,N) +
669 max(max(M,N),5*min(M,N)*min(M,N)+4*min(M,N)).
670 If JOBZ = 'S' or 'A'
671 LWORK >= min(M,N)*(7+4*min(M,N))
672 For good performance, LWORK should generally be larger.
673 If LWORK = -1 but other input arguments are legal, WORK(1)
674 returns the optimal LWORK.
675
676 \param[out] IWORK
677 \parblock
678 IWORK is INTEGER array, dimension (8*min(M,N))
679 \endparblock
680
681 \returns
682 \parblock
683 = 0: successful exit.<br />
684 < 0: if INFO = -i, the i-th argument had an illegal value.<br />
685 > 0: SBDSDC did not converge, updating process failed.
686 \endparblock
687
688\ingroup template_lapack
689
690*/
691template <typename dataT>
692MXLAPACK_INT gesdd( char JOBZ,
693 MXLAPACK_INT M,
694 MXLAPACK_INT N,
695 dataT *A,
696 MXLAPACK_INT LDA,
697 dataT *S,
698 dataT *U,
699 MXLAPACK_INT LDU,
700 dataT *VT,
701 MXLAPACK_INT LDVT,
702 dataT *WORK,
703 MXLAPACK_INT LWORK,
704 MXLAPACK_INT *IWORK,
705 MXLAPACK_INT INFO )
706{
707 return -1;
708}
709
710/*
711//Declarations of the lapack calls
712extern "C"
713{
714 void sgesdd_(MKL_CONST_PTR char * JOBZ, MKL_CONST_PTR MXLAPACK_INT *M, MKL_CONST_PTR MXLAPACK_INT *N, float *A,
715MKL_CONST_PTR MXLAPACK_INT *LDA, float *S, float * U, MKL_CONST_PTR MXLAPACK_INT *LDU, float * VT, MKL_CONST_PTR
716MXLAPACK_INT *LDVT, float *WORK, MKL_CONST_PTR MXLAPACK_INT * LWORK, MXLAPACK_INT * IWORK, MXLAPACK_INT *INFO);
717
718 void dgesdd_(MKL_CONST_PTR char * JOBZ, MKL_CONST_PTR MXLAPACK_INT *M, MKL_CONST_PTR MXLAPACK_INT *N, double *A,
719MKL_CONST_PTR MXLAPACK_INT *LDA, double *S, double * U, MKL_CONST_PTR MXLAPACK_INT *LDU, double * VT, MKL_CONST_PTR
720MXLAPACK_INT *LDVT, double *WORK, MKL_CONST_PTR MXLAPACK_INT * LWORK, MXLAPACK_INT * IWORK, MXLAPACK_INT *INFO);
721
722}*/
723
724// float specialization of gesdd
725template <>
726MXLAPACK_INT gesdd<float>( char JOBZ,
727 MXLAPACK_INT M,
728 MXLAPACK_INT N,
729 float *A,
730 MXLAPACK_INT LDA,
731 float *S,
732 float *U,
733 MXLAPACK_INT LDU,
734 float *VT,
735 MXLAPACK_INT LDVT,
736 float *WORK,
737 MXLAPACK_INT LWORK,
738 MXLAPACK_INT *IWORK,
739 MXLAPACK_INT INFO );
740
741// double specialization of gesdd
742template <>
743MXLAPACK_INT gesdd<double>( char JOBZ,
744 MXLAPACK_INT M,
745 MXLAPACK_INT N,
746 double *A,
747 MXLAPACK_INT LDA,
748 double *S,
749 double *U,
750 MXLAPACK_INT LDU,
751 double *VT,
752 MXLAPACK_INT LDVT,
753 double *WORK,
754 MXLAPACK_INT LWORK,
755 MXLAPACK_INT *IWORK,
756 MXLAPACK_INT INFO );
757
758} // namespace math
759} // namespace mx
760
761//
762
763#endif // math_templateLapack_hpp
mx::math::six_fifths
constexpr floatT six_fifths()
Return 6/5 in the specified precision.
Definition constants.hpp:404
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:184
mx::math::syevr
MXLAPACK_INT syevr(char JOBZ, char RANGE, char UPLO, MXLAPACK_INT N, dataT *A, MXLAPACK_INT LDA, dataT VL, dataT VU, MXLAPACK_INT IL, MXLAPACK_INT IU, dataT ABSTOL, MXLAPACK_INT *M, dataT *W, dataT *Z, MXLAPACK_INT LDZ, MXLAPACK_INT *ISUPPZ, dataT *WORK, MXLAPACK_INT LWORK, MXLAPACK_INT *IWORK, MXLAPACK_INT LIWORK)
Compute selected eigenvalues and, optionally, eigenvectors of a real symmetric matrix.
Definition templateLapack.hpp:243
mx::math::lamch
dataT lamch(char CMACH)
Determine machine parameters.
Definition templateLapack.hpp:105
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:481
mx
The mxlib c++ namespace.
Definition mxError.hpp:106
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.