/usr/include/dlib/matrix/lapack/gees.h is in libdlib-dev 18.18-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 | // 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 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 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_
|