Skip to content

D
dlib
Commit 16e910c3 authored by Davis King's avatar Davis King
Browse files
Options

Added a bunch of LAPACK bindings. 

--HG--
extra : convert_revision : svn%3Afdd8eb12-d10e-0410-9acb-85c331704f74/trunk%403817
parent fb542d70
Hide whitespace changes
Inline Side-by-side
Showing with 2186 additions and 0 deletions
dlib/matrix/lapack/fortran_id.h 0 → 100644
View file @ 16e910c3
#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(__GNUC__) || defined(__ICC) || defined(__sgi) || defined(__COMO__) || defined(__KCC)
#define DLIB_BIND_FORTRAN_LOWERCASE_UNDERSCORE
#elif defined(__IBMCPP__) || defined(_MSC_VER)
#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__)
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
dlib/matrix/lapack/gees.h 0 → 100644
View file @ 16e910c3
#ifndef DLIB_LAPACk_ES_H__
#define DLIB_LAPACk_ES_H__
#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 nil = 0;
integer sdim = 0;
DLIB_FORTRAN_ID(dgees)(&jobvs, &sort, fnil, &n,
a, &lda, &sdim, wr,
wi, vs, &ldvs, work,
&lwork, &nil, &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 nil = 0;
integer sdim = 0;
DLIB_FORTRAN_ID(sgees)(&jobvs, &sort, fnil, &n,
a, &lda, &sdim, wr,
wi, vs, &ldvs, work,
&lwork, &nil, &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 NR5,
long NC1, long NC2, long NC3, long NC4, long NC5,
typename MM
>
int gees (
const char jobz,
matrix<T,NR1,NC1,MM,column_major_layout>& a,
matrix<T,NR2,NC2,MM,column_major_layout>& wr,
matrix<T,NR3,NC3,MM,column_major_layout>& wi,
matrix<T,NR4,NC4,MM,column_major_layout>& vs,
matrix<T,NR5,NC5,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(1,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(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_H__
dlib/matrix/lapack/geev.h 0 → 100644
View file @ 16e910c3
#ifndef DLIB_LAPACk_GEEV_H__
#define DLIB_LAPACk_GEEV_H__
#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 NR6,
long NC1, long NC2, long NC3, long NC4, long NC5, long NC6,
typename MM
>
int geev (
const char jobvl,
const char jobvr,
matrix<T,NR1,NC1,MM,column_major_layout>& a,
matrix<T,NR2,NC2,MM,column_major_layout>& wr,
matrix<T,NR3,NC3,MM,column_major_layout>& wi,
matrix<T,NR4,NC4,MM,column_major_layout>& vl,
matrix<T,NR5,NC5,MM,column_major_layout>& vr,
matrix<T,NR6,NC6,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(1,1);
if (jobvr == 'V')
vr.set_size(n,n);
else
vr.set_size(1,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(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_H__
dlib/matrix/lapack/geqrf.h 0 → 100644
View file @ 16e910c3
#ifndef DLIB_LAPACk_GEQRF_H__
#define DLIB_LAPACk_GEQRF_H__
#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 NR3,
long NC1, long NC2, long NC3,
typename MM
>
int geqrf (
matrix<T,NR1,NC1,MM,column_major_layout>& a,
matrix<T,NR2,NC2,MM,column_major_layout>& tau,
matrix<T,NR3,NC3,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(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_H__
dlib/matrix/lapack/gesdd.h 0 → 100644
View file @ 16e910c3
#ifndef DLIB_LAPACk_SDD_H__
#define DLIB_LAPACk_SDD_H__
#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 NR5, long NR6,
long NC1, long NC2, long NC3, long NC4, long NC5, long NC6,
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,NR5,NC5,MM,column_major_layout>& work,
matrix<integer,NR6,NC6,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
{
u.set_size(1,1);
vt.set_size(1,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;
if (work.size() < work_size)
work.set_size(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 NR5, long NR6,
long NC1, long NC2, long NC3, long NC4, long NC5, long NC6,
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>& vt,
matrix<T,NR4,NC4,MM,row_major_layout>& u,
matrix<T,NR5,NC5,MM,row_major_layout>& work,
matrix<integer,NR6,NC6,MM,row_major_layout>& iwork
)
{
// Row major order matrices are transposed from LAPACK's point of view.
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
{
u.set_size(1,1);
vt.set_size(1,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;
if (work.size() < work_size)
work.set_size(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_H__
dlib/matrix/lapack/gesvd.h 0 → 100644
View file @ 16e910c3
#ifndef DLIB_LAPACk_SVD_H__
#define DLIB_LAPACk_SVD_H__
#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 NR5,
long NC1, long NC2, long NC3, long NC4, long NC5,
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,NR5,NC5,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(1,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(1,1);
// 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(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 NR5,
long NC1, long NC2, long NC3, long NC4, long NC5,
typename MM
>
int gesvd (
const char jobu,
const 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>& vt,
matrix<T,NR4,NC4,MM,row_major_layout>& u,
matrix<T,NR5,NC5,MM,row_major_layout>& work
)
{
// Row major order matrices are transposed from LAPACK's point of view.
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(1,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(1,1);
// 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(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_H__
dlib/matrix/lapack/getrf.h 0 → 100644
View file @ 16e910c3
#ifndef DLIB_LAPACk_GETRF_H__
#define DLIB_LAPACk_GETRF_H__
#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<long,NR2,NC2,MM,layout>& ipiv
)
{
const long m = a.nr();
const long n = a.nc();
matrix<integer,NR2,NC2,MM,column_major_layout> ipiv_temp(std::min(m,n), 1);
// compute the actual decomposition
int info = binding::getrf(m, n, &a(0,0), a.nr(), &ipiv_temp(0,0));
// Turn the P vector into a more useful form. This way we will have the identity
// a == rowm(L*U, ipiv). The permutation vector that comes out of LAPACK is somewhat
// different.
ipiv = trans(range(0, a.nr()-1));
for (long i = ipiv_temp.size()-1; i >= 0; --i)
{
// -1 because FORTRAN is indexed starting with 1 instead of 0
std::swap(ipiv(i), ipiv(ipiv_temp(i)-1));
}
return info;
}
// ------------------------------------------------------------------------------------
}
}
// ----------------------------------------------------------------------------------------
#endif // DLIB_LAPACk_GETRF_H__
dlib/matrix/lapack/potrf.h 0 → 100644
View file @ 16e910c3
#ifndef DLIB_LAPACk_POTRF_H__
#define DLIB_LAPACk_POTRF_H__
#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 getrf (
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 getrf (
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_H__
dlib/matrix/lapack/syev.h 0 → 100644
View file @ 16e910c3
#ifndef DLIB_LAPACk_EV_H__
#define DLIB_LAPACk_EV_H__
#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 NR3,
long NC1, long NC2, long NC3,
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,NR3,NC3,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(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;
}
// ------------------------------------------------------------------------------------
}
}
// ----------------------------------------------------------------------------------------
#endif // DLIB_LAPACk_EV_H__
dlib/matrix/lapack/syevr.h 0 → 100644
View file @ 16e910c3
#ifndef DLIB_LAPACk_EVR_H__
#define DLIB_LAPACk_EVR_H__
#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 NR5, long NR6,
long NC1, long NC2, long NC3, long NC4, long NC5, long NC6,
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,NR5,NC5,MM,column_major_layout>& work,
matrix<integer,NR6,NC6,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(1,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(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;
}
// ------------------------------------------------------------------------------------
}
}
// ----------------------------------------------------------------------------------------
#endif // DLIB_LAPACk_EVR_H__
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment