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