diff options
Diffstat (limited to 'ml/dlib/dlib/matrix/lapack')
| -rw-r--r-- | ml/dlib/dlib/matrix/lapack/fortran_id.h | 62 | ||||
| -rw-r--r-- | ml/dlib/dlib/matrix/lapack/gees.h | 264 | ||||
| -rw-r--r-- | ml/dlib/dlib/matrix/lapack/geev.h | 234 | ||||
| -rw-r--r-- | ml/dlib/dlib/matrix/lapack/geqrf.h | 168 | ||||
| -rw-r--r-- | ml/dlib/dlib/matrix/lapack/gesdd.h | 364 | ||||
| -rw-r--r-- | ml/dlib/dlib/matrix/lapack/gesvd.h | 323 | ||||
| -rw-r--r-- | ml/dlib/dlib/matrix/lapack/getrf.h | 132 | ||||
| -rw-r--r-- | ml/dlib/dlib/matrix/lapack/ormqr.h | 224 | ||||
| -rw-r--r-- | ml/dlib/dlib/matrix/lapack/pbtrf.h | 178 | ||||
| -rw-r--r-- | ml/dlib/dlib/matrix/lapack/potrf.h | 174 | ||||
| -rw-r--r-- | ml/dlib/dlib/matrix/lapack/syev.h | 218 | ||||
| -rw-r--r-- | ml/dlib/dlib/matrix/lapack/syevr.h | 445 |
12 files changed, 2786 insertions, 0 deletions
diff --git a/ml/dlib/dlib/matrix/lapack/fortran_id.h b/ml/dlib/dlib/matrix/lapack/fortran_id.h new file mode 100644 index 000000000..8027ea34f --- /dev/null +++ b/ml/dlib/dlib/matrix/lapack/fortran_id.h @@ -0,0 +1,62 @@ +// Copyright (C) 2010 Davis E. King (davis@dlib.net) +// License: Boost Software License See LICENSE.txt for the full license. +#ifndef DLIB_BOOST_NUMERIC_BINDINGS_TRAITS_FORTRAN_H +#define DLIB_BOOST_NUMERIC_BINDINGS_TRAITS_FORTRAN_H + +// ---------------------------------------------------------------------------------------- +// ---------------------------------------------------------------------------------------- +// FORTRAN BINDING STUFF FROM BOOST +// ---------------------------------------------------------------------------------------- +// ---------------------------------------------------------------------------------------- + +// Permission to copy, use, modify, sell and +// distribute this software is granted provided this copyright notice appears +// in all copies. This software is provided "as is" without express or implied +// warranty, and with no claim as to its suitability for any purpose. +// Copyright (C) 2002, 2003 Si-Lab b.v.b.a., Toon Knapen and Kresimir Fresl + + +// First we need to know what the conventions for linking +// C with Fortran is on this platform/toolset +#if defined(LAPACK_FORCE_UNDERSCORE) +#define DLIB_BIND_FORTRAN_LOWERCASE_UNDERSCORE +#elif defined(LAPACK_FORCE_NOUNDERSCORE) +#define DLIB_BIND_FORTRAN_LOWERCASE +#elif defined(__GNUC__) || defined(__ICC) || defined(__sgi) || defined(__COMO__) || defined(__KCC) +#define DLIB_BIND_FORTRAN_LOWERCASE_UNDERSCORE +#elif defined(__IBMCPP__) || defined(_MSC_VER) || defined(__BORLANDC__) +#define DLIB_BIND_FORTRAN_LOWERCASE +#else +#error do not know how to link with fortran for the given platform +#endif + +// Next we define macros to convert our symbols to +// the current convention +#if defined(DLIB_BIND_FORTRAN_LOWERCASE_UNDERSCORE) +#define DLIB_FORTRAN_ID( id ) id##_ +#elif defined(DLIB_BIND_FORTRAN_LOWERCASE) +#define DLIB_FORTRAN_ID( id ) id +#else +#error do not know how to bind to fortran calling convention +#endif + + + +namespace dlib +{ + namespace lapack + { + // stuff from f2c used to define what exactly is an integer in fortran +#if (defined(__alpha__) || defined(__sparc64__) || defined(__x86_64__) || defined(__ia64__)) && !defined(MATLAB_MEX_FILE) + typedef int integer; + typedef unsigned int uinteger; +#else + typedef long int integer; + typedef unsigned long int uinteger; +#endif + + } +} + +#endif // DLIB_BOOST_NUMERIC_BINDINGS_TRAITS_FORTRAN_H + diff --git a/ml/dlib/dlib/matrix/lapack/gees.h b/ml/dlib/dlib/matrix/lapack/gees.h new file mode 100644 index 000000000..a8ee63ff1 --- /dev/null +++ b/ml/dlib/dlib/matrix/lapack/gees.h @@ -0,0 +1,264 @@ +// Copyright (C) 2010 Davis E. King (davis@dlib.net) +// License: Boost Software License See LICENSE.txt for the full license. +#ifndef DLIB_LAPACk_ES_Hh_ +#define DLIB_LAPACk_ES_Hh_ + +#include "fortran_id.h" +#include "../matrix.h" + +namespace dlib +{ + namespace lapack + { + namespace binding + { +#if defined(__alpha__) || defined(__sparc64__) || defined(__x86_64__) || defined(__ia64__) + typedef int logical; +#else + typedef long int logical; +#endif + typedef logical (*L_fp)(...); + + extern "C" + { + void DLIB_FORTRAN_ID(dgees) (char *jobvs, char *sort, L_fp select, integer *n, + double *a, integer *lda, integer *sdim, double *wr, + double *wi, double *vs, integer *ldvs, double *work, + integer *lwork, logical *bwork, integer *info); + + void DLIB_FORTRAN_ID(sgees) (char *jobvs, char *sort, L_fp select, integer *n, + float *a, integer *lda, integer *sdim, float *wr, + float *wi, float *vs, integer *ldvs, float *work, + integer *lwork, logical *bwork, integer *info); + + } + + inline int gees (char jobvs, integer n, + double *a, integer lda, double *wr, + double *wi, double *vs, integer ldvs, double *work, + integer lwork) + { + // No sorting allowed + integer info = 0; + char sort = 'N'; + L_fp fnil = 0; + logical bwork = 0; + integer sdim = 0; + DLIB_FORTRAN_ID(dgees)(&jobvs, &sort, fnil, &n, + a, &lda, &sdim, wr, + wi, vs, &ldvs, work, + &lwork, &bwork, &info); + return info; + } + + + inline int gees (char jobvs, integer n, + float *a, integer lda, float *wr, + float *wi, float *vs, integer ldvs, float *work, + integer lwork) + { + // No sorting allowed + integer info = 0; + char sort = 'N'; + L_fp fnil = 0; + logical bwork = 0; + integer sdim = 0; + DLIB_FORTRAN_ID(sgees)(&jobvs, &sort, fnil, &n, + a, &lda, &sdim, wr, + wi, vs, &ldvs, work, + &lwork, &bwork, &info); + return info; + } + + + } + + // ------------------------------------------------------------------------------------ + +/* -- LAPACK driver routine (version 3.1) -- */ +/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ +/* November 2006 */ + +/* .. Scalar Arguments .. */ +/* .. */ +/* .. Array Arguments .. */ +/* .. */ +/* .. Function Arguments .. */ +/* .. */ + +/* Purpose */ +/* ======= */ + +/* DGEES computes for an N-by-N real nonsymmetric matrix A, the */ +/* eigenvalues, the real Schur form T, and, optionally, the matrix of */ +/* Schur vectors Z. This gives the Schur factorization A = Z*T*(Z**T). */ + +/* Optionally, it also orders the eigenvalues on the diagonal of the */ +/* real Schur form so that selected eigenvalues are at the top left. */ +/* The leading columns of Z then form an orthonormal basis for the */ +/* invariant subspace corresponding to the selected eigenvalues. */ + +/* A matrix is in real Schur form if it is upper quasi-triangular with */ +/* 1-by-1 and 2-by-2 blocks. 2-by-2 blocks will be standardized in the */ +/* form */ +/* [ a b ] */ +/* [ c a ] */ + +/* where b*c < 0. The eigenvalues of such a block are a +- sqrt(bc). */ + +/* Arguments */ +/* ========= */ + +/* JOBVS (input) CHARACTER*1 */ +/* = 'N': Schur vectors are not computed; */ +/* = 'V': Schur vectors are computed. */ + +/* SORT (input) CHARACTER*1 */ +/* Specifies whether or not to order the eigenvalues on the */ +/* diagonal of the Schur form. */ +/* = 'N': Eigenvalues are not ordered; */ +/* = 'S': Eigenvalues are ordered (see SELECT). */ + +/* SELECT (external procedure) LOGICAL FUNCTION of two DOUBLE PRECISION arguments */ +/* SELECT must be declared EXTERNAL in the calling subroutine. */ +/* If SORT = 'S', SELECT is used to select eigenvalues to sort */ +/* to the top left of the Schur form. */ +/* If SORT = 'N', SELECT is not referenced. */ +/* An eigenvalue WR(j)+sqrt(-1)*WI(j) is selected if */ +/* SELECT(WR(j),WI(j)) is true; i.e., if either one of a complex */ +/* conjugate pair of eigenvalues is selected, then both complex */ +/* eigenvalues are selected. */ +/* Note that a selected complex eigenvalue may no longer */ +/* satisfy SELECT(WR(j),WI(j)) = .TRUE. after ordering, since */ +/* ordering may change the value of complex eigenvalues */ +/* (especially if the eigenvalue is ill-conditioned); in this */ +/* case INFO is set to N+2 (see INFO below). */ + +/* N (input) INTEGER */ +/* The order of the matrix A. N >= 0. */ + +/* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) */ +/* On entry, the N-by-N matrix A. */ +/* On exit, A has been overwritten by its real Schur form T. */ + +/* LDA (input) INTEGER */ +/* The leading dimension of the array A. LDA >= max(1,N). */ + +/* SDIM (output) INTEGER */ +/* If SORT = 'N', SDIM = 0. */ +/* If SORT = 'S', SDIM = number of eigenvalues (after sorting) */ +/* for which SELECT is true. (Complex conjugate */ +/* pairs for which SELECT is true for either */ +/* eigenvalue count as 2.) */ + +/* WR (output) DOUBLE PRECISION array, dimension (N) */ +/* WI (output) DOUBLE PRECISION array, dimension (N) */ +/* WR and WI contain the real and imaginary parts, */ +/* respectively, of the computed eigenvalues in the same order */ +/* that they appear on the diagonal of the output Schur form T. */ +/* Complex conjugate pairs of eigenvalues will appear */ +/* consecutively with the eigenvalue having the positive */ +/* imaginary part first. */ + +/* VS (output) DOUBLE PRECISION array, dimension (LDVS,N) */ +/* If JOBVS = 'V', VS contains the orthogonal matrix Z of Schur */ +/* vectors. */ +/* If JOBVS = 'N', VS is not referenced. */ + +/* LDVS (input) INTEGER */ +/* The leading dimension of the array VS. LDVS >= 1; if */ +/* JOBVS = 'V', LDVS >= N. */ + +/* WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) */ +/* On exit, if INFO = 0, WORK(1) contains the optimal LWORK. */ + +/* LWORK (input) INTEGER */ +/* The dimension of the array WORK. LWORK >= max(1,3*N). */ +/* For good performance, LWORK must generally be larger. */ + +/* If LWORK = -1, then a workspace query is assumed; the routine */ +/* only calculates the optimal size of the WORK array, returns */ +/* this value as the first entry of the WORK array, and no error */ +/* message related to LWORK is issued by XERBLA. */ + +/* BWORK (workspace) LOGICAL array, dimension (N) */ +/* Not referenced if SORT = 'N'. */ + +/* INFO (output) INTEGER */ +/* = 0: successful exit */ +/* < 0: if INFO = -i, the i-th argument had an illegal value. */ +/* > 0: if INFO = i, and i is */ +/* <= N: the QR algorithm failed to compute all the */ +/* eigenvalues; elements 1:ILO-1 and i+1:N of WR and WI */ +/* contain those eigenvalues which have converged; if */ +/* JOBVS = 'V', VS contains the matrix which reduces A */ +/* to its partially converged Schur form. */ +/* = N+1: the eigenvalues could not be reordered because some */ +/* eigenvalues were too close to separate (the problem */ +/* is very ill-conditioned); */ +/* = N+2: after reordering, roundoff changed values of some */ +/* complex eigenvalues so that leading eigenvalues in */ +/* the Schur form no longer satisfy SELECT=.TRUE. This */ +/* could also be caused by underflow due to scaling. */ + + // ------------------------------------------------------------------------------------ + + template < + typename T, + long NR1, long NR2, long NR3, long NR4, + long NC1, long NC2, long NC3, long NC4, + typename MM, + typename layout + > + int gees ( + const char jobz, + matrix<T,NR1,NC1,MM,column_major_layout>& a, + matrix<T,NR2,NC2,MM,layout>& wr, + matrix<T,NR3,NC3,MM,layout>& wi, + matrix<T,NR4,NC4,MM,column_major_layout>& vs + ) + { + matrix<T,0,1,MM,column_major_layout> work; + + const long n = a.nr(); + + wr.set_size(n,1); + wi.set_size(n,1); + + if (jobz == 'V') + vs.set_size(n,n); + else + vs.set_size(NR4?NR4:1, NC4?NC4:1); + + // figure out how big the workspace needs to be. + T work_size = 1; + int info = binding::gees(jobz, n, + &a(0,0), a.nr(), &wr(0,0), + &wi(0,0), &vs(0,0), vs.nr(), &work_size, + -1); + + if (info != 0) + return info; + + if (work.size() < work_size) + work.set_size(static_cast<long>(work_size), 1); + + // compute the actual decomposition + info = binding::gees(jobz, n, + &a(0,0), a.nr(), &wr(0,0), + &wi(0,0), &vs(0,0), vs.nr(), &work(0,0), + work.size()); + + return info; + } + + // ------------------------------------------------------------------------------------ + + } + +} + +// ---------------------------------------------------------------------------------------- + +#endif // DLIB_LAPACk_ES_Hh_ + diff --git a/ml/dlib/dlib/matrix/lapack/geev.h b/ml/dlib/dlib/matrix/lapack/geev.h new file mode 100644 index 000000000..d8fdc4af5 --- /dev/null +++ b/ml/dlib/dlib/matrix/lapack/geev.h @@ -0,0 +1,234 @@ +// Copyright (C) 2010 Davis E. King (davis@dlib.net) +// License: Boost Software License See LICENSE.txt for the full license. +#ifndef DLIB_LAPACk_GEEV_Hh_ +#define DLIB_LAPACk_GEEV_Hh_ + +#include "fortran_id.h" +#include "../matrix.h" + +namespace dlib +{ + namespace lapack + { + namespace binding + { + extern "C" + { + void DLIB_FORTRAN_ID(dgeev) (char *jobvl, char *jobvr, integer *n, double * a, + integer *lda, double *wr, double *wi, double *vl, + integer *ldvl, double *vr, integer *ldvr, double *work, + integer *lwork, integer *info); + + void DLIB_FORTRAN_ID(sgeev) (char *jobvl, char *jobvr, integer *n, float * a, + integer *lda, float *wr, float *wi, float *vl, + integer *ldvl, float *vr, integer *ldvr, float *work, + integer *lwork, integer *info); + + } + + inline int geev (char jobvl, char jobvr, integer n, double *a, + integer lda, double *wr, double *wi, double *vl, + integer ldvl, double *vr, integer ldvr, double *work, + integer lwork) + { + integer info = 0; + DLIB_FORTRAN_ID(dgeev)(&jobvl, &jobvr, &n, a, + &lda, wr, wi, vl, + &ldvl, vr, &ldvr, work, + &lwork, &info); + return info; + } + + inline int geev (char jobvl, char jobvr, integer n, float *a, + integer lda, float *wr, float *wi, float *vl, + integer ldvl, float *vr, integer ldvr, float *work, + integer lwork) + { + integer info = 0; + DLIB_FORTRAN_ID(sgeev)(&jobvl, &jobvr, &n, a, + &lda, wr, wi, vl, + &ldvl, vr, &ldvr, work, + &lwork, &info); + return info; + } + + + } + + // ------------------------------------------------------------------------------------ + +/* -- LAPACK driver routine (version 3.1) -- */ +/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ +/* November 2006 */ + +/* .. Scalar Arguments .. */ +/* .. */ +/* .. Array Arguments .. */ +/* .. */ + +/* Purpose */ +/* ======= */ + +/* DGEEV computes for an N-by-N real nonsymmetric matrix A, the */ +/* eigenvalues and, optionally, the left and/or right eigenvectors. */ + +/* The right eigenvector v(j) of A satisfies */ +/* A * v(j) = lambda(j) * v(j) */ +/* where lambda(j) is its eigenvalue. */ +/* The left eigenvector u(j) of A satisfies */ +/* u(j)**H * A = lambda(j) * u(j)**H */ +/* where u(j)**H denotes the conjugate transpose of u(j). */ + +/* The computed eigenvectors are normalized to have Euclidean norm */ +/* equal to 1 and largest component real. */ + +/* Arguments */ +/* ========= */ + +/* JOBVL (input) CHARACTER*1 */ +/* = 'N': left eigenvectors of A are not computed; */ +/* = 'V': left eigenvectors of A are computed. */ + +/* JOBVR (input) CHARACTER*1 */ +/* = 'N': right eigenvectors of A are not computed; */ +/* = 'V': right eigenvectors of A are computed. */ + +/* N (input) INTEGER */ +/* The order of the matrix A. N >= 0. */ + +/* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) */ +/* On entry, the N-by-N matrix A. */ +/* On exit, A has been overwritten. */ + +/* LDA (input) INTEGER */ +/* The leading dimension of the array A. LDA >= max(1,N). */ + +/* WR (output) DOUBLE PRECISION array, dimension (N) */ +/* WI (output) DOUBLE PRECISION array, dimension (N) */ +/* WR and WI contain the real and imaginary parts, */ +/* respectively, of the computed eigenvalues. Complex */ +/* conjugate pairs of eigenvalues appear consecutively */ +/* with the eigenvalue having the positive imaginary part */ +/* first. */ + +/* VL (output) DOUBLE PRECISION array, dimension (LDVL,N) */ +/* If JOBVL = 'V', the left eigenvectors u(j) are stored one */ +/* after another in the columns of VL, in the same order */ +/* as their eigenvalues. */ +/* If JOBVL = 'N', VL is not referenced. */ +/* If the j-th eigenvalue is real, then u(j) = VL(:,j), */ +/* the j-th column of VL. */ +/* If the j-th and (j+1)-st eigenvalues form a complex */ +/* conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and */ +/* u(j+1) = VL(:,j) - i*VL(:,j+1). */ + +/* LDVL (input) INTEGER */ +/* The leading dimension of the array VL. LDVL >= 1; if */ +/* JOBVL = 'V', LDVL >= N. */ + +/* VR (output) DOUBLE PRECISION array, dimension (LDVR,N) */ +/* If JOBVR = 'V', the right eigenvectors v(j) are stored one */ +/* after another in the columns of VR, in the same order */ +/* as their eigenvalues. */ +/* If JOBVR = 'N', VR is not referenced. */ +/* If the j-th eigenvalue is real, then v(j) = VR(:,j), */ +/* the j-th column of VR. */ +/* If the j-th and (j+1)-st eigenvalues form a complex */ +/* conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and */ +/* v(j+1) = VR(:,j) - i*VR(:,j+1). */ + +/* LDVR (input) INTEGER */ +/* The leading dimension of the array VR. LDVR >= 1; if */ +/* JOBVR = 'V', LDVR >= N. */ + +/* WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) */ +/* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */ + +/* LWORK (input) INTEGER */ +/* The dimension of the array WORK. LWORK >= max(1,3*N), and */ +/* if JOBVL = 'V' or JOBVR = 'V', LWORK >= 4*N. For good */ +/* performance, LWORK must generally be larger. */ + +/* If LWORK = -1, then a workspace query is assumed; the routine */ +/* only calculates the optimal size of the WORK array, returns */ +/* this value as the first entry of the WORK array, and no error */ +/* message related to LWORK is issued by XERBLA. */ + +/* INFO (output) INTEGER */ +/* = 0: successful exit */ +/* < 0: if INFO = -i, the i-th argument had an illegal value. */ +/* > 0: if INFO = i, the QR algorithm failed to compute all the */ +/* eigenvalues, and no eigenvectors have been computed; */ +/* elements i+1:N of WR and WI contain eigenvalues which */ +/* have converged. */ + + // ------------------------------------------------------------------------------------ + + template < + typename T, + long NR1, long NR2, long NR3, long NR4, long NR5, + long NC1, long NC2, long NC3, long NC4, long NC5, + typename MM, + typename layout + > + int geev ( + const char jobvl, + const char jobvr, + matrix<T,NR1,NC1,MM,column_major_layout>& a, + matrix<T,NR2,NC2,MM,layout>& wr, + matrix<T,NR3,NC3,MM,layout>& wi, + matrix<T,NR4,NC4,MM,column_major_layout>& vl, + matrix<T,NR5,NC5,MM,column_major_layout>& vr + ) + { + matrix<T,0,1,MM,column_major_layout> work; + + const long n = a.nr(); + + wr.set_size(n,1); + wi.set_size(n,1); + + if (jobvl == 'V') + vl.set_size(n,n); + else + vl.set_size(NR4?NR4:1, NC4?NC4:1); + + if (jobvr == 'V') + vr.set_size(n,n); + else + vr.set_size(NR5?NR5:1, NC5?NC5:1); + + + // figure out how big the workspace needs to be. + T work_size = 1; + int info = binding::geev(jobvl, jobvr, n, &a(0,0), + a.nr(), &wr(0,0), &wi(0,0), &vl(0,0), + vl.nr(), &vr(0,0), vr.nr(), &work_size, + -1); + + if (info != 0) + return info; + + if (work.size() < work_size) + work.set_size(static_cast<long>(work_size), 1); + + // compute the actual decomposition + info = binding::geev(jobvl, jobvr, n, &a(0,0), + a.nr(), &wr(0,0), &wi(0,0), &vl(0,0), + vl.nr(), &vr(0,0), vr.nr(), &work(0,0), + work.size()); + + return info; + } + + // ------------------------------------------------------------------------------------ + + } + +} + +// ---------------------------------------------------------------------------------------- + +#endif // DLIB_LAPACk_GEEV_Hh_ + + diff --git a/ml/dlib/dlib/matrix/lapack/geqrf.h b/ml/dlib/dlib/matrix/lapack/geqrf.h new file mode 100644 index 000000000..c1f8fc050 --- /dev/null +++ b/ml/dlib/dlib/matrix/lapack/geqrf.h @@ -0,0 +1,168 @@ +// Copyright (C) 2010 Davis E. King (davis@dlib.net) +// License: Boost Software License See LICENSE.txt for the full license. +#ifndef DLIB_LAPACk_GEQRF_Hh_ +#define DLIB_LAPACk_GEQRF_Hh_ + +#include "fortran_id.h" +#include "../matrix.h" + +namespace dlib +{ + namespace lapack + { + namespace binding + { + extern "C" + { + void DLIB_FORTRAN_ID(dgeqrf) (integer *m, integer *n, double *a, integer * + lda, double *tau, double *work, integer *lwork, + integer *info); + + void DLIB_FORTRAN_ID(sgeqrf) (integer *m, integer *n, float *a, integer * + lda, float *tau, float *work, integer *lwork, + integer *info); + } + + inline int geqrf (integer m, integer n, double *a, integer lda, + double *tau, double *work, integer lwork) + { + integer info = 0; + DLIB_FORTRAN_ID(dgeqrf)(&m, &n, a, &lda, + tau, work, &lwork, &info); + return info; + } + + inline int geqrf (integer m, integer n, float *a, integer lda, + float *tau, float *work, integer lwork) + { + integer info = 0; + DLIB_FORTRAN_ID(sgeqrf)(&m, &n, a, &lda, + tau, work, &lwork, &info); + return info; + } + + + } + + // ------------------------------------------------------------------------------------ + +/* -- LAPACK routine (version 3.1) -- */ +/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ +/* November 2006 */ + +/* .. Scalar Arguments .. */ +/* .. */ +/* .. Array Arguments .. */ +/* .. */ + +/* Purpose */ +/* ======= */ + +/* DGEQRF computes a QR factorization of a real M-by-N matrix A: */ +/* A = Q * R. */ + +/* Arguments */ +/* ========= */ + +/* M (input) INTEGER */ +/* The number of rows of the matrix A. M >= 0. */ + +/* N (input) INTEGER */ +/* The number of columns of the matrix A. N >= 0. */ + +/* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) */ +/* On entry, the M-by-N matrix A. */ +/* On exit, the elements on and above the diagonal of the array */ +/* contain the min(M,N)-by-N upper trapezoidal matrix R (R is */ +/* upper triangular if m >= n); the elements below the diagonal, */ +/* with the array TAU, represent the orthogonal matrix Q as a */ +/* product of min(m,n) elementary reflectors (see Further */ +/* Details). */ + +/* LDA (input) INTEGER */ +/* The leading dimension of the array A. LDA >= max(1,M). */ + +/* TAU (output) DOUBLE PRECISION array, dimension (min(M,N)) */ +/* The scalar factors of the elementary reflectors (see Further */ +/* Details). */ + +/* WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) */ +/* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */ + +/* LWORK (input) INTEGER */ +/* The dimension of the array WORK. LWORK >= max(1,N). */ +/* For optimum performance LWORK >= N*NB, where NB is */ +/* the optimal blocksize. */ + +/* If LWORK = -1, then a workspace query is assumed; the routine */ +/* only calculates the optimal size of the WORK array, returns */ +/* this value as the first entry of the WORK array, and no error */ +/* message related to LWORK is issued by XERBLA. */ + +/* INFO (output) INTEGER */ +/* = 0: successful exit */ +/* < 0: if INFO = -i, the i-th argument had an illegal value */ + +/* Further Details */ +/* =============== */ + +/* The matrix Q is represented as a product of elementary reflectors */ + +/* Q = H(1) H(2) . . . H(k), where k = min(m,n). */ + +/* Each H(i) has the form */ + +/* H(i) = I - tau * v * v' */ + +/* where tau is a real scalar, and v is a real vector with */ +/* v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i), */ +/* and tau in TAU(i). */ + + + // ------------------------------------------------------------------------------------ + + template < + typename T, + long NR1, long NR2, + long NC1, long NC2, + typename MM + > + int geqrf ( + matrix<T,NR1,NC1,MM,column_major_layout>& a, + matrix<T,NR2,NC2,MM,column_major_layout>& tau + ) + { + matrix<T,0,1,MM,column_major_layout> work; + + tau.set_size(std::min(a.nr(), a.nc()), 1); + + // figure out how big the workspace needs to be. + T work_size = 1; + int info = binding::geqrf(a.nr(), a.nc(), &a(0,0), a.nr(), + &tau(0,0), &work_size, -1); + + if (info != 0) + return info; + + if (work.size() < work_size) + work.set_size(static_cast<long>(work_size), 1); + + // compute the actual decomposition + info = binding::geqrf(a.nr(), a.nc(), &a(0,0), a.nr(), + &tau(0,0), &work(0,0), work.size()); + + return info; + } + + // ------------------------------------------------------------------------------------ + + } + +} + +// ---------------------------------------------------------------------------------------- + +#endif // DLIB_LAPACk_GEQRF_Hh_ + + + diff --git a/ml/dlib/dlib/matrix/lapack/gesdd.h b/ml/dlib/dlib/matrix/lapack/gesdd.h new file mode 100644 index 000000000..e6b4d26e1 --- /dev/null +++ b/ml/dlib/dlib/matrix/lapack/gesdd.h @@ -0,0 +1,364 @@ +// Copyright (C) 2010 Davis E. King (davis@dlib.net) +// License: Boost Software License See LICENSE.txt for the full license. +#ifndef DLIB_LAPACk_SDD_Hh_ +#define DLIB_LAPACk_SDD_Hh_ + +#include "fortran_id.h" +#include "../matrix.h" + +namespace dlib +{ + namespace lapack + { + namespace binding + { + extern "C" + { + void DLIB_FORTRAN_ID(dgesdd) (char const* jobz, + const integer* m, const integer* n, double* a, const integer* lda, + double* s, double* u, const integer* ldu, + double* vt, const integer* ldvt, + double* work, const integer* lwork, integer* iwork, integer* info); + + void DLIB_FORTRAN_ID(sgesdd) (char const* jobz, + const integer* m, const integer* n, float* a, const integer* lda, + float* s, float* u, const integer* ldu, + float* vt, const integer* ldvt, + float* work, const integer* lwork, integer* iwork, integer* info); + + } + + inline integer gesdd (const char jobz, + const integer m, const integer n, double* a, const integer lda, + double* s, double* u, const integer ldu, + double* vt, const integer ldvt, + double* work, const integer lwork, integer* iwork) + { + integer info = 0; + DLIB_FORTRAN_ID(dgesdd)(&jobz, &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work, &lwork, iwork, &info); + return info; + } + + inline integer gesdd (const char jobz, + const integer m, const integer n, float* a, const integer lda, + float* s, float* u, const integer ldu, + float* vt, const integer ldvt, + float* work, const integer lwork, integer* iwork) + { + integer info = 0; + DLIB_FORTRAN_ID(sgesdd)(&jobz, &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work, &lwork, iwork, &info); + return info; + } + } + + // ------------------------------------------------------------------------------------ + +/* -- LAPACK driver routine (version 3.1) -- */ +/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ +/* November 2006 */ + +/* .. Scalar Arguments .. */ +/* .. */ +/* .. Array Arguments .. */ +/* .. */ + +/* Purpose */ +/* ======= */ + +/* DGESDD computes the singular value decomposition (SVD) of a real */ +/* M-by-N matrix A, optionally computing the left and right singular */ +/* vectors. If singular vectors are desired, it uses a */ +/* divide-and-conquer algorithm. */ + +/* The SVD is written */ + +/* A = U * SIGMA * transpose(V) */ + +/* where SIGMA is an M-by-N matrix which is zero except for its */ +/* min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and */ +/* V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA */ +/* are the singular values of A; they are real and non-negative, and */ +/* are returned in descending order. The first min(m,n) columns of */ +/* U and V are the left and right singular vectors of A. */ + +/* Note that the routine returns VT = V**T, not V. */ + +/* The divide and conquer algorithm makes very mild assumptions about */ +/* floating point arithmetic. It will work on machines with a guard */ +/* digit in add/subtract, or on those binary machines without guard */ +/* digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or */ +/* Cray-2. It could conceivably fail on hexadecimal or decimal machines */ +/* without guard digits, but we know of none. */ + +/* Arguments */ +/* ========= */ + +/* JOBZ (input) CHARACTER*1 */ +/* Specifies options for computing all or part of the matrix U: */ +/* = 'A': all M columns of U and all N rows of V**T are */ +/* returned in the arrays U and VT; */ +/* = 'S': the first min(M,N) columns of U and the first */ +/* min(M,N) rows of V**T are returned in the arrays U */ +/* and VT; */ +/* = 'O': If M >= N, the first N columns of U are overwritten */ +/* on the array A and all rows of V**T are returned in */ +/* the array VT; */ +/* otherwise, all columns of U are returned in the */ +/* array U and the first M rows of V**T are overwritten */ +/* in the array A; */ +/* = 'N': no columns of U or rows of V**T are computed. */ + +/* M (input) INTEGER */ +/* The number of rows of the input matrix A. M >= 0. */ + +/* N (input) INTEGER */ +/* The number of columns of the input matrix A. N >= 0. */ + +/* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) */ +/* On entry, the M-by-N matrix A. */ +/* On exit, */ +/* if JOBZ = 'O', A is overwritten with the first N columns */ +/* of U (the left singular vectors, stored */ +/* columnwise) if M >= N; */ +/* A is overwritten with the first M rows */ +/* of V**T (the right singular vectors, stored */ +/* rowwise) otherwise. */ +/* if JOBZ .ne. 'O', the contents of A are destroyed. */ + +/* LDA (input) INTEGER */ +/* The leading dimension of the array A. LDA >= max(1,M). */ + +/* S (output) DOUBLE PRECISION array, dimension (min(M,N)) */ +/* The singular values of A, sorted so that S(i) >= S(i+1). */ + +/* U (output) DOUBLE PRECISION array, dimension (LDU,UCOL) */ +/* UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N; */ +/* UCOL = min(M,N) if JOBZ = 'S'. */ +/* If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M */ +/* orthogonal matrix U; */ +/* if JOBZ = 'S', U contains the first min(M,N) columns of U */ +/* (the left singular vectors, stored columnwise); */ +/* if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced. */ + +/* LDU (input) INTEGER */ +/* The leading dimension of the array U. LDU >= 1; if */ +/* JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M. */ + +/* VT (output) DOUBLE PRECISION array, dimension (LDVT,N) */ +/* If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the */ +/* N-by-N orthogonal matrix V**T; */ +/* if JOBZ = 'S', VT contains the first min(M,N) rows of */ +/* V**T (the right singular vectors, stored rowwise); */ +/* if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced. */ + +/* LDVT (input) INTEGER */ +/* The leading dimension of the array VT. LDVT >= 1; if */ +/* JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N; */ +/* if JOBZ = 'S', LDVT >= min(M,N). */ + +/* WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) */ +/* On exit, if INFO = 0, WORK(1) returns the optimal LWORK; */ + +/* LWORK (input) INTEGER */ +/* The dimension of the array WORK. LWORK >= 1. */ +/* If JOBZ = 'N', */ +/* LWORK >= 3*min(M,N) + max(max(M,N),7*min(M,N)). */ +/* If JOBZ = 'O', */ +/* LWORK >= 3*min(M,N)*min(M,N) + */ +/* max(max(M,N),5*min(M,N)*min(M,N)+4*min(M,N)). */ +/* If JOBZ = 'S' or 'A' */ +/* LWORK >= 3*min(M,N)*min(M,N) + */ +/* max(max(M,N),4*min(M,N)*min(M,N)+4*min(M,N)). */ +/* For good performance, LWORK should generally be larger. */ +/* If LWORK = -1 but other input arguments are legal, WORK(1) */ +/* returns the optimal LWORK. */ + +/* IWORK (workspace) INTEGER array, dimension (8*min(M,N)) */ + +/* INFO (output) INTEGER */ +/* = 0: successful exit. */ +/* < 0: if INFO = -i, the i-th argument had an illegal value. */ +/* > 0: DBDSDC did not converge, updating process failed. */ + +/* Further Details */ +/* =============== */ + +/* Based on contributions by */ +/* Ming Gu and Huan Ren, Computer Science Division, University of */ +/* California at Berkeley, USA */ + + // ------------------------------------------------------------------------------------ + + template < + typename T, + long NR1, long NR2, long NR3, long NR4, + long NC1, long NC2, long NC3, long NC4, + typename MM + > + int gesdd ( + const char jobz, + matrix<T,NR1,NC1,MM,column_major_layout>& a, + matrix<T,NR2,NC2,MM,column_major_layout>& s, + matrix<T,NR3,NC3,MM,column_major_layout>& u, + matrix<T,NR4,NC4,MM,column_major_layout>& vt + ) + { + matrix<T,0,1,MM,column_major_layout> work; + matrix<integer,0,1,MM,column_major_layout> iwork; + + const long m = a.nr(); + const long n = a.nc(); + s.set_size(std::min(m,n), 1); + + // make sure the iwork memory block is big enough + if (iwork.size() < 8*std::min(m,n)) + iwork.set_size(8*std::min(m,n), 1); + + if (jobz == 'A') + { + u.set_size(m,m); + vt.set_size(n,n); + } + else if (jobz == 'S') + { + u.set_size(m, std::min(m,n)); + vt.set_size(std::min(m,n), n); + } + else if (jobz == 'O') + { + DLIB_CASSERT(false, "jobz == 'O' not supported"); + } + else + { + u.set_size(NR3?NR3:1, NC3?NC3:1); + vt.set_size(NR4?NR4:1, NC4?NC4:1); + } + + // figure out how big the workspace needs to be. + T work_size = 1; + int info = binding::gesdd(jobz, a.nr(), a.nc(), &a(0,0), a.nr(), + &s(0,0), &u(0,0), u.nr(), &vt(0,0), vt.nr(), + &work_size, -1, &iwork(0,0)); + + if (info != 0) + return info; + + // There is a bug in an older version of LAPACK in Debian etch + // that causes the gesdd to return the wrong value for work_size + // when jobz == 'N'. So verify the value of work_size. + if (jobz == 'N') + { + using std::min; + using std::max; + const T min_work_size = 3*min(m,n) + max(max(m,n),7*min(m,n)); + if (work_size < min_work_size) + work_size = min_work_size; + } + + if (work.size() < work_size) + work.set_size(static_cast<long>(work_size), 1); + + // compute the actual SVD + info = binding::gesdd(jobz, a.nr(), a.nc(), &a(0,0), a.nr(), + &s(0,0), &u(0,0), u.nr(), &vt(0,0), vt.nr(), + &work(0,0), work.size(), &iwork(0,0)); + + return info; + } + + // ------------------------------------------------------------------------------------ + + template < + typename T, + long NR1, long NR2, long NR3, long NR4, + long NC1, long NC2, long NC3, long NC4, + typename MM + > + int gesdd ( + const char jobz, + matrix<T,NR1,NC1,MM,row_major_layout>& a, + matrix<T,NR2,NC2,MM,row_major_layout>& s, + matrix<T,NR3,NC3,MM,row_major_layout>& u_, + matrix<T,NR4,NC4,MM,row_major_layout>& vt_ + ) + { + matrix<T,0,1,MM,row_major_layout> work; + matrix<integer,0,1,MM,row_major_layout> iwork; + + // Row major order matrices are transposed from LAPACK's point of view. + matrix<T,NR4,NC4,MM,row_major_layout>& u = vt_; + matrix<T,NR3,NC3,MM,row_major_layout>& vt = u_; + + + const long m = a.nc(); + const long n = a.nr(); + s.set_size(std::min(m,n), 1); + + // make sure the iwork memory block is big enough + if (iwork.size() < 8*std::min(m,n)) + iwork.set_size(8*std::min(m,n), 1); + + if (jobz == 'A') + { + u.set_size(m,m); + vt.set_size(n,n); + } + else if (jobz == 'S') + { + u.set_size(std::min(m,n), m); + vt.set_size(n, std::min(m,n)); + } + else if (jobz == 'O') + { + DLIB_CASSERT(false, "jobz == 'O' not supported"); + } + else + { + u.set_size(NR4?NR4:1, NC4?NC4:1); + vt.set_size(NR3?NR3:1, NC3?NC3:1); + } + + // figure out how big the workspace needs to be. + T work_size = 1; + int info = binding::gesdd(jobz, m, n, &a(0,0), a.nc(), + &s(0,0), &u(0,0), u.nc(), &vt(0,0), vt.nc(), + &work_size, -1, &iwork(0,0)); + + if (info != 0) + return info; + + // There is a bug in an older version of LAPACK in Debian etch + // that causes the gesdd to return the wrong value for work_size + // when jobz == 'N'. So verify the value of work_size. + if (jobz == 'N') + { + using std::min; + using std::max; + const T min_work_size = 3*min(m,n) + max(max(m,n),7*min(m,n)); + if (work_size < min_work_size) + work_size = min_work_size; + } + + + if (work.size() < work_size) + work.set_size(static_cast<long>(work_size), 1); + + // compute the actual SVD + info = binding::gesdd(jobz, m, n, &a(0,0), a.nc(), + &s(0,0), &u(0,0), u.nc(), &vt(0,0), vt.nc(), + &work(0,0), work.size(), &iwork(0,0)); + + return info; + } + + // ------------------------------------------------------------------------------------ + + } + +} + +// ---------------------------------------------------------------------------------------- + +#endif // DLIB_LAPACk_SDD_Hh_ + + diff --git a/ml/dlib/dlib/matrix/lapack/gesvd.h b/ml/dlib/dlib/matrix/lapack/gesvd.h new file mode 100644 index 000000000..e00654db6 --- /dev/null +++ b/ml/dlib/dlib/matrix/lapack/gesvd.h @@ -0,0 +1,323 @@ +// Copyright (C) 2010 Davis E. King (davis@dlib.net) +// License: Boost Software License See LICENSE.txt for the full license. +#ifndef DLIB_LAPACk_SVD_Hh_ +#define DLIB_LAPACk_SVD_Hh_ + +#include "fortran_id.h" +#include "../matrix.h" + +namespace dlib +{ + namespace lapack + { + namespace binding + { + extern "C" + { + void DLIB_FORTRAN_ID(dgesvd) (const char* jobu, const char* jobvt, + const integer* m, const integer* n, double* a, const integer* lda, + double* s, double* u, const integer* ldu, + double* vt, const integer* ldvt, + double* work, const integer* lwork, integer* info); + + void DLIB_FORTRAN_ID(sgesvd) (const char* jobu, const char* jobvt, + const integer* m, const integer* n, float* a, const integer* lda, + float* s, float* u, const integer* ldu, + float* vt, const integer* ldvt, + float* work, const integer* lwork, integer* info); + + } + + inline integer gesvd (const char jobu, const char jobvt, + const integer m, const integer n, double* a, const integer lda, + double* s, double* u, const integer ldu, + double* vt, const integer ldvt, + double* work, const integer lwork) + { + integer info = 0; + DLIB_FORTRAN_ID(dgesvd)(&jobu, &jobvt, &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work, &lwork, &info); + return info; + } + + inline integer gesvd (const char jobu, const char jobvt, + const integer m, const integer n, float* a, const integer lda, + float* s, float* u, const integer ldu, + float* vt, const integer ldvt, + float* work, const integer lwork) + { + integer info = 0; + DLIB_FORTRAN_ID(sgesvd)(&jobu, &jobvt, &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work, &lwork, &info); + return info; + } + + } + + // ------------------------------------------------------------------------------------ + +/* -- LAPACK driver routine (version 3.1) -- */ +/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ +/* November 2006 */ + +/* .. Scalar Arguments .. */ +/* .. */ +/* .. Array Arguments .. */ +/* .. */ + +/* Purpose */ +/* ======= */ + +/* DGESVD computes the singular value decomposition (SVD) of a real */ +/* M-by-N matrix A, optionally computing the left and/or right singular */ +/* vectors. The SVD is written */ + +/* A = U * SIGMA * transpose(V) */ + +/* where SIGMA is an M-by-N matrix which is zero except for its */ +/* min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and */ +/* V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA */ +/* are the singular values of A; they are real and non-negative, and */ +/* are returned in descending order. The first min(m,n) columns of */ +/* U and V are the left and right singular vectors of A. */ + +/* Note that the routine returns V**T, not V. */ + +/* Arguments */ +/* ========= */ + +/* JOBU (input) CHARACTER*1 */ +/* Specifies options for computing all or part of the matrix U: */ +/* = 'A': all M columns of U are returned in array U: */ +/* = 'S': the first min(m,n) columns of U (the left singular */ +/* vectors) are returned in the array U; */ +/* = 'O': the first min(m,n) columns of U (the left singular */ +/* vectors) are overwritten on the array A; */ +/* = 'N': no columns of U (no left singular vectors) are */ +/* computed. */ + +/* JOBVT (input) CHARACTER*1 */ +/* Specifies options for computing all or part of the matrix */ +/* V**T: */ +/* = 'A': all N rows of V**T are returned in the array VT; */ +/* = 'S': the first min(m,n) rows of V**T (the right singular */ +/* vectors) are returned in the array VT; */ +/* = 'O': the first min(m,n) rows of V**T (the right singular */ +/* vectors) are overwritten on the array A; */ +/* = 'N': no rows of V**T (no right singular vectors) are */ +/* computed. */ + +/* JOBVT and JOBU cannot both be 'O'. */ + +/* M (input) INTEGER */ +/* The number of rows of the input matrix A. M >= 0. */ + +/* N (input) INTEGER */ +/* The number of columns of the input matrix A. N >= 0. */ + +/* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) */ +/* On entry, the M-by-N matrix A. */ +/* On exit, */ +/* if JOBU = 'O', A is overwritten with the first min(m,n) */ +/* columns of U (the left singular vectors, */ +/* stored columnwise); */ +/* if JOBVT = 'O', A is overwritten with the first min(m,n) */ +/* rows of V**T (the right singular vectors, */ +/* stored rowwise); */ +/* if JOBU .ne. 'O' and JOBVT .ne. 'O', the contents of A */ +/* are destroyed. */ + +/* LDA (input) INTEGER */ +/* The leading dimension of the array A. LDA >= max(1,M). */ + +/* S (output) DOUBLE PRECISION array, dimension (min(M,N)) */ +/* The singular values of A, sorted so that S(i) >= S(i+1). */ + +/* U (output) DOUBLE PRECISION array, dimension (LDU,UCOL) */ +/* (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'. */ +/* If JOBU = 'A', U contains the M-by-M orthogonal matrix U; */ +/* if JOBU = 'S', U contains the first min(m,n) columns of U */ +/* (the left singular vectors, stored columnwise); */ +/* if JOBU = 'N' or 'O', U is not referenced. */ + +/* LDU (input) INTEGER */ +/* The leading dimension of the array U. LDU >= 1; if */ +/* JOBU = 'S' or 'A', LDU >= M. */ + +/* VT (output) DOUBLE PRECISION array, dimension (LDVT,N) */ +/* If JOBVT = 'A', VT contains the N-by-N orthogonal matrix */ +/* V**T; */ +/* if JOBVT = 'S', VT contains the first min(m,n) rows of */ +/* V**T (the right singular vectors, stored rowwise); */ +/* if JOBVT = 'N' or 'O', VT is not referenced. */ + +/* LDVT (input) INTEGER */ +/* The leading dimension of the array VT. LDVT >= 1; if */ +/* JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N). */ + +/* WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) */ +/* On exit, if INFO = 0, WORK(1) returns the optimal LWORK; */ +/* if INFO > 0, WORK(2:MIN(M,N)) contains the unconverged */ +/* superdiagonal elements of an upper bidiagonal matrix B */ +/* whose diagonal is in S (not necessarily sorted). B */ +/* satisfies A = U * B * VT, so it has the same singular values */ +/* as A, and singular vectors related by U and VT. */ + +/* LWORK (input) INTEGER */ +/* The dimension of the array WORK. */ +/* LWORK >= MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(M,N)). */ +/* For good performance, LWORK should generally be larger. */ + +/* If LWORK = -1, then a workspace query is assumed; the routine */ +/* only calculates the optimal size of the WORK array, returns */ +/* this value as the first entry of the WORK array, and no error */ +/* message related to LWORK is issued by XERBLA. */ + +/* INFO (output) INTEGER */ +/* = 0: successful exit. */ +/* < 0: if INFO = -i, the i-th argument had an illegal value. */ +/* > 0: if DBDSQR did not converge, INFO specifies how many */ +/* superdiagonals of an intermediate bidiagonal form B */ +/* did not converge to zero. See the description of WORK */ +/* above for details. */ + + // ------------------------------------------------------------------------------------ + + template < + typename T, + long NR1, long NR2, long NR3, long NR4, + long NC1, long NC2, long NC3, long NC4, + typename MM + > + int gesvd ( + const char jobu, + const char jobvt, + matrix<T,NR1,NC1,MM,column_major_layout>& a, + matrix<T,NR2,NC2,MM,column_major_layout>& s, + matrix<T,NR3,NC3,MM,column_major_layout>& u, + matrix<T,NR4,NC4,MM,column_major_layout>& vt + ) + { + matrix<T,0,1,MM,column_major_layout> work; + + const long m = a.nr(); + const long n = a.nc(); + s.set_size(std::min(m,n), 1); + + if (jobu == 'A') + u.set_size(m,m); + else if (jobu == 'S') + u.set_size(m, std::min(m,n)); + else + u.set_size(NR3?NR3:1, NC3?NC3:1); + + if (jobvt == 'A') + vt.set_size(n,n); + else if (jobvt == 'S') + vt.set_size(std::min(m,n), n); + else + vt.set_size(NR4?NR4:1, NC4?NC4:1); + + + if (jobu == 'O' || jobvt == 'O') + { + DLIB_CASSERT(false, "job == 'O' not supported"); + } + + + // figure out how big the workspace needs to be. + T work_size = 1; + int info = binding::gesvd(jobu, jobvt, a.nr(), a.nc(), &a(0,0), a.nr(), + &s(0,0), &u(0,0), u.nr(), &vt(0,0), vt.nr(), + &work_size, -1); + + if (info != 0) + return info; + + if (work.size() < work_size) + work.set_size(static_cast<long>(work_size), 1); + + // compute the actual SVD + info = binding::gesvd(jobu, jobvt, a.nr(), a.nc(), &a(0,0), a.nr(), + &s(0,0), &u(0,0), u.nr(), &vt(0,0), vt.nr(), + &work(0,0), work.size()); + + return info; + } + + // ------------------------------------------------------------------------------------ + + template < + typename T, + long NR1, long NR2, long NR3, long NR4, + long NC1, long NC2, long NC3, long NC4, + typename MM + > + int gesvd ( + char jobu, + char jobvt, + matrix<T,NR1,NC1,MM,row_major_layout>& a, + matrix<T,NR2,NC2,MM,row_major_layout>& s, + matrix<T,NR3,NC3,MM,row_major_layout>& u_, + matrix<T,NR4,NC4,MM,row_major_layout>& vt_ + ) + { + matrix<T,0,1,MM,row_major_layout> work; + + // Row major order matrices are transposed from LAPACK's point of view. + matrix<T,NR4,NC4,MM,row_major_layout>& u = vt_; + matrix<T,NR3,NC3,MM,row_major_layout>& vt = u_; + std::swap(jobu, jobvt); + + const long m = a.nc(); + const long n = a.nr(); + s.set_size(std::min(m,n), 1); + + if (jobu == 'A') + u.set_size(m,m); + else if (jobu == 'S') + u.set_size(std::min(m,n), m); + else + u.set_size(NR4?NR4:1, NC4?NC4:1); + + if (jobvt == 'A') + vt.set_size(n,n); + else if (jobvt == 'S') + vt.set_size(n, std::min(m,n)); + else + vt.set_size(NR3?NR3:1, NC3?NC3:1); + + if (jobu == 'O' || jobvt == 'O') + { + DLIB_CASSERT(false, "job == 'O' not supported"); + } + + + // figure out how big the workspace needs to be. + T work_size = 1; + int info = binding::gesvd(jobu, jobvt, m, n, &a(0,0), a.nc(), + &s(0,0), &u(0,0), u.nc(), &vt(0,0), vt.nc(), + &work_size, -1); + + if (info != 0) + return info; + + if (work.size() < work_size) + work.set_size(static_cast<long>(work_size), 1); + + // compute the actual SVD + info = binding::gesvd(jobu, jobvt, m, n, &a(0,0), a.nc(), + &s(0,0), &u(0,0), u.nc(), &vt(0,0), vt.nc(), + &work(0,0), work.size()); + + return info; + } + + // ------------------------------------------------------------------------------------ + + } + +} + +// ---------------------------------------------------------------------------------------- + +#endif // DLIB_LAPACk_SVD_Hh_ + diff --git a/ml/dlib/dlib/matrix/lapack/getrf.h b/ml/dlib/dlib/matrix/lapack/getrf.h new file mode 100644 index 000000000..a1f0b139d --- /dev/null +++ b/ml/dlib/dlib/matrix/lapack/getrf.h @@ -0,0 +1,132 @@ +// Copyright (C) 2010 Davis E. King (davis@dlib.net) +// License: Boost Software License See LICENSE.txt for the full license. +#ifndef DLIB_LAPACk_GETRF_Hh_ +#define DLIB_LAPACk_GETRF_Hh_ + +#include "fortran_id.h" +#include "../matrix.h" + +namespace dlib +{ + namespace lapack + { + namespace binding + { + extern "C" + { + void DLIB_FORTRAN_ID(dgetrf) (integer* m, integer *n, double *a, + integer* lda, integer *ipiv, integer *info); + + void DLIB_FORTRAN_ID(sgetrf) (integer* m, integer *n, float *a, + integer* lda, integer *ipiv, integer *info); + + } + + inline int getrf (integer m, integer n, double *a, + integer lda, integer *ipiv) + { + integer info = 0; + DLIB_FORTRAN_ID(dgetrf)(&m, &n, a, &lda, ipiv, &info); + return info; + } + + inline int getrf (integer m, integer n, float *a, + integer lda, integer *ipiv) + { + integer info = 0; + DLIB_FORTRAN_ID(sgetrf)(&m, &n, a, &lda, ipiv, &info); + return info; + } + + + } + + // ------------------------------------------------------------------------------------ + + +/* -- LAPACK routine (version 3.1) -- */ +/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ +/* November 2006 */ + +/* .. Scalar Arguments .. */ +/* .. */ +/* .. Array Arguments .. */ +/* .. */ + +/* Purpose */ +/* ======= */ + +/* DGETRF computes an LU factorization of a general M-by-N matrix A */ +/* using partial pivoting with row interchanges. */ + +/* The factorization has the form */ +/* A = P * L * U */ +/* where P is a permutation matrix, L is lower triangular with unit */ +/* diagonal elements (lower trapezoidal if m > n), and U is upper */ +/* triangular (upper trapezoidal if m < n). */ + +/* This is the right-looking Level 3 BLAS version of the algorithm. */ + +/* Arguments */ +/* ========= */ + +/* M (input) INTEGER */ +/* The number of rows of the matrix A. M >= 0. */ + +/* N (input) INTEGER */ +/* The number of columns of the matrix A. N >= 0. */ + +/* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) */ +/* On entry, the M-by-N matrix to be factored. */ +/* On exit, the factors L and U from the factorization */ +/* A = P*L*U; the unit diagonal elements of L are not stored. */ + +/* LDA (input) INTEGER */ +/* The leading dimension of the array A. LDA >= max(1,M). */ + +/* IPIV (output) INTEGER array, dimension (min(M,N)) */ +/* The pivot indices; for 1 <= i <= min(M,N), row i of the */ +/* matrix was interchanged with row IPIV(i). */ + +/* INFO (output) INTEGER */ +/* = 0: successful exit */ +/* < 0: if INFO = -i, the i-th argument had an illegal value */ +/* > 0: if INFO = i, U(i,i) is exactly zero. The factorization */ +/* has been completed, but the factor U is exactly */ +/* singular, and division by zero will occur if it is used */ +/* to solve a system of equations. */ + + + // ------------------------------------------------------------------------------------ + + template < + typename T, + long NR1, long NR2, + long NC1, long NC2, + typename MM, + typename layout + > + int getrf ( + matrix<T,NR1,NC1,MM,column_major_layout>& a, + matrix<integer,NR2,NC2,MM,layout>& ipiv + ) + { + const long m = a.nr(); + const long n = a.nc(); + + ipiv.set_size(std::min(m,n), 1); + + // compute the actual decomposition + return binding::getrf(m, n, &a(0,0), a.nr(), &ipiv(0,0)); + } + + // ------------------------------------------------------------------------------------ + + } + +} + +// ---------------------------------------------------------------------------------------- + +#endif // DLIB_LAPACk_GETRF_Hh_ + diff --git a/ml/dlib/dlib/matrix/lapack/ormqr.h b/ml/dlib/dlib/matrix/lapack/ormqr.h new file mode 100644 index 000000000..ab66ff4d2 --- /dev/null +++ b/ml/dlib/dlib/matrix/lapack/ormqr.h @@ -0,0 +1,224 @@ +// Copyright (C) 2010 Davis E. King (davis@dlib.net) +// License: Boost Software License See LICENSE.txt for the full license. +#ifndef DLIB_LAPACk_ORMQR_Hh_ +#define DLIB_LAPACk_ORMQR_Hh_ + +#include "fortran_id.h" +#include "../matrix.h" + +namespace dlib +{ + namespace lapack + { + namespace binding + { + extern "C" + { + void DLIB_FORTRAN_ID(dormqr) (char *side, char *trans, integer *m, integer *n, + integer *k, const double *a, integer *lda, const double *tau, + double * c_, integer *ldc, double *work, integer *lwork, + integer *info); + + void DLIB_FORTRAN_ID(sormqr) (char *side, char *trans, integer *m, integer *n, + integer *k, const float *a, integer *lda, const float *tau, + float * c_, integer *ldc, float *work, integer *lwork, + integer *info); + + } + + inline int ormqr (char side, char trans, integer m, integer n, + integer k, const double *a, integer lda, const double *tau, + double *c_, integer ldc, double *work, integer lwork) + { + integer info = 0; + DLIB_FORTRAN_ID(dormqr)(&side, &trans, &m, &n, + &k, a, &lda, tau, + c_, &ldc, work, &lwork, &info); + return info; + } + + inline int ormqr (char side, char trans, integer m, integer n, + integer k, const float *a, integer lda, const float *tau, + float *c_, integer ldc, float *work, integer lwork) + { + integer info = 0; + DLIB_FORTRAN_ID(sormqr)(&side, &trans, &m, &n, + &k, a, &lda, tau, + c_, &ldc, work, &lwork, &info); + return info; + } + + + + } + + // ------------------------------------------------------------------------------------ + +/* -- LAPACK routine (version 3.1) -- */ +/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ +/* November 2006 */ + +/* .. Scalar Arguments .. */ +/* .. */ +/* .. Array Arguments .. */ +/* .. */ + +/* Purpose */ +/* ======= */ + +/* DORMQR overwrites the general real M-by-N matrix C with */ + +/* SIDE = 'L' SIDE = 'R' */ +/* TRANS = 'N': Q * C C * Q */ +/* TRANS = 'T': Q**T * C C * Q**T */ + +/* where Q is a real orthogonal matrix defined as the product of k */ +/* elementary reflectors */ + +/* Q = H(1) H(2) . . . H(k) */ + +/* as returned by DGEQRF. Q is of order M if SIDE = 'L' and of order N */ +/* if SIDE = 'R'. */ + +/* Arguments */ +/* ========= */ + +/* SIDE (input) CHARACTER*1 */ +/* = 'L': apply Q or Q**T from the Left; */ +/* = 'R': apply Q or Q**T from the Right. */ + +/* TRANS (input) CHARACTER*1 */ +/* = 'N': No transpose, apply Q; */ +/* = 'T': Transpose, apply Q**T. */ + +/* M (input) INTEGER */ +/* The number of rows of the matrix C. M >= 0. */ + +/* N (input) INTEGER */ +/* The number of columns of the matrix C. N >= 0. */ + +/* K (input) INTEGER */ +/* The number of elementary reflectors whose product defines */ +/* the matrix Q. */ +/* If SIDE = 'L', M >= K >= 0; */ +/* if SIDE = 'R', N >= K >= 0. */ + +/* A (input) DOUBLE PRECISION array, dimension (LDA,K) */ +/* The i-th column must contain the vector which defines the */ +/* elementary reflector H(i), for i = 1,2,...,k, as returned by */ +/* DGEQRF in the first k columns of its array argument A. */ +/* A is modified by the routine but restored on exit. */ + +/* LDA (input) INTEGER */ +/* The leading dimension of the array A. */ +/* If SIDE = 'L', LDA >= max(1,M); */ +/* if SIDE = 'R', LDA >= max(1,N). */ + +/* TAU (input) DOUBLE PRECISION array, dimension (K) */ +/* TAU(i) must contain the scalar factor of the elementary */ +/* reflector H(i), as returned by DGEQRF. */ + +/* C (input/output) DOUBLE PRECISION array, dimension (LDC,N) */ +/* On entry, the M-by-N matrix C. */ +/* On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q. */ + +/* LDC (input) INTEGER */ +/* The leading dimension of the array C. LDC >= max(1,M). */ + +/* WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) */ +/* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */ + +/* LWORK (input) INTEGER */ +/* The dimension of the array WORK. */ +/* If SIDE = 'L', LWORK >= max(1,N); */ +/* if SIDE = 'R', LWORK >= max(1,M). */ +/* For optimum performance LWORK >= N*NB if SIDE = 'L', and */ +/* LWORK >= M*NB if SIDE = 'R', where NB is the optimal */ +/* blocksize. */ + +/* If LWORK = -1, then a workspace query is assumed; the routine */ +/* only calculates the optimal size of the WORK array, returns */ +/* this value as the first entry of the WORK array, and no error */ +/* message related to LWORK is issued by XERBLA. */ + +/* INFO (output) INTEGER */ +/* = 0: successful exit */ +/* < 0: if INFO = -i, the i-th argument had an illegal value */ + + // ------------------------------------------------------------------------------------ + + template < + typename T, + long NR1, long NR2, long NR3, + long NC1, long NC2, long NC3, + typename MM, + typename C_LAYOUT + > + int ormqr ( + char side, + char trans, + const matrix<T,NR1,NC1,MM,column_major_layout>& a, + const matrix<T,NR2,NC2,MM,column_major_layout>& tau, + matrix<T,NR3,NC3,MM,C_LAYOUT>& c + ) + { + long m = c.nr(); + long n = c.nc(); + const long k = a.nc(); + long ldc; + if (is_same_type<C_LAYOUT,column_major_layout>::value) + { + ldc = c.nr(); + } + else + { + // Since lapack expects c to be in column major layout we have to + // do something to make this work. Since a row major layout matrix + // will look just like a transposed C we can just swap a few things around. + + ldc = c.nc(); + swap(m,n); + + if (side == 'L') + side = 'R'; + else + side = 'L'; + + if (trans == 'T') + trans = 'N'; + else + trans = 'T'; + } + + matrix<T,0,1,MM,column_major_layout> work; + + // figure out how big the workspace needs to be. + T work_size = 1; + int info = binding::ormqr(side, trans, m, n, + k, &a(0,0), a.nr(), &tau(0,0), + &c(0,0), ldc, &work_size, -1); + + if (info != 0) + return info; + + if (work.size() < work_size) + work.set_size(static_cast<long>(work_size), 1); + + // compute the actual result + info = binding::ormqr(side, trans, m, n, + k, &a(0,0), a.nr(), &tau(0,0), + &c(0,0), ldc, &work(0,0), work.size()); + + return info; + } + + // ------------------------------------------------------------------------------------ + + } + +} + +// ---------------------------------------------------------------------------------------- + +#endif // DLIB_LAPACk_ORMQR_Hh_ + diff --git a/ml/dlib/dlib/matrix/lapack/pbtrf.h b/ml/dlib/dlib/matrix/lapack/pbtrf.h new file mode 100644 index 000000000..23bcc127b --- /dev/null +++ b/ml/dlib/dlib/matrix/lapack/pbtrf.h @@ -0,0 +1,178 @@ +// Copyright (C) 2010 Davis E. King (davis@dlib.net) +// License: Boost Software License See LICENSE.txt for the full license. +#ifndef DLIB_LAPACk_BDC_Hh_ +#define DLIB_LAPACk_BDC_Hh_ + +#include "fortran_id.h" +#include "../matrix.h" +namespace dlib +{ + namespace lapack + { + namespace binding + { + extern "C" + { + void DLIB_FORTRAN_ID(dpbtrf) (const char* uplo, const integer* n, const integer* kd, + double* ab, const integer* ldab, integer* info); + + void DLIB_FORTRAN_ID(spbtrf) (const char* uplo, const integer* n, const integer* kd, + float* ab, const integer* ldab, integer* info); + + } + + inline integer pbtrf (const char uplo, const integer n, const integer kd, + double* ab, const integer ldab) + { + integer info = 0; + DLIB_FORTRAN_ID(dpbtrf)(&uplo, &n, &kd, ab, &ldab, &info); + return info; + } + + inline integer pbtrf (const char uplo, const integer n, const integer kd, + float* ab, const integer ldab) + { + integer info = 0; + DLIB_FORTRAN_ID(spbtrf)(&uplo, &n, &kd, ab, &ldab, &info); + return info; + } + } + + // ------------------------------------------------------------------------------------ +/* DPBTRF(l) LAPACK routine (version 1.1) DPBTRF(l) + +NAME + DPBTRF - compute the Cholesky factorization of a real symmetric positive + definite band matrix A + +SYNOPSIS + + SUBROUTINE DPBTRF( UPLO, N, KD, AB, LDAB, INFO ) + + CHARACTER UPLO + + INTEGER INFO, KD, LDAB, N + + DOUBLE PRECISION AB( LDAB, * ) + +PURPOSE + DPBTRF computes the Cholesky factorization of a real symmetric positive + definite band matrix A. + + The factorization has the form + A = U**T * U, if UPLO = 'U', or + A = L * L**T, if UPLO = 'L', + where U is an upper triangular matrix and L is lower triangular. + +ARGUMENTS + + UPLO (input) CHARACTER*1 + = 'U': Upper triangle of A is stored; + = 'L': Lower triangle of A is stored. + + N (input) INTEGER + The order of the matrix A. N >= 0. + + KD (input) INTEGER + The number of superdiagonals of the matrix A if UPLO = 'U', or the + number of subdiagonals if UPLO = 'L'. KD >= 0. + + AB (input/output) DOUBLE PRECISION array, dimension (LDAB,N) + On entry, the upper or lower triangle of the symmetric band matrix + A, stored in the first KD+1 rows of the array. The j-th column of + A is stored in the j-th column of the array AB as follows: if UPLO + = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j; if UPLO = + 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd). + + On exit, if INFO = 0, the triangular factor U or L from the Chole- + sky factorization A = U**T*U or A = L*L**T of the band matrix A, in + the same storage format as A. + + LDAB (input) INTEGER + The leading dimension of the array AB. LDAB >= KD+1. + + INFO (output) INTEGER + = 0: successful exit + < 0: if INFO = -i, the i-th argument had an illegal value + > 0: if INFO = i, the leading minor of order i is not positive + definite, and the factorization could not be completed. + +FURTHER DETAILS + The band storage scheme is illustrated by the following example, when N = + 6, KD = 2, and UPLO = 'U': + + On entry: On exit: + + * * a13 a24 a35 a46 * * u13 u24 u35 u46 + * a12 a23 a34 a45 a56 * u12 u23 u34 u45 u56 + a11 a22 a33 a44 a55 a66 u11 u22 u33 u44 u55 u66 + + Similarly, if UPLO = 'L' the format of A is as follows: + + On entry: On exit: + + a11 a22 a33 a44 a55 a66 l11 l22 l33 l44 l55 l66 + a21 a32 a43 a54 a65 * l21 l32 l43 l54 l65 * + a31 a42 a53 a64 * * l31 l42 l53 l64 * * + + Array elements marked * are not used by the routine. + + Contributed by + Peter Mayes and Giuseppe Radicati, IBM ECSEC, Rome, March 23, 1989 */ + + // ------------------------------------------------------------------------------------ + + template < + typename T, + long NR1, long NC1, + typename MM + > + int pbtrf ( + char uplo, matrix<T,NR1,NC1,MM,column_major_layout>& ab + ) + { + const long ldab = ab.nr(); + const long n = ab.nc(); + const long kd = ldab - 1; // assume fully packed + + int info = binding::pbtrf(uplo, n, kd, &ab(0,0), ldab); + + return info; + } + + // ------------------------------------------------------------------------------------ + + + template < + typename T, + long NR1, long NC1, + typename MM + > + int pbtrf ( + char uplo, matrix<T,NR1,NC1,MM,row_major_layout>& ab + ) + { + const long ldab = ab.nr(); + const long n = ab.nc(); + const long kd = ldab - 1; // assume fully packed + + matrix<T,NC1,NR1,MM,row_major_layout> abt = trans(ab); + + int info = binding::pbtrf(uplo, n, kd, &abt(0,0), ldab); + + ab = trans(abt); + + return info; + } + + // ------------------------------------------------------------------------------------ + + } + +} + +// ---------------------------------------------------------------------------------------- + +#endif // DLIB_LAPACk_BDC_Hh_ + + diff --git a/ml/dlib/dlib/matrix/lapack/potrf.h b/ml/dlib/dlib/matrix/lapack/potrf.h new file mode 100644 index 000000000..b9d6a7cc8 --- /dev/null +++ b/ml/dlib/dlib/matrix/lapack/potrf.h @@ -0,0 +1,174 @@ +// Copyright (C) 2010 Davis E. King (davis@dlib.net) +// License: Boost Software License See LICENSE.txt for the full license. +#ifndef DLIB_LAPACk_POTRF_Hh_ +#define DLIB_LAPACk_POTRF_Hh_ + +#include "fortran_id.h" +#include "../matrix.h" + +namespace dlib +{ + namespace lapack + { + namespace binding + { + extern "C" + { + void DLIB_FORTRAN_ID(dpotrf) (char *uplo, integer *n, double *a, + integer* lda, integer *info); + + void DLIB_FORTRAN_ID(spotrf) (char *uplo, integer *n, float *a, + integer* lda, integer *info); + + } + + inline int potrf (char uplo, integer n, double *a, integer lda) + { + integer info = 0; + DLIB_FORTRAN_ID(dpotrf)(&uplo, &n, a, &lda, &info); + return info; + } + + inline int potrf (char uplo, integer n, float *a, integer lda) + { + integer info = 0; + DLIB_FORTRAN_ID(spotrf)(&uplo, &n, a, &lda, &info); + return info; + } + + } + + // ------------------------------------------------------------------------------------ + +/* -- LAPACK routine (version 3.1) -- */ +/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ +/* November 2006 */ + +/* .. Scalar Arguments .. */ +/* .. */ +/* .. Array Arguments .. */ +/* .. */ + +/* Purpose */ +/* ======= */ + +/* DPOTRF computes the Cholesky factorization of a real symmetric */ +/* positive definite matrix A. */ + +/* The factorization has the form */ +/* A = U**T * U, if UPLO = 'U', or */ +/* A = L * L**T, if UPLO = 'L', */ +/* where U is an upper triangular matrix and L is lower triangular. */ + +/* This is the block version of the algorithm, calling Level 3 BLAS. */ + +/* Arguments */ +/* ========= */ + +/* UPLO (input) CHARACTER*1 */ +/* = 'U': Upper triangle of A is stored; */ +/* = 'L': Lower triangle of A is stored. */ + +/* N (input) INTEGER */ +/* The order of the matrix A. N >= 0. */ + +/* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) */ +/* On entry, the symmetric matrix A. If UPLO = 'U', the leading */ +/* N-by-N upper triangular part of A contains the upper */ +/* triangular part of the matrix A, and the strictly lower */ +/* triangular part of A is not referenced. If UPLO = 'L', the */ +/* leading N-by-N lower triangular part of A contains the lower */ +/* triangular part of the matrix A, and the strictly upper */ +/* triangular part of A is not referenced. */ + +/* On exit, if INFO = 0, the factor U or L from the Cholesky */ +/* factorization A = U**T*U or A = L*L**T. */ + +/* LDA (input) INTEGER */ +/* The leading dimension of the array A. LDA >= max(1,N). */ + +/* INFO (output) INTEGER */ +/* = 0: successful exit */ +/* < 0: if INFO = -i, the i-th argument had an illegal value */ +/* > 0: if INFO = i, the leading minor of order i is not */ +/* positive definite, and the factorization could not be */ +/* completed. */ + + + // ------------------------------------------------------------------------------------ + + template < + typename T, + long NR1, + long NC1, + typename MM + > + int potrf ( + char uplo, + matrix<T,NR1,NC1,MM,column_major_layout>& a + ) + { + // compute the actual decomposition + int info = binding::potrf(uplo, a.nr(), &a(0,0), a.nr()); + + // If it fails part way though the factorization then make sure + // the end of the matrix gets properly initialized with zeros. + if (info > 0) + { + if (uplo == 'L') + set_colm(a, range(info-1, a.nc()-1)) = 0; + else + set_rowm(a, range(info-1, a.nr()-1)) = 0; + } + + return info; + } + + // ------------------------------------------------------------------------------------ + + template < + typename T, + long NR1, + long NC1, + typename MM + > + int potrf ( + char uplo, + matrix<T,NR1,NC1,MM,row_major_layout>& a + ) + { + // since we are working on a row major order matrix we need to ask + // LAPACK for the transpose of whatever the user asked for. + + if (uplo == 'L') + uplo = 'U'; + else + uplo = 'L'; + + // compute the actual decomposition + int info = binding::potrf(uplo, a.nr(), &a(0,0), a.nr()); + + // If it fails part way though the factorization then make sure + // the end of the matrix gets properly initialized with zeros. + if (info > 0) + { + if (uplo == 'U') + set_colm(a, range(info-1, a.nc()-1)) = 0; + else + set_rowm(a, range(info-1, a.nr()-1)) = 0; + } + + return info; + } + + // ------------------------------------------------------------------------------------ + + } + +} + +// ---------------------------------------------------------------------------------------- + +#endif // DLIB_LAPACk_POTRF_Hh_ + + diff --git a/ml/dlib/dlib/matrix/lapack/syev.h b/ml/dlib/dlib/matrix/lapack/syev.h new file mode 100644 index 000000000..0c9fd251a --- /dev/null +++ b/ml/dlib/dlib/matrix/lapack/syev.h @@ -0,0 +1,218 @@ +// Copyright (C) 2010 Davis E. King (davis@dlib.net) +// License: Boost Software License See LICENSE.txt for the full license. +#ifndef DLIB_LAPACk_EV_Hh_ +#define DLIB_LAPACk_EV_Hh_ + +#include "fortran_id.h" +#include "../matrix.h" + +namespace dlib +{ + namespace lapack + { + namespace binding + { + extern "C" + { + void DLIB_FORTRAN_ID(dsyev) (char *jobz, char *uplo, integer *n, double *a, + integer *lda, double *w, double *work, integer *lwork, + integer *info); + + void DLIB_FORTRAN_ID(ssyev) (char *jobz, char *uplo, integer *n, float *a, + integer *lda, float *w, float *work, integer *lwork, + integer *info); + + } + + inline int syev (char jobz, char uplo, integer n, double *a, + integer lda, double *w, double *work, integer lwork) + { + integer info = 0; + DLIB_FORTRAN_ID(dsyev)(&jobz, &uplo, &n, a, + &lda, w, work, &lwork, &info); + return info; + } + + inline int syev (char jobz, char uplo, integer n, float *a, + integer lda, float *w, float *work, integer lwork) + { + integer info = 0; + DLIB_FORTRAN_ID(ssyev)(&jobz, &uplo, &n, a, + &lda, w, work, &lwork, &info); + return info; + } + + + } + + // ------------------------------------------------------------------------------------ + +/* -- LAPACK driver routine (version 3.1) -- */ +/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ +/* November 2006 */ + +/* .. Scalar Arguments .. */ +/* .. */ +/* .. Array Arguments .. */ +/* .. */ + +/* Purpose */ +/* ======= */ + +/* DSYEV computes all eigenvalues and, optionally, eigenvectors of a */ +/* real symmetric matrix A. */ + +/* Arguments */ +/* ========= */ + +/* JOBZ (input) CHARACTER*1 */ +/* = 'N': Compute eigenvalues only; */ +/* = 'V': Compute eigenvalues and eigenvectors. */ + +/* UPLO (input) CHARACTER*1 */ +/* = 'U': Upper triangle of A is stored; */ +/* = 'L': Lower triangle of A is stored. */ + +/* N (input) INTEGER */ +/* The order of the matrix A. N >= 0. */ + +/* A (input/output) DOUBLE PRECISION array, dimension (LDA, N) */ +/* On entry, the symmetric matrix A. If UPLO = 'U', the */ +/* leading N-by-N upper triangular part of A contains the */ +/* upper triangular part of the matrix A. If UPLO = 'L', */ +/* the leading N-by-N lower triangular part of A contains */ +/* the lower triangular part of the matrix A. */ +/* On exit, if JOBZ = 'V', then if INFO = 0, A contains the */ +/* orthonormal eigenvectors of the matrix A. */ +/* If JOBZ = 'N', then on exit the lower triangle (if UPLO='L') */ +/* or the upper triangle (if UPLO='U') of A, including the */ +/* diagonal, is destroyed. */ + +/* LDA (input) INTEGER */ +/* The leading dimension of the array A. LDA >= max(1,N). */ + +/* W (output) DOUBLE PRECISION array, dimension (N) */ +/* If INFO = 0, the eigenvalues in ascending order. */ + +/* WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) */ +/* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */ + +/* LWORK (input) INTEGER */ +/* The length of the array WORK. LWORK >= max(1,3*N-1). */ +/* For optimal efficiency, LWORK >= (NB+2)*N, */ +/* where NB is the blocksize for DSYTRD returned by ILAENV. */ + +/* If LWORK = -1, then a workspace query is assumed; the routine */ +/* only calculates the optimal size of the WORK array, returns */ +/* this value as the first entry of the WORK array, and no error */ +/* message related to LWORK is issued by XERBLA. */ + +/* INFO (output) INTEGER */ +/* = 0: successful exit */ +/* < 0: if INFO = -i, the i-th argument had an illegal value */ +/* > 0: if INFO = i, the algorithm failed to converge; i */ +/* off-diagonal elements of an intermediate tridiagonal */ +/* form did not converge to zero. */ + + + // ------------------------------------------------------------------------------------ + + template < + typename T, + long NR1, long NR2, + long NC1, long NC2, + typename MM + > + int syev ( + const char jobz, + const char uplo, + matrix<T,NR1,NC1,MM,column_major_layout>& a, + matrix<T,NR2,NC2,MM,column_major_layout>& w + ) + { + matrix<T,0,1,MM,column_major_layout> work; + + const long n = a.nr(); + + w.set_size(n,1); + + + // figure out how big the workspace needs to be. + T work_size = 1; + int info = binding::syev(jobz, uplo, n, &a(0,0), + a.nr(), &w(0,0), &work_size, -1); + + if (info != 0) + return info; + + if (work.size() < work_size) + work.set_size(static_cast<long>(work_size), 1); + + // compute the actual decomposition + info = binding::syev(jobz, uplo, n, &a(0,0), + a.nr(), &w(0,0), &work(0,0), work.size()); + + return info; + } + + // ------------------------------------------------------------------------------------ + + template < + typename T, + long NR1, long NR2, + long NC1, long NC2, + typename MM + > + int syev ( + char jobz, + char uplo, + matrix<T,NR1,NC1,MM,row_major_layout>& a, + matrix<T,NR2,NC2,MM,row_major_layout>& w + ) + { + matrix<T,0,1,MM,row_major_layout> work; + + if (uplo == 'L') + uplo = 'U'; + else + uplo = 'L'; + + const long n = a.nr(); + + w.set_size(n,1); + + + // figure out how big the workspace needs to be. + T work_size = 1; + int info = binding::syev(jobz, uplo, n, &a(0,0), + a.nc(), &w(0,0), &work_size, -1); + + if (info != 0) + return info; + + if (work.size() < work_size) + work.set_size(static_cast<long>(work_size), 1); + + // compute the actual decomposition + info = binding::syev(jobz, uplo, n, &a(0,0), + a.nc(), &w(0,0), &work(0,0), work.size()); + + + a = trans(a); + + return info; + } + + // ------------------------------------------------------------------------------------ + + } + +} + +// ---------------------------------------------------------------------------------------- + +#endif // DLIB_LAPACk_EV_Hh_ + + + + diff --git a/ml/dlib/dlib/matrix/lapack/syevr.h b/ml/dlib/dlib/matrix/lapack/syevr.h new file mode 100644 index 000000000..65190b3d8 --- /dev/null +++ b/ml/dlib/dlib/matrix/lapack/syevr.h @@ -0,0 +1,445 @@ +// Copyright (C) 2010 Davis E. King (davis@dlib.net) +// License: Boost Software License See LICENSE.txt for the full license. +#ifndef DLIB_LAPACk_EVR_Hh_ +#define DLIB_LAPACk_EVR_Hh_ + +#include "fortran_id.h" +#include "../matrix.h" + +namespace dlib +{ + namespace lapack + { + namespace binding + { + extern "C" + { + void DLIB_FORTRAN_ID(dsyevr) (char *jobz, char *range, char *uplo, integer *n, + double *a, integer *lda, double *vl, double *vu, integer * il, + integer *iu, double *abstol, integer *m, double *w, + double *z_, integer *ldz, integer *isuppz, double *work, + integer *lwork, integer *iwork, integer *liwork, integer *info); + + void DLIB_FORTRAN_ID(ssyevr) (char *jobz, char *range, char *uplo, integer *n, + float *a, integer *lda, float *vl, float *vu, integer * il, + integer *iu, float *abstol, integer *m, float *w, + float *z_, integer *ldz, integer *isuppz, float *work, + integer *lwork, integer *iwork, integer *liwork, integer *info); + } + + inline int syevr (char jobz, char range, char uplo, integer n, + double* a, integer lda, double vl, double vu, integer il, + integer iu, double abstol, integer *m, double *w, + double *z, integer ldz, integer *isuppz, double *work, + integer lwork, integer *iwork, integer liwork) + { + integer info = 0; + DLIB_FORTRAN_ID(dsyevr)(&jobz, &range, &uplo, &n, + a, &lda, &vl, &vu, &il, + &iu, &abstol, m, w, + z, &ldz, isuppz, work, + &lwork, iwork, &liwork, &info); + return info; + } + + inline int syevr (char jobz, char range, char uplo, integer n, + float* a, integer lda, float vl, float vu, integer il, + integer iu, float abstol, integer *m, float *w, + float *z, integer ldz, integer *isuppz, float *work, + integer lwork, integer *iwork, integer liwork) + { + integer info = 0; + DLIB_FORTRAN_ID(ssyevr)(&jobz, &range, &uplo, &n, + a, &lda, &vl, &vu, &il, + &iu, &abstol, m, w, + z, &ldz, isuppz, work, + &lwork, iwork, &liwork, &info); + return info; + } + + } + + // ------------------------------------------------------------------------------------ + + /* + +* -- LAPACK driver routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + CHARACTER JOBZ, RANGE, UPLO + INTEGER IL, INFO, IU, LDA, LDZ, LIWORK, LWORK, M, N + DOUBLE PRECISION ABSTOL, VL, VU +* .. +* .. Array Arguments .. + INTEGER ISUPPZ( * ), IWORK( * ) + DOUBLE PRECISION A( LDA, * ), W( * ), WORK( * ), Z( LDZ, * ) +* .. +* +* Purpose +* ======= +* +* DSYEVR computes selected eigenvalues and, optionally, eigenvectors +* of a real symmetric matrix A. Eigenvalues and eigenvectors can be +* selected by specifying either a range of values or a range of +* indices for the desired eigenvalues. +* +* DSYEVR first reduces the matrix A to tridiagonal form T with a call +* to DSYTRD. Then, whenever possible, DSYEVR calls DSTEMR to compute +* the eigenspectrum using Relatively Robust Representations. DSTEMR +* computes eigenvalues by the dqds algorithm, while orthogonal +* eigenvectors are computed from various "good" L D L^T representations +* (also known as Relatively Robust Representations). Gram-Schmidt +* orthogonalization is avoided as far as possible. More specifically, +* the various steps of the algorithm are as follows. +* +* For each unreduced block (submatrix) of T, +* (a) Compute T - sigma I = L D L^T, so that L and D +* define all the wanted eigenvalues to high relative accuracy. +* This means that small relative changes in the entries of D and L +* cause only small relative changes in the eigenvalues and +* eigenvectors. The standard (unfactored) representation of the +* tridiagonal matrix T does not have this property in general. +* (b) Compute the eigenvalues to suitable accuracy. +* If the eigenvectors are desired, the algorithm attains full +* accuracy of the computed eigenvalues only right before +* the corresponding vectors have to be computed, see steps c) and d). +* (c) For each cluster of close eigenvalues, select a new +* shift close to the cluster, find a new factorization, and refine +* the shifted eigenvalues to suitable accuracy. +* (d) For each eigenvalue with a large enough relative separation compute +* the corresponding eigenvector by forming a rank revealing twisted +* factorization. Go back to (c) for any clusters that remain. +* +* The desired accuracy of the output can be specified by the input +* parameter ABSTOL. +* +* For more details, see DSTEMR's documentation and: +* - Inderjit S. Dhillon and Beresford N. Parlett: "Multiple representations +* to compute orthogonal eigenvectors of symmetric tridiagonal matrices," +* Linear Algebra and its Applications, 387(1), pp. 1-28, August 2004. +* - Inderjit Dhillon and Beresford Parlett: "Orthogonal Eigenvectors and +* Relative Gaps," SIAM Journal on Matrix Analysis and Applications, Vol. 25, +* 2004. Also LAPACK Working Note 154. +* - Inderjit Dhillon: "A new O(n^2) algorithm for the symmetric +* tridiagonal eigenvalue/eigenvector problem", +* Computer Science Division Technical Report No. UCB/CSD-97-971, +* UC Berkeley, May 1997. +* +* +* Note 1 : DSYEVR calls DSTEMR when the full spectrum is requested +* on machines which conform to the ieee-754 floating point standard. +* DSYEVR calls DSTEBZ and SSTEIN on non-ieee machines and +* when partial spectrum requests are made. +* +* Normal execution of DSTEMR may create NaNs and infinities and +* hence may abort due to a floating point exception in environments +* which do not handle NaNs and infinities in the ieee standard default +* manner. +* +* Arguments +* ========= +* +* JOBZ (input) CHARACTER*1 +* = 'N': Compute eigenvalues only; +* = 'V': Compute eigenvalues and eigenvectors. +* +* RANGE (input) CHARACTER*1 +* = 'A': all eigenvalues will be found. +* = 'V': all eigenvalues in the half-open interval (VL,VU] +* will be found. +* = 'I': the IL-th through IU-th eigenvalues will be found. +********** For RANGE = 'V' or 'I' and IU - IL < N - 1, DSTEBZ and +********** DSTEIN are called +* +* UPLO (input) CHARACTER*1 +* = 'U': Upper triangle of A is stored; +* = 'L': Lower triangle of A is stored. +* +* N (input) INTEGER +* The order of the matrix A. N >= 0. +* +* A (input/output) DOUBLE PRECISION array, dimension (LDA, N) +* On entry, the symmetric matrix A. If UPLO = 'U', the +* leading N-by-N upper triangular part of A contains the +* upper triangular part of the matrix A. If UPLO = 'L', +* the leading N-by-N lower triangular part of A contains +* the lower triangular part of the matrix A. +* On exit, the lower triangle (if UPLO='L') or the upper +* triangle (if UPLO='U') of A, including the diagonal, is +* destroyed. +* +* LDA (input) INTEGER +* The leading dimension of the array A. LDA >= max(1,N). +* +* VL (input) DOUBLE PRECISION +* VU (input) DOUBLE PRECISION +* If RANGE='V', the lower and upper bounds of the interval to +* be searched for eigenvalues. VL < VU. +* Not referenced if RANGE = 'A' or 'I'. +* +* IL (input) INTEGER +* IU (input) INTEGER +* If RANGE='I', the indices (in ascending order) of the +* smallest and largest eigenvalues to be returned. +* 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. +* Not referenced if RANGE = 'A' or 'V'. +* +* ABSTOL (input) DOUBLE PRECISION +* The absolute error tolerance for the eigenvalues. +* An approximate eigenvalue is accepted as converged +* when it is determined to lie in an interval [a,b] +* of width less than or equal to +* +* ABSTOL + EPS * max( |a|,|b| ) , +* +* where EPS is the machine precision. If ABSTOL is less than +* or equal to zero, then EPS*|T| will be used in its place, +* where |T| is the 1-norm of the tridiagonal matrix obtained +* by reducing A to tridiagonal form. +* +* See "Computing Small Singular Values of Bidiagonal Matrices +* with Guaranteed High Relative Accuracy," by Demmel and +* Kahan, LAPACK Working Note #3. +* +* If high relative accuracy is important, set ABSTOL to +* DLAMCH( 'Safe minimum' ). Doing so will guarantee that +* eigenvalues are computed to high relative accuracy when +* possible in future releases. The current code does not +* make any guarantees about high relative accuracy, but +* future releases will. See J. Barlow and J. Demmel, +* "Computing Accurate Eigensystems of Scaled Diagonally +* Dominant Matrices", LAPACK Working Note #7, for a discussion +* of which matrices define their eigenvalues to high relative +* accuracy. +* +* M (output) INTEGER +* The total number of eigenvalues found. 0 <= M <= N. +* If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1. +* +* W (output) DOUBLE PRECISION array, dimension (N) +* The first M elements contain the selected eigenvalues in +* ascending order. +* +* Z (output) DOUBLE PRECISION array, dimension (LDZ, max(1,M)) +* If JOBZ = 'V', then if INFO = 0, the first M columns of Z +* contain the orthonormal eigenvectors of the matrix A +* corresponding to the selected eigenvalues, with the i-th +* column of Z holding the eigenvector associated with W(i). +* If JOBZ = 'N', then Z is not referenced. +* Note: the user must ensure that at least max(1,M) columns are +* supplied in the array Z; if RANGE = 'V', the exact value of M +* is not known in advance and an upper bound must be used. +* Supplying N columns is always safe. +* +* LDZ (input) INTEGER +* The leading dimension of the array Z. LDZ >= 1, and if +* JOBZ = 'V', LDZ >= max(1,N). +* +* ISUPPZ (output) INTEGER array, dimension ( 2*max(1,M) ) +* The support of the eigenvectors in Z, i.e., the indices +* indicating the nonzero elements in Z. The i-th eigenvector +* is nonzero only in elements ISUPPZ( 2*i-1 ) through +* ISUPPZ( 2*i ). +********** Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1 +* +* WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) +* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. +* +* LWORK (input) INTEGER +* The dimension of the array WORK. LWORK >= max(1,26*N). +* For optimal efficiency, LWORK >= (NB+6)*N, +* where NB is the max of the blocksize for DSYTRD and DORMTR +* returned by ILAENV. +* +* If LWORK = -1, then a workspace query is assumed; the routine +* only calculates the optimal size of the WORK array, returns +* this value as the first entry of the WORK array, and no error +* message related to LWORK is issued by XERBLA. +* +* IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK)) +* On exit, if INFO = 0, IWORK(1) returns the optimal LWORK. +* +* LIWORK (input) INTEGER +* The dimension of the array IWORK. LIWORK >= max(1,10*N). +* +* If LIWORK = -1, then a workspace query is assumed; the +* routine only calculates the optimal size of the IWORK array, +* returns this value as the first entry of the IWORK array, and +* no error message related to LIWORK is issued by XERBLA. +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value +* > 0: Internal error +* +* Further Details +* =============== +* +* Based on contributions by +* Inderjit Dhillon, IBM Almaden, USA +* Osni Marques, LBNL/NERSC, USA +* Ken Stanley, Computer Science Division, University of +* California at Berkeley, USA +* Jason Riedy, Computer Science Division, University of +* California at Berkeley, USA +* +* ===================================================================== + + */ + + // ------------------------------------------------------------------------------------ + + template < + typename T, + long NR1, long NR2, long NR3, long NR4, + long NC1, long NC2, long NC3, long NC4, + typename MM + > + int syevr ( + const char jobz, + const char range, + const char uplo, + matrix<T,NR1,NC1,MM,column_major_layout>& a, + const double vl, + const double vu, + const integer il, + const integer iu, + const double abstol, + integer& num_eigenvalues_found, + matrix<T,NR2,NC2,MM,column_major_layout>& w, + matrix<T,NR3,NC3,MM,column_major_layout>& z, + matrix<integer,NR4,NC4,MM,column_major_layout>& isuppz + ) + { + matrix<T,0,1,MM,column_major_layout> work; + matrix<integer,0,1,MM,column_major_layout> iwork; + + const long n = a.nr(); + + w.set_size(n,1); + + isuppz.set_size(2*n, 1); + + if (jobz == 'V') + { + z.set_size(n,n); + } + else + { + z.set_size(NR3?NR3:1, NC3?NC3:1); + } + + // figure out how big the workspace needs to be. + T work_size = 1; + integer iwork_size = 1; + int info = binding::syevr(jobz, range, uplo, n, &a(0,0), + a.nr(), vl, vu, il, iu, abstol, &num_eigenvalues_found, + &w(0,0), &z(0,0), z.nr(), &isuppz(0,0), &work_size, -1, + &iwork_size, -1); + + if (info != 0) + return info; + + if (work.size() < work_size) + work.set_size(static_cast<long>(work_size), 1); + if (iwork.size() < iwork_size) + iwork.set_size(iwork_size, 1); + + // compute the actual decomposition + info = binding::syevr(jobz, range, uplo, n, &a(0,0), + a.nr(), vl, vu, il, iu, abstol, &num_eigenvalues_found, + &w(0,0), &z(0,0), z.nr(), &isuppz(0,0), &work(0,0), work.size(), + &iwork(0,0), iwork.size()); + + + return info; + } + + // ------------------------------------------------------------------------------------ + + template < + typename T, + long NR1, long NR2, long NR3, long NR4, + long NC1, long NC2, long NC3, long NC4, + typename MM + > + int syevr ( + const char jobz, + const char range, + char uplo, + matrix<T,NR1,NC1,MM,row_major_layout>& a, + const double vl, + const double vu, + const integer il, + const integer iu, + const double abstol, + integer& num_eigenvalues_found, + matrix<T,NR2,NC2,MM,row_major_layout>& w, + matrix<T,NR3,NC3,MM,row_major_layout>& z, + matrix<integer,NR4,NC4,MM,row_major_layout>& isuppz + ) + { + matrix<T,0,1,MM,row_major_layout> work; + matrix<integer,0,1,MM,row_major_layout> iwork; + + if (uplo == 'L') + uplo = 'U'; + else + uplo = 'L'; + + const long n = a.nr(); + + w.set_size(n,1); + + isuppz.set_size(2*n, 1); + + if (jobz == 'V') + { + z.set_size(n,n); + } + else + { + z.set_size(NR3?NR3:1, NC3?NC3:1); + } + + // figure out how big the workspace needs to be. + T work_size = 1; + integer iwork_size = 1; + int info = binding::syevr(jobz, range, uplo, n, &a(0,0), + a.nc(), vl, vu, il, iu, abstol, &num_eigenvalues_found, + &w(0,0), &z(0,0), z.nc(), &isuppz(0,0), &work_size, -1, + &iwork_size, -1); + + if (info != 0) + return info; + + if (work.size() < work_size) + work.set_size(static_cast<long>(work_size), 1); + if (iwork.size() < iwork_size) + iwork.set_size(iwork_size, 1); + + // compute the actual decomposition + info = binding::syevr(jobz, range, uplo, n, &a(0,0), + a.nc(), vl, vu, il, iu, abstol, &num_eigenvalues_found, + &w(0,0), &z(0,0), z.nc(), &isuppz(0,0), &work(0,0), work.size(), + &iwork(0,0), iwork.size()); + + z = trans(z); + + return info; + } + + // ------------------------------------------------------------------------------------ + + } + +} + +// ---------------------------------------------------------------------------------------- + +#endif // DLIB_LAPACk_EVR_Hh_ + + + |