/*
* Copyright (c) 2002-2006 Samit Basu
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*/
#include "EigenDecompose.hpp"
#include "LAPACK.hpp"
#include <stdlib.h>
#include <stdio.h>
#include "Malloc.hpp"
#define MAX(a,b) ((a) > (b) ? (a) : (b))
void floatEigenDecomposeSymmetric(int n, float *v, float *d, float *a,
bool eigenvectors) {
// SUBROUTINE SSYEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO )
//* Purpose
//* =======
//*
//* SSYEV 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.
char JOBZ;
if (eigenvectors)
JOBZ = 'V';
else
JOBZ = 'N';
//*
//* UPLO (input) CHARACTER*1
//* = 'U': Upper triangle of A is stored;
//* = 'L': Lower triangle of A is stored.
//*
char UPLO = 'U';
//* N (input) INTEGER
//* The order of the matrix A. N >= 0.
//*
int N = n;
//* A (input/output) REAL 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.
//*
float *Ain = a;
//* LDA (input) INTEGER
//* The leading dimension of the array A. LDA >= max(1,N).
//*
int LDA = n;
//* W (output) REAL array, dimension (N)
//* If INFO = 0, the eigenvalues in ascending order.
//*
//* WORK (workspace/output) REAL array, dimension (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 SSYTRD 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.
//
int INFO;
float WORKSZE;
int LWORK;
LWORK = -1;
ssyev_(&JOBZ, &UPLO, &N, Ain, &LDA, d, &WORKSZE, &LWORK, &INFO);
LWORK = (int) WORKSZE;
float *WORK = (float*) Malloc(LWORK*sizeof(float));
ssyev_(&JOBZ, &UPLO, &N, Ain, &LDA, d, WORK, &LWORK, &INFO);
Free(WORK);
if (eigenvectors)
memcpy(v,a,n*n*sizeof(float));
}
void floatEigenDecompose(int n, float *v, float *d, float *a,
bool eigenvectors, bool balance) {
// SUBROUTINE SGEEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, WR, WI,
// $ VL, LDVL, VR, LDVR, ILO, IHI, SCALE, ABNRM,
// $ RCONDE, RCONDV, WORK, LWORK, IWORK, INFO )
//* Purpose
//* =======
//*
//* SGEEVX computes for an N-by-N real nonsymmetric matrix A, the
//* eigenvalues and, optionally, the left and/or right eigenvectors.
//*
//* Optionally also, it computes a balancing transformation to improve
//* the conditioning of the eigenvalues and eigenvectors (ILO, IHI,
//* SCALE, and ABNRM), reciprocal condition numbers for the eigenvalues
//* (RCONDE), and reciprocal condition numbers for the right
//* eigenvectors (RCONDV).
//*
//* 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.
//*
//* Balancing a matrix means permuting the rows and columns to make it
//* more nearly upper triangular, and applying a diagonal similarity
//* transformation D * A * D**(-1), where D is a diagonal matrix, to
//* make its rows and columns closer in norm and the condition numbers
//* of its eigenvalues and eigenvectors smaller. The computed
//* reciprocal condition numbers correspond to the balanced matrix.
//* Permuting rows and columns will not change the condition numbers
//* (in exact arithmetic) but diagonal scaling will. For further
//* explanation of balancing, see section 4.10.2 of the LAPACK
//* Users' Guide.
//*
//* Arguments
//* =========
//*
//* BALANC (input) CHARACTER*1
//* Indicates how the input matrix should be diagonally scaled
//* and/or permuted to improve the conditioning of its
//* eigenvalues.
//* = 'N': Do not diagonally scale or permute;
//* = 'P': Perform permutations to make the matrix more nearly
//* upper triangular. Do not diagonally scale;
//* = 'S': Diagonally scale the matrix, i.e. replace A by
//* D*A*D**(-1), where D is a diagonal matrix chosen
//* to make the rows and columns of A more equal in
//* norm. Do not permute;
//* = 'B': Both diagonally scale and permute A.
//*
//* Computed reciprocal condition numbers will be for the matrix
//* after balancing and/or permuting. Permuting does not change
//* condition numbers (in exact arithmetic), but balancing does.
//*
char BALANC;
if (balance)
BALANC = 'B';
else
BALANC = 'N';
//* JOBVL (input) CHARACTER*1
//* = 'N': left eigenvectors of A are not computed;
//* = 'V': left eigenvectors of A are computed.
//* If SENSE = 'E' or 'B', JOBVL must = 'V'.
char JOBVL = 'N';
//* JOBVR (input) CHARACTER*1
//* = 'N': right eigenvectors of A are not computed;
//* = 'V': right eigenvectors of A are computed.
//* If SENSE = 'E' or 'B', JOBVR must = 'V'.
char JOBVR;
if (eigenvectors)
JOBVR = 'V';
else
JOBVR = 'N';
//* SENSE (input) CHARACTER*1
//* Determines which reciprocal condition numbers are computed.
//* = 'N': None are computed;
//* = 'E': Computed for eigenvalues only;
//* = 'V': Computed for right eigenvectors only;
//* = 'B': Computed for eigenvalues and right eigenvectors.
//*
//* If SENSE = 'E' or 'B', both left and right eigenvectors
//* must also be computed (JOBVL = 'V' and JOBVR = 'V').
char SENSE = 'N';
//* N (input) INTEGER
//* The order of the matrix A. N >= 0.
int N = n;
//* A (input/output) REAL array, dimension (LDA,N)
//* On entry, the N-by-N matrix A.
//* On exit, A has been overwritten. If JOBVL = 'V' or
//* JOBVR = 'V', A contains the real Schur form of the balanced
//* version of the input matrix A.
float *Ain = a;
//* LDA (input) INTEGER
//* The leading dimension of the array A. LDA >= max(1,N).
int LDA = n;
//* WR (output) REAL array, dimension (N)
//* WI (output) REAL array, dimension (N)
//* WR and WI contain the real and imaginary parts,
//* respectively, of the computed eigenvalues. Complex
//* conjugate pairs of eigenvalues will appear consecutively
//* with the eigenvalue having the positive imaginary part
//* first.
float *WR = (float*) Malloc(n*sizeof(float));
float *WI = (float*) Malloc(n*sizeof(float));
//* VL (output) REAL 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).
float *VL = NULL;
//* LDVL (input) INTEGER
//* The leading dimension of the array VL. LDVL >= 1; if
//* JOBVL = 'V', LDVL >= N.
int LDVL = 1;
//* VR (output) REAL 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).
float *VR = v;
//* LDVR (input) INTEGER
//* The leading dimension of the array VR. LDVR >= 1, and if
//* JOBVR = 'V', LDVR >= N.
int LDVR = n;
//* ILO,IHI (output) INTEGER
//* ILO and IHI are integer values determined when A was
//* balanced. The balanced A(i,j) = 0 if I > J and
//* J = 1,...,ILO-1 or I = IHI+1,...,N.
int ILO;
int IHI;
//* SCALE (output) REAL array, dimension (N)
//* Details of the permutations and scaling factors applied
//* when balancing A. If P(j) is the index of the row and column
//* interchanged with row and column j, and D(j) is the scaling
//* factor applied to row and column j, then
//* SCALE(J) = P(J), for J = 1,...,ILO-1
//* = D(J), for J = ILO,...,IHI
//* = P(J) for J = IHI+1,...,N.
//* The order in which the interchanges are made is N to IHI+1,
//* then 1 to ILO-1.
float *SCALE = (float*) Malloc(n*sizeof(float));
//* ABNRM (output) REAL
//* The one-norm of the balanced matrix (the maximum
//* of the sum of absolute values of elements of any column).
float ABNRM;
//* RCONDE (output) REAL array, dimension (N)
//* RCONDE(j) is the reciprocal condition number of the j-th
//* eigenvalue.
float *RCONDE = (float*) Malloc(n*sizeof(float));
//* RCONDV (output) REAL array, dimension (N)
//* RCONDV(j) is the reciprocal condition number of the j-th
//* right eigenvector.
float *RCONDV = (float*) Malloc(n*sizeof(float));
//* WORK (workspace/output) REAL array, dimension (LWORK)
//* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
// double *WORK = (double*) Malloc(250*sizeof(double));
//*
//* LWORK (input) INTEGER
//* The dimension of the array WORK. If SENSE = 'N' or 'E',
//* LWORK >= max(1,2*N), and if JOBVL = 'V' or JOBVR = 'V',
//* LWORK >= 3*N. If SENSE = 'V' or 'B', LWORK >= N*(N+6).
//* For good performance, LWORK must generally be larger.
//*
//* If LWORK = -1, a workspace query is assumed. The optimal
//* size for the WORK array is calculated and stored in WORK(1),
//* and no other work except argument checking is performed.
int LWORK;
//* IWORK (workspace) INTEGER array, dimension (2*N-2)
//* If SENSE = 'N' or 'E', not referenced.
int *IWORK = (int*) Malloc((2*n-2)*sizeof(int));
//* 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 or condition numbers
//* have been computed; elements 1:ILO-1 and i+1:N of WR
//* and WI contain eigenvalues which have converged.
int INFO;
float WORKSZE;
LWORK = -1;
sgeevx_( &BALANC, &JOBVL, &JOBVR, &SENSE, &N, Ain, &LDA, WR, WI,
VL, &LDVL, VR, &LDVR, &ILO, &IHI, SCALE, &ABNRM,
RCONDE, RCONDV, &WORKSZE, &LWORK, IWORK, &INFO );
LWORK = (int) WORKSZE;
float *WORK = (float*) Malloc(LWORK*sizeof(float));
sgeevx_( &BALANC, &JOBVL, &JOBVR, &SENSE, &N, Ain, &LDA, WR, WI,
VL, &LDVL, VR, &LDVR, &ILO, &IHI, SCALE, &ABNRM,
RCONDE, RCONDV, WORK, &LWORK, IWORK, &INFO );
for (int i=0;i<N;i++) {
d[2*i] = WR[i];
d[2*i+1] = WI[i];
}
Free(WORK);
Free(WR);
Free(WI);
Free(SCALE);
Free(RCONDE);
Free(RCONDV);
Free(IWORK);
}
void floatGenEigenDecompose(int n, float *v, float *d, float *a,
float *b, bool eigenvectors) {
// SUBROUTINE SGGEV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI,
// $ BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO )
//* Purpose
//* =======
//*
//* SGGEV computes for a pair of N-by-N real nonsymmetric matrices (A,B)
//* the generalized eigenvalues, and optionally, the left and/or right
//* generalized eigenvectors.
//*
//* A generalized eigenvalue for a pair of matrices (A,B) is a scalar
//* lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
//* singular. It is usually represented as the pair (alpha,beta), as
//* there is a reasonable interpretation for beta=0, and even for both
//* being zero.
//*
//* The right eigenvector v(j) corresponding to the eigenvalue lambda(j)
//* of (A,B) satisfies
//*
//* A * v(j) = lambda(j) * B * v(j).
//*
//* The left eigenvector u(j) corresponding to the eigenvalue lambda(j)
//* of (A,B) satisfies
//*
//* u(j)**H * A = lambda(j) * u(j)**H * B .
//*
//* where u(j)**H is the conjugate-transpose of u(j).
//*
//*
//* Arguments
//* =========
//*
//* JOBVL (input) CHARACTER*1
//* = 'N': do not compute the left generalized eigenvectors;
//* = 'V': compute the left generalized eigenvectors.
//*
char JOBVL = 'N';
//* JOBVR (input) CHARACTER*1
//* = 'N': do not compute the right generalized eigenvectors;
//* = 'V': compute the right generalized eigenvectors.
//*
char JOBVR;
if (eigenvectors)
JOBVR = 'V';
else
JOBVR = 'N';
//* N (input) INTEGER
//* The order of the matrices A, B, VL, and VR. N >= 0.
//*
int N = n;
//* A (input/output) REAL array, dimension (LDA, N)
//* On entry, the matrix A in the pair (A,B).
//* On exit, A has been overwritten.
//*
float *A = a;
//* LDA (input) INTEGER
//* The leading dimension of A. LDA >= max(1,N).
//*
int LDA = n;
//* B (input/output) REAL array, dimension (LDB, N)
//* On entry, the matrix B in the pair (A,B).
//* On exit, B has been overwritten.
//*
float *B = b;
//* LDB (input) INTEGER
//* The leading dimension of B. LDB >= max(1,N).
//*
int LDB = n;
//* ALPHAR (output) REAL array, dimension (N)
//* ALPHAI (output) REAL array, dimension (N)
//* BETA (output) REAL array, dimension (N)
//* On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
//* be the generalized eigenvalues. If ALPHAI(j) is zero, then
//* the j-th eigenvalue is real; if positive, then the j-th and
//* (j+1)-st eigenvalues are a complex conjugate pair, with
//* ALPHAI(j+1) negative.
//*
float *ALPHAR = (float*) Malloc(n*sizeof(float));
float *ALPHAI = (float*) Malloc(n*sizeof(float));
float *BETA = (float*) Malloc(n*sizeof(float));
//* Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j)
//* may easily over- or underflow, and BETA(j) may even be zero.
//* Thus, the user should avoid naively computing the ratio
//* alpha/beta. However, ALPHAR and ALPHAI will be always less
//* than and usually comparable with norm(A) in magnitude, and
//* BETA always less than and usually comparable with norm(B).
//*
//* VL (output) REAL 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 the j-th eigenvalue is real, then
//* u(j) = VL(:,j), the j-th column of VL. If the j-th and
//* (j+1)-th 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).
//* Each eigenvector will be scaled so the largest component have
//* abs(real part)+abs(imag. part)=1.
//* Not referenced if JOBVL = 'N'.
//*
float *VL = NULL;
//* LDVL (input) INTEGER
//* The leading dimension of the matrix VL. LDVL >= 1, and
//* if JOBVL = 'V', LDVL >= N.
//*
int LDVL = 1;
//* VR (output) REAL 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 the j-th eigenvalue is real, then
//* v(j) = VR(:,j), the j-th column of VR. If the j-th and
//* (j+1)-th 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).
//* Each eigenvector will be scaled so the largest component have
//* abs(real part)+abs(imag. part)=1.
//* Not referenced if JOBVR = 'N'.
//*
float *VR = v;
//* LDVR (input) INTEGER
//* The leading dimension of the matrix VR. LDVR >= 1, and
//* if JOBVR = 'V', LDVR >= N.
//*
int LDVR = n;
//* WORK (workspace/output) REAL array, dimension (LWORK)
//* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
//*
//* LWORK (input) INTEGER
//* The dimension of the array WORK. LWORK >= max(1,8*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.
//* = 1,...,N:
//* The QZ iteration failed. No eigenvectors have been
//* calculated, but ALPHAR(j), ALPHAI(j), and BETA(j)
//* should be correct for j=INFO+1,...,N.
//* > N: =N+1: other than QZ iteration failed in SHGEQZ.
//* =N+2: error return from STGEVC.
//
float WORKSZE;
int LWORK = -1;
int INFO;
sggev_( &JOBVL, &JOBVR, &N, A, &LDA, B, &LDB, ALPHAR, ALPHAI,
BETA, VL, &LDVL, VR, &LDVR, &WORKSZE, &LWORK, &INFO );
LWORK = (int) WORKSZE;
float *WORK = (float*) Malloc(LWORK*sizeof(float));
sggev_( &JOBVL, &JOBVR, &N, A, &LDA, B, &LDB, ALPHAR, ALPHAI,
BETA, VL, &LDVL, VR, &LDVR, WORK, &LWORK, &INFO );
int i;
for (i=0;i<n;i++) {
d[2*i] = ALPHAR[i]/BETA[i];
d[2*i+1] = ALPHAI[i]/BETA[i];
}
Free(ALPHAR);
Free(BETA);
Free(ALPHAI);
Free(WORK);
}
bool floatGenEigenDecomposeSymmetric(int n, float *v, float *d,
float *a, float *b, bool eigenvectors) {
// SUBROUTINE SSYGV( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK,
// $ LWORK, INFO )
//* Purpose
//* =======
//*
//* SSYGV computes all the eigenvalues, and optionally, the eigenvectors
//* of a real generalized symmetric-definite eigenproblem, of the form
//* A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x.
//* Here A and B are assumed to be symmetric and B is also
//* positive definite.
//*
//* Arguments
//* =========
//*
//* ITYPE (input) INTEGER
//* Specifies the problem type to be solved:
//* = 1: A*x = (lambda)*B*x
//* = 2: A*B*x = (lambda)*x
//* = 3: B*A*x = (lambda)*x
//*
int ITYPE = 1;
//* JOBZ (input) CHARACTER*1
//* = 'N': Compute eigenvalues only;
//* = 'V': Compute eigenvalues and eigenvectors.
//*
char JOBZ;
if (eigenvectors)
JOBZ = 'V';
else
JOBZ = 'N';
//* UPLO (input) CHARACTER*1
//* = 'U': Upper triangles of A and B are stored;
//* = 'L': Lower triangles of A and B are stored.
//*
char UPLO = 'U';
//* N (input) INTEGER
//* The order of the matrices A and B. N >= 0.
//*
int N = n;
//* A (input/output) REAL 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
//* matrix Z of eigenvectors. The eigenvectors are normalized
//* as follows:
//* if ITYPE = 1 or 2, Z**T*B*Z = I;
//* if ITYPE = 3, Z**T*inv(B)*Z = I.
//* If JOBZ = 'N', then on exit the upper triangle (if UPLO='U')
//* or the lower triangle (if UPLO='L') of A, including the
//* diagonal, is destroyed.
//*
float *A = a;
//* LDA (input) INTEGER
//* The leading dimension of the array A. LDA >= max(1,N).
//*
int LDA = n;
//* B (input/output) REAL array, dimension (LDB, N)
//* On entry, the symmetric positive definite matrix B.
//* If UPLO = 'U', the leading N-by-N upper triangular part of B
//* contains the upper triangular part of the matrix B.
//* If UPLO = 'L', the leading N-by-N lower triangular part of B
//* contains the lower triangular part of the matrix B.
//*
//* On exit, if INFO <= N, the part of B containing the matrix is
//* overwritten by the triangular factor U or L from the Cholesky
//* factorization B = U**T*U or B = L*L**T.
//*
float *B = b;
//* LDB (input) INTEGER
//* The leading dimension of the array B. LDB >= max(1,N).
//*
int LDB = n;
//* W (output) REAL array, dimension (N)
//* If INFO = 0, the eigenvalues in ascending order.
//*
float *W = d;
//* WORK (workspace/output) REAL array, dimension (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 SSYTRD 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: SPOTRF or SSYEV returned an error code:
//* <= N: if INFO = i, SSYEV failed to converge;
//* i off-diagonal elements of an intermediate
//* tridiagonal form did not converge to zero;
//* > N: if INFO = N + i, for 1 <= i <= N, then the leading
//* minor of order i of B is not positive definite.
//* The factorization of B could not be completed and
//* no eigenvalues or eigenvectors were computed.
float WORKSIZE;
int LWORK = -1;
int INFO;
ssygv_( &ITYPE, &JOBZ, &UPLO, &N, A, &LDA, B, &LDB, W, &WORKSIZE,
&LWORK, &INFO );
LWORK = (int) WORKSIZE;
float *WORK = (float*) Malloc(LWORK*sizeof(float));
ssygv_( &ITYPE, &JOBZ, &UPLO, &N, A, &LDA, B, &LDB, W, WORK,
&LWORK, &INFO );
Free(WORK);
if (INFO>N) return false;
if (eigenvectors)
memcpy(v,a,n*n*sizeof(float));
return true;
}
void doubleEigenDecomposeSymmetric(int n, double *v, double *d,
double *a, bool eigenvectors) {
// SUBROUTINE DSYEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO )
//* 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.
char JOBZ;
if (eigenvectors)
JOBZ = 'V';
else
JOBZ = 'N';
//*
//* UPLO (input) CHARACTER*1
//* = 'U': Upper triangle of A is stored;
//* = 'L': Lower triangle of A is stored.
//*
char UPLO = 'U';
//* N (input) INTEGER
//* The order of the matrix A. N >= 0.
//*
int N = n;
//* 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.
//*
double *Ain = a;
//* LDA (input) INTEGER
//* The leading dimension of the array A. LDA >= max(1,N).
//*
int LDA = n;
//* W (output) DOUBLE PRECISION array, dimension (N)
//* If INFO = 0, the eigenvalues in ascending order.
//*
//* WORK (workspace/output) DOUBLE PRECISION array, dimension (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 SSYTRD 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.
//
int INFO;
double WORKSZE;
int LWORK;
LWORK = -1;
dsyev_(&JOBZ, &UPLO, &N, Ain, &LDA, d, &WORKSZE, &LWORK, &INFO);
LWORK = (int) WORKSZE;
double *WORK = (double*) Malloc(LWORK*sizeof(double));
dsyev_(&JOBZ, &UPLO, &N, Ain, &LDA, d, WORK, &LWORK, &INFO);
Free(WORK);
if (eigenvectors)
memcpy(v,a,n*n*sizeof(double));
}
void doubleEigenDecompose(int n, double *v, double *d, double *a,
bool eigenvectors, bool balance) {
// SUBROUTINE DGEEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, WR, WI,
// $ VL, LDVL, VR, LDVR, ILO, IHI, SCALE, ABNRM,
// $ RCONDE, RCONDV, WORK, LWORK, IWORK, INFO )
//* Purpose
//* =======
//*
//* DGEEVX computes for an N-by-N real nonsymmetric matrix A, the
//* eigenvalues and, optionally, the left and/or right eigenvectors.
//*
//* Optionally also, it computes a balancing transformation to improve
//* the conditioning of the eigenvalues and eigenvectors (ILO, IHI,
//* SCALE, and ABNRM), reciprocal condition numbers for the eigenvalues
//* (RCONDE), and reciprocal condition numbers for the right
//* eigenvectors (RCONDV).
//*
//* 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.
//*
//* Balancing a matrix means permuting the rows and columns to make it
//* more nearly upper triangular, and applying a diagonal similarity
//* transformation D * A * D**(-1), where D is a diagonal matrix, to
//* make its rows and columns closer in norm and the condition numbers
//* of its eigenvalues and eigenvectors smaller. The computed
//* reciprocal condition numbers correspond to the balanced matrix.
//* Permuting rows and columns will not change the condition numbers
//* (in exact arithmetic) but diagonal scaling will. For further
//* explanation of balancing, see section 4.10.2 of the LAPACK
//* Users' Guide.
//*
//* Arguments
//* =========
//*
//* BALANC (input) CHARACTER*1
//* Indicates how the input matrix should be diagonally scaled
//* and/or permuted to improve the conditioning of its
//* eigenvalues.
//* = 'N': Do not diagonally scale or permute;
//* = 'P': Perform permutations to make the matrix more nearly
//* upper triangular. Do not diagonally scale;
//* = 'S': Diagonally scale the matrix, i.e. replace A by
//* D*A*D**(-1), where D is a diagonal matrix chosen
//* to make the rows and columns of A more equal in
//* norm. Do not permute;
//* = 'B': Both diagonally scale and permute A.
//*
//* Computed reciprocal condition numbers will be for the matrix
//* after balancing and/or permuting. Permuting does not change
//* condition numbers (in exact arithmetic), but balancing does.
//*
char BALANC;
if (balance)
BALANC = 'B';
else
BALANC = 'N';
//* JOBVL (input) CHARACTER*1
//* = 'N': left eigenvectors of A are not computed;
//* = 'V': left eigenvectors of A are computed.
//* If SENSE = 'E' or 'B', JOBVL must = 'V'.
char JOBVL = 'N';
//* JOBVR (input) CHARACTER*1
//* = 'N': right eigenvectors of A are not computed;
//* = 'V': right eigenvectors of A are computed.
//* If SENSE = 'E' or 'B', JOBVR must = 'V'.
char JOBVR;
if (eigenvectors)
JOBVR = 'V';
else
JOBVR = 'N';
//* SENSE (input) CHARACTER*1
//* Determines which reciprocal condition numbers are computed.
//* = 'N': None are computed;
//* = 'E': Computed for eigenvalues only;
//* = 'V': Computed for right eigenvectors only;
//* = 'B': Computed for eigenvalues and right eigenvectors.
//*
//* If SENSE = 'E' or 'B', both left and right eigenvectors
//* must also be computed (JOBVL = 'V' and JOBVR = 'V').
char SENSE = 'N';
//* N (input) INTEGER
//* The order of the matrix A. N >= 0.
int N = n;
//* A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
//* On entry, the N-by-N matrix A.
//* On exit, A has been overwritten. If JOBVL = 'V' or
//* JOBVR = 'V', A contains the real Schur form of the balanced
//* version of the input matrix A.
double *Ain = a;
//* LDA (input) INTEGER
//* The leading dimension of the array A. LDA >= max(1,N).
int LDA = 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 will appear consecutively
//* with the eigenvalue having the positive imaginary part
//* first.
double *WR = (double*) Malloc(n*sizeof(double));
double *WI = (double*) Malloc(n*sizeof(double));
//* 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).
double *VL = NULL;
//* LDVL (input) INTEGER
//* The leading dimension of the array VL. LDVL >= 1; if
//* JOBVL = 'V', LDVL >= N.
int LDVL = 1;
//* 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).
double *VR = v;
//* LDVR (input) INTEGER
//* The leading dimension of the array VR. LDVR >= 1, and if
//* JOBVR = 'V', LDVR >= N.
int LDVR = n;
//* ILO,IHI (output) INTEGER
//* ILO and IHI are integer values determined when A was
//* balanced. The balanced A(i,j) = 0 if I > J and
//* J = 1,...,ILO-1 or I = IHI+1,...,N.
int ILO;
int IHI;
//* SCALE (output) DOUBLE PRECISION array, dimension (N)
//* Details of the permutations and scaling factors applied
//* when balancing A. If P(j) is the index of the row and column
//* interchanged with row and column j, and D(j) is the scaling
//* factor applied to row and column j, then
//* SCALE(J) = P(J), for J = 1,...,ILO-1
//* = D(J), for J = ILO,...,IHI
//* = P(J) for J = IHI+1,...,N.
//* The order in which the interchanges are made is N to IHI+1,
//* then 1 to ILO-1.
double *SCALE = (double*) Malloc(n*sizeof(double));
//* ABNRM (output) DOUBLE PRECISION
//* The one-norm of the balanced matrix (the maximum
//* of the sum of absolute values of elements of any column).
double ABNRM;
//* RCONDE (output) DOUBLE PRECISION array, dimension (N)
//* RCONDE(j) is the reciprocal condition number of the j-th
//* eigenvalue.
double *RCONDE = (double*) Malloc(n*sizeof(double));
//* RCONDV (output) DOUBLE PRECISION array, dimension (N)
//* RCONDV(j) is the reciprocal condition number of the j-th
//* right eigenvector.
double *RCONDV = (double*) Malloc(n*sizeof(double));
//* WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
//* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
// double *WORK = (double*) Malloc(250*sizeof(double));
//*
//* LWORK (input) INTEGER
//* The dimension of the array WORK. If SENSE = 'N' or 'E',
//* LWORK >= max(1,2*N), and if JOBVL = 'V' or JOBVR = 'V',
//* LWORK >= 3*N. If SENSE = 'V' or 'B', LWORK >= N*(N+6).
//* For good performance, LWORK must generally be larger.
//*
//* If LWORK = -1, a workspace query is assumed. The optimal
//* size for the WORK array is calculated and stored in WORK(1),
//* and no other work except argument checking is performed.
int maxN = (N < 6) ? 6 : N;
int LWORK = maxN*maxN*2;
//* IWORK (workspace) INTEGER array, dimension (2*N-2)
//* If SENSE = 'N' or 'E', not referenced.
int *IWORK = (int*) Malloc((2*n-2)*sizeof(int));
//* 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 or condition numbers
//* have been computed; elements 1:ILO-1 and i+1:N of WR
//* and WI contain eigenvalues which have converged.
int INFO;
double WORKSZE;
// LWORK = -1;
// dgeevx_( &BALANC, &JOBVL, &JOBVR, &SENSE, &N, Ain, &LDA, WR, WI,
// VL, &LDVL, VR, &LDVR, &ILO, &IHI, SCALE, &ABNRM,
// RCONDE, RCONDV, &WORKSZE, &LWORK, IWORK, &INFO );
// LWORK = (int) WORKSZE;
double *WORK = (double*) Malloc(LWORK*sizeof(double));
dgeevx_( &BALANC, &JOBVL, &JOBVR, &SENSE, &N, Ain, &LDA, WR, WI,
VL, &LDVL, VR, &LDVR, &ILO, &IHI, SCALE, &ABNRM,
RCONDE, RCONDV, WORK, &LWORK, IWORK, &INFO );
for (int i=0;i<N;i++) {
d[2*i] = WR[i];
d[2*i+1] = WI[i];
}
Free(WORK);
Free(WR);
Free(WI);
Free(SCALE);
Free(RCONDE);
Free(RCONDV);
Free(IWORK);
}
void doubleGenEigenDecompose(int n, double *v, double *d, double *a,
double *b, bool eigenvectors) {
char JOBVL = 'N';
char JOBVR;
if (eigenvectors)
JOBVR = 'V';
else
JOBVR = 'N';
int N = n;
double *A = a;
int LDA = n;
double *B = b;
int LDB = n;
double *ALPHAR = (double*) Malloc(n*sizeof(double));
double *ALPHAI = (double*) Malloc(n*sizeof(double));
double *BETA = (double*) Malloc(n*sizeof(double));
double *VL = NULL;
int LDVL = 1;
double *VR = v;
int LDVR = n;
double WORKSZE;
int LWORK = -1;
int INFO;
dggev_( &JOBVL, &JOBVR, &N, A, &LDA, B, &LDB, ALPHAR, ALPHAI,
BETA, VL, &LDVL, VR, &LDVR, &WORKSZE, &LWORK, &INFO );
LWORK = (int) WORKSZE;
double *WORK = (double*) Malloc(LWORK*sizeof(double));
dggev_( &JOBVL, &JOBVR, &N, A, &LDA, B, &LDB, ALPHAR, ALPHAI,
BETA, VL, &LDVL, VR, &LDVR, WORK, &LWORK, &INFO );
int i;
for (i=0;i<n;i++) {
d[2*i] = ALPHAR[i]/BETA[i];
d[2*i+1] = ALPHAI[i]/BETA[i];
}
Free(ALPHAR);
Free(BETA);
Free(ALPHAI);
Free(WORK);
}
bool doubleGenEigenDecomposeSymmetric(int n, double *v, double *d,
double *a, double *b, bool eigenvectors) {
int ITYPE = 1;
char JOBZ;
if (eigenvectors)
JOBZ = 'V';
else
JOBZ = 'N';
char UPLO = 'U';
int N = n;
double *A = a;
int LDA = n;
double *B = b;
int LDB = n;
double *W = d;
double WORKSIZE;
int LWORK = -1;
int INFO;
dsygv_( &ITYPE, &JOBZ, &UPLO, &N, A, &LDA, B, &LDB, W, &WORKSIZE,
&LWORK, &INFO );
LWORK = (int) WORKSIZE;
double *WORK = (double*) Malloc(LWORK*sizeof(double));
dsygv_( &ITYPE, &JOBZ, &UPLO, &N, A, &LDA, B, &LDB, W, WORK,
&LWORK, &INFO );
Free(WORK);
if (INFO>N) return false;
if (eigenvectors)
memcpy(v,a,n*n*sizeof(double));
return true;
}
void complexEigenDecomposeSymmetric(int n, float *v, float *d,
float *a, bool eigenvectors) {
// SUBROUTINE CHEEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK,
// $ INFO )
//* Purpose
//* =======
//*
//* CHEEV computes all eigenvalues and, optionally, eigenvectors of a
//* complex Hermitian matrix A.
//*
//* Arguments
//* =========
//*
//* JOBZ (input) CHARACTER*1
//* = 'N': Compute eigenvalues only;
//* = 'V': Compute eigenvalues and eigenvectors.
//*
char JOBZ;
if (eigenvectors)
JOBZ = 'V';
else
JOBZ = 'N';
//* UPLO (input) CHARACTER*1
//* = 'U': Upper triangle of A is stored;
//* = 'L': Lower triangle of A is stored.
//*
char UPLO = 'U';
//* N (input) INTEGER
//* The order of the matrix A. N >= 0.
//*
int N = n;
//* A (input/output) COMPLEX array, dimension (LDA, N)
//* On entry, the Hermitian 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).
//*
int LDA = n;
//* W (output) REAL array, dimension (N)
//* If INFO = 0, the eigenvalues in ascending order.
//*
//* WORK (workspace/output) COMPLEX array, dimension (LWORK)
//* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
//*
//* LWORK (input) INTEGER
//* The length of the array WORK. LWORK >= max(1,2*N-1).
//* For optimal efficiency, LWORK >= (NB+1)*N,
//* where NB is the blocksize for CHETRD 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.
//*
//* RWORK (workspace) REAL array, dimension (max(1, 3*N-2))
//*
float *RWORK = (float*) Malloc(MAX(1,3*N-2)*sizeof(float));
//* 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.
//*
int LWORK;
int INFO;
float WORKSZE[2];
LWORK = -1;
cheev_(&JOBZ,&UPLO,&N,a,&LDA,d,WORKSZE,&LWORK,RWORK,&INFO);
LWORK = (int) WORKSZE[0];
float *WORK = (float*) Malloc(2*LWORK*sizeof(float));
cheev_(&JOBZ,&UPLO,&N,a,&LDA,d,WORK,&LWORK,RWORK,&INFO);
Free(WORK);
Free(RWORK);
if (eigenvectors)
memcpy(v,a,2*n*n*sizeof(float));
}
void complexEigenDecompose(int n, float *v, float *d, float *a,
bool eigenvectors, bool balance) {
// SUBROUTINE CGEEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, W, VL,
// $ LDVL, VR, LDVR, ILO, IHI, SCALE, ABNRM, RCONDE,
// $ RCONDV, WORK, LWORK, RWORK, INFO )
//* Purpose
//* =======
//*
//* ZGEEVX computes for an N-by-N complex nonsymmetric matrix A, the
//* eigenvalues and, optionally, the left and/or right eigenvectors.
//*
//* Optionally also, it computes a balancing transformation to improve
//* the conditioning of the eigenvalues and eigenvectors (ILO, IHI,
//* SCALE, and ABNRM), reciprocal condition numbers for the eigenvalues
//* (RCONDE), and reciprocal condition numbers for the right
//* eigenvectors (RCONDV).
//*
//* 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.
//*
//* Balancing a matrix means permuting the rows and columns to make it
//* more nearly upper triangular, and applying a diagonal similarity
//* transformation D * A * D**(-1), where D is a diagonal matrix, to
//* make its rows and columns closer in norm and the condition numbers
//* of its eigenvalues and eigenvectors smaller. The computed
//* reciprocal condition numbers correspond to the balanced matrix.
//* Permuting rows and columns will not change the condition numbers
//* (in exact arithmetic) but diagonal scaling will. For further
//* explanation of balancing, see section 4.10.2 of the LAPACK
//* Users' Guide.
//*
//* Arguments
//* =========
//*
//* BALANC (input) CHARACTER*1
//* Indicates how the input matrix should be diagonally scaled
//* and/or permuted to improve the conditioning of its
//* eigenvalues.
//* = 'N': Do not diagonally scale or permute;
//* = 'P': Perform permutations to make the matrix more nearly
//* upper triangular. Do not diagonally scale;
//* = 'S': Diagonally scale the matrix, ie. replace A by
//* D*A*D**(-1), where D is a diagonal matrix chosen
//* to make the rows and columns of A more equal in
//* norm. Do not permute;
//* = 'B': Both diagonally scale and permute A.
//*
//* Computed reciprocal condition numbers will be for the matrix
//* after balancing and/or permuting. Permuting does not change
//* condition numbers (in exact arithmetic), but balancing does.
char BALANC;
if (balance)
BALANC = 'B';
else
BALANC = 'N';
//* JOBVL (input) CHARACTER*1
//* = 'N': left eigenvectors of A are not computed;
//* = 'V': left eigenvectors of A are computed.
//* If SENSE = 'E' or 'B', JOBVL must = 'V'.
char JOBVL = 'N';
//* JOBVR (input) CHARACTER*1
//* = 'N': right eigenvectors of A are not computed;
//* = 'V': right eigenvectors of A are computed.
//* If SENSE = 'E' or 'B', JOBVR must = 'V'.
char JOBVR;
if (eigenvectors)
JOBVR = 'V';
else
JOBVR = 'N';
//* SENSE (input) CHARACTER*1
//* Determines which reciprocal condition numbers are computed.
//* = 'N': None are computed;
//* = 'E': Computed for eigenvalues only;
//* = 'V': Computed for right eigenvectors only;
//* = 'B': Computed for eigenvalues and right eigenvectors.
//*
//* If SENSE = 'E' or 'B', both left and right eigenvectors
//* must also be computed (JOBVL = 'V' and JOBVR = 'V').
char SENSE = 'N';
//* N (input) INTEGER
//* The order of the matrix A. N >= 0.
int N = n;
//* A (input/output) REAL array, dimension (LDA,N)
//* On entry, the N-by-N matrix A.
//* On exit, A has been overwritten. If JOBVL = 'V' or
//* JOBVR = 'V', A contains the real Schur form of the balanced
//* version of the input matrix A.
float *Ain = a;
//* LDA (input) INTEGER
//* The leading dimension of the array A. LDA >= max(1,N).
int LDA = n;
//* W (output) COMPLEX array, dimension (N)
//* W contains the computed eigenvalues.
float *W = d;
//* VL (output) COMPLEX 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.
//* u(j) = VL(:,j), the j-th column of VL.
//
float *VL = NULL;
//* LDVL (input) INTEGER
//* The leading dimension of the array VL. LDVL >= 1; if
//* JOBVL = 'V', LDVL >= N.
int LDVL = 1;
//* VR (output) COMPLEX 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.
//* v(j) = VR(:,j), the j-th column of VR.
float *VR = v;
//* LDVR (input) INTEGER
//* The leading dimension of the array VR. LDVR >= 1, and if
//* JOBVR = 'V', LDVR >= N.
int LDVR = n;
//* ILO,IHI (output) INTEGER
//* ILO and IHI are integer values determined when A was
//* balanced. The balanced A(i,j) = 0 if I > J and
//* J = 1,...,ILO-1 or I = IHI+1,...,N.
int ILO;
int IHI;
//* SCALE (output) REAL array, dimension (N)
//* Details of the permutations and scaling factors applied
//* when balancing A. If P(j) is the index of the row and column
//* interchanged with row and column j, and D(j) is the scaling
//* factor applied to row and column j, then
//* SCALE(J) = P(J), for J = 1,...,ILO-1
//* = D(J), for J = ILO,...,IHI
//* = P(J) for J = IHI+1,...,N.
//* The order in which the interchanges are made is N to IHI+1,
//* then 1 to ILO-1.
float *SCALE = (float*) Malloc(n*sizeof(float));
//* ABNRM (output) REAL
//* The one-norm of the balanced matrix (the maximum
//* of the sum of absolute values of elements of any column).
float ABNRM;
//* RCONDE (output) REAL array, dimension (N)
//* RCONDE(j) is the reciprocal condition number of the j-th
//* eigenvalue.
float *RCONDE = (float*) Malloc(n*sizeof(float));
//* RCONDV (output) REAL array, dimension (N)
//* RCONDV(j) is the reciprocal condition number of the j-th
//* right eigenvector.
float *RCONDV = (float*) Malloc(n*sizeof(float));
//* WORK (workspace/output) COMPLEX array, dimension (LWORK)
//* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
//* LWORK (input) INTEGER
//* The dimension of the array WORK. If SENSE = 'N' or 'E',
//* LWORK >= max(1,2*N), and if SENSE = 'V' or 'B',
//* LWORK >= N*N+2*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.
int LWORK;
//* RWORK (workspace) REAL array, dimension (2*N)
float *RWORK = (float*) Malloc(2*n*sizeof(float));
//* 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 or condition numbers
//* have been computed; elements 1:ILO-1 and i+1:N of W
//* contain eigenvalues which have converged.
int INFO;
float WORKSZE[2];
LWORK = -1;
cgeevx_( &BALANC, &JOBVL, &JOBVR, &SENSE, &N, Ain, &LDA, W,
VL, &LDVL, VR, &LDVR, &ILO, &IHI, SCALE, &ABNRM,
RCONDE, RCONDV, WORKSZE, &LWORK, RWORK, &INFO );
LWORK = (int) WORKSZE[0];
float *WORK = (float*) Malloc(2*LWORK*sizeof(float));
cgeevx_( &BALANC, &JOBVL, &JOBVR, &SENSE, &N, Ain, &LDA, W,
VL, &LDVL, VR, &LDVR, &ILO, &IHI, SCALE, &ABNRM,
RCONDE, RCONDV, WORK, &LWORK, RWORK, &INFO );
Free(WORK);
Free(SCALE);
Free(RCONDE);
Free(RCONDV);
Free(RWORK);
}
void local_c_div(float *c, float *a, float *b) {
double ratio, den;
double abr, abi, cr;
if( (abr = b[0]) < 0.)
abr = - abr;
if( (abi = b[1]) < 0.)
abi = - abi;
if( abr <= abi ) {
if(abi == 0) {
float af, bf;
af = bf = abr;
if (a[1] != 0 || a[0] != 0)
af = 1.;
c[1] = c[0] = af / bf;
return;
}
ratio = (double)b[0] / b[1] ;
den = b[1] * (1 + ratio*ratio);
cr = (a[0]*ratio + a[1]) / den;
c[1] = (a[1]*ratio - a[0]) / den;
} else {
ratio = (double)b[1] / b[0] ;
den = b[0] * (1 + ratio*ratio);
cr = (a[0] + a[1]*ratio) / den;
c[1] = (a[1] - a[0]*ratio) / den;
}
c[0] = cr;
}
void complexGenEigenDecompose(int n, float *v, float *d, float *a,
float *b, bool eigenvectors) {
// SUBROUTINE CGGEV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA,
// $ VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO )
//* Purpose
//* =======
//*
//* CGGEV computes for a pair of N-by-N complex nonsymmetric matrices
//* (A,B), the generalized eigenvalues, and optionally, the left and/or
//* right generalized eigenvectors.
//*
//* A generalized eigenvalue for a pair of matrices (A,B) is a scalar
//* lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
//* singular. It is usually represented as the pair (alpha,beta), as
//* there is a reasonable interpretation for beta=0, and even for both
//* being zero.
//*
//* The right generalized eigenvector v(j) corresponding to the
//* generalized eigenvalue lambda(j) of (A,B) satisfies
//*
//* A * v(j) = lambda(j) * B * v(j).
//*
//* The left generalized eigenvector u(j) corresponding to the
//* generalized eigenvalues lambda(j) of (A,B) satisfies
//*
//* u(j)**H * A = lambda(j) * u(j)**H * B
//*
//* where u(j)**H is the conjugate-transpose of u(j).
//*
//* Arguments
//* =========
//*
//* JOBVL (input) CHARACTER*1
//* = 'N': do not compute the left generalized eigenvectors;
//* = 'V': compute the left generalized eigenvectors.
char JOBVL = 'N';
//*
//* JOBVR (input) CHARACTER*1
//* = 'N': do not compute the right generalized eigenvectors;
//* = 'V': compute the right generalized eigenvectors.
//*
char JOBVR;
if (eigenvectors)
JOBVR = 'V';
else
JOBVR = 'N';
//* N (input) INTEGER
//* The order of the matrices A, B, VL, and VR. N >= 0.
//*
int N = n;
//* A (input/output) COMPLEX array, dimension (LDA, N)
//* On entry, the matrix A in the pair (A,B).
//* On exit, A has been overwritten.
//*
float *A = a;
//* LDA (input) INTEGER
//* The leading dimension of A. LDA >= max(1,N).
//*
int LDA = n;
//* B (input/output) COMPLEX array, dimension (LDB, N)
//* On entry, the matrix B in the pair (A,B).
//* On exit, B has been overwritten.
//*
float *B = b;
//* LDB (input) INTEGER
//* The leading dimension of B. LDB >= max(1,N).
//*
int LDB = N;
//* ALPHA (output) COMPLEX array, dimension (N)
//* BETA (output) COMPLEX array, dimension (N)
//* On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the
//* generalized eigenvalues.
//*
//* Note: the quotients ALPHA(j)/BETA(j) may easily over- or
//* underflow, and BETA(j) may even be zero. Thus, the user
//* should avoid naively computing the ratio alpha/beta.
//* However, ALPHA will be always less than and usually
//* comparable with norm(A) in magnitude, and BETA always less
//* than and usually comparable with norm(B).
//*
float *ALPHA = (float*) Malloc(2*n*sizeof(float));
float *BETA = (float*) Malloc(2*n*sizeof(float));
//* VL (output) COMPLEX array, dimension (LDVL,N)
//* If JOBVL = 'V', the left generalized eigenvectors u(j) are
//* stored one after another in the columns of VL, in the same
//* order as their eigenvalues.
//* Each eigenvector will be scaled so the largest component
//* will have abs(real part) + abs(imag. part) = 1.
//* Not referenced if JOBVL = 'N'.
//*
float *VL = NULL;
//* LDVL (input) INTEGER
//* The leading dimension of the matrix VL. LDVL >= 1, and
//* if JOBVL = 'V', LDVL >= N.
//*
int LDVL = n;
//* VR (output) COMPLEX array, dimension (LDVR,N)
//* If JOBVR = 'V', the right generalized eigenvectors v(j) are
//* stored one after another in the columns of VR, in the same
//* order as their eigenvalues.
//* Each eigenvector will be scaled so the largest component
//* will have abs(real part) + abs(imag. part) = 1.
//* Not referenced if JOBVR = 'N'.
//*
float *VR = v;
//* LDVR (input) INTEGER
//* The leading dimension of the matrix VR. LDVR >= 1, and
//* if JOBVR = 'V', LDVR >= N.
//*
int LDVR = n;
//* WORK (workspace/output) COMPLEX array, dimension (LWORK)
//* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
//*
//* LWORK (input) INTEGER
//* The dimension of the array WORK. LWORK >= max(1,2*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.
//*
//* RWORK (workspace/output) REAL array, dimension (8*N)
//*
float *RWORK = (float*) Malloc(8*n*sizeof(float));
//* INFO (output) INTEGER
//* = 0: successful exit
//* < 0: if INFO = -i, the i-th argument had an illegal value.
//* =1,...,N:
//* The QZ iteration failed. No eigenvectors have been
//* calculated, but ALPHA(j) and BETA(j) should be
//* correct for j=INFO+1,...,N.
//* > N: =N+1: other then QZ iteration failed in SHGEQZ,
//* =N+2: error return from STGEVC.
//
float WORKSIZE[2];
int LWORK = -1;
int INFO;
cggev_( &JOBVL, &JOBVR, &N, A, &LDA, B, &LDB, ALPHA, BETA,
VL, &LDVL, VR, &LDVR, &WORKSIZE[0], &LWORK, RWORK, &INFO );
LWORK = (int) WORKSIZE[0];
float *WORK = (float*) Malloc(LWORK*2*sizeof(float));
cggev_( &JOBVL, &JOBVR, &N, A, &LDA, B, &LDB, ALPHA, BETA,
VL, &LDVL, VR, &LDVR, WORK, &LWORK, RWORK, &INFO );
int i;
for (i=0;i<n;i++)
local_c_div(d+2*i,ALPHA+2*i,BETA+2*i);
Free(ALPHA);
Free(BETA);
Free(RWORK);
Free(WORK);
}
bool complexGenEigenDecomposeSymmetric(int n, float *v, float *d,
float *a, float *b,
bool eigenvectors) {
// SUBROUTINE CHEGV( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK,
// $ LWORK, RWORK, INFO )
//* Purpose
//* =======
//*
//* CHEGV computes all the eigenvalues, and optionally, the eigenvectors
//* of a complex generalized Hermitian-definite eigenproblem, of the form
//* A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x.
//* Here A and B are assumed to be Hermitian and B is also
//* positive definite.
//*
//* Arguments
//* =========
//*
//* ITYPE (input) INTEGER
//* Specifies the problem type to be solved:
//* = 1: A*x = (lambda)*B*x
//* = 2: A*B*x = (lambda)*x
//* = 3: B*A*x = (lambda)*x
//*
int ITYPE = 1;
//* JOBZ (input) CHARACTER*1
//* = 'N': Compute eigenvalues only;
//* = 'V': Compute eigenvalues and eigenvectors.
//*
char JOBZ;
if (eigenvectors)
JOBZ = 'V';
else
JOBZ = 'N';
//* UPLO (input) CHARACTER*1
//* = 'U': Upper triangles of A and B are stored;
//* = 'L': Lower triangles of A and B are stored.
//*
char UPLO = 'U';
//* N (input) INTEGER
//* The order of the matrices A and B. N >= 0.
//*
int N = n;
//* A (input/output) COMPLEX array, dimension (LDA, N)
//* On entry, the Hermitian 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
//* matrix Z of eigenvectors. The eigenvectors are normalized
//* as follows:
//* if ITYPE = 1 or 2, Z**H*B*Z = I;
//* if ITYPE = 3, Z**H*inv(B)*Z = I.
//* If JOBZ = 'N', then on exit the upper triangle (if UPLO='U')
//* or the lower triangle (if UPLO='L') of A, including the
//* diagonal, is destroyed.
//*
float *A = a;
//* LDA (input) INTEGER
//* The leading dimension of the array A. LDA >= max(1,N).
//*
int LDA = n;
//* B (input/output) COMPLEX array, dimension (LDB, N)
//* On entry, the Hermitian positive definite matrix B.
//* If UPLO = 'U', the leading N-by-N upper triangular part of B
//* contains the upper triangular part of the matrix B.
//* If UPLO = 'L', the leading N-by-N lower triangular part of B
//* contains the lower triangular part of the matrix B.
//*
//* On exit, if INFO <= N, the part of B containing the matrix is
//* overwritten by the triangular factor U or L from the Cholesky
//* factorization B = U**H*U or B = L*L**H.
//*
float *B = b;
//* LDB (input) INTEGER
//* The leading dimension of the array B. LDB >= max(1,N).
//*
int LDB = n;
//* W (output) REAL array, dimension (N)
//* If INFO = 0, the eigenvalues in ascending order.
//*
float *W = d;
//* WORK (workspace/output) COMPLEX array, dimension (LWORK)
//* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
//*
//* LWORK (input) INTEGER
//* The length of the array WORK. LWORK >= max(1,2*N-1).
//* For optimal efficiency, LWORK >= (NB+1)*N,
//* where NB is the blocksize for CHETRD 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.
//*
//* RWORK (workspace) REAL array, dimension (max(1, 3*N-2))
float *RWORK = (float*) Malloc(MAX(1,3*N-2)*sizeof(float));
//*
//* INFO (output) INTEGER
//* = 0: successful exit
//* < 0: if INFO = -i, the i-th argument had an illegal value
//* > 0: CPOTRF or CHEEV returned an error code:
//* <= N: if INFO = i, CHEEV failed to converge;
//* i off-diagonal elements of an intermediate
//* tridiagonal form did not converge to zero;
//* > N: if INFO = N + i, for 1 <= i <= N, then the leading
//* minor of order i of B is not positive definite.
//* The factorization of B could not be completed and
//* no eigenvalues or eigenvectors were computed.
int INFO;
int LWORK;
LWORK = MAX(1,2*N-1);
float *WORK = (float*) Malloc(2*LWORK*sizeof(float));
chegv_(&ITYPE, &JOBZ, &UPLO, &N, A, &LDA, B, &LDB, d, WORK,
&LWORK, RWORK, &INFO );
if (INFO>N) return false;
if (eigenvectors)
memcpy(v,a,2*n*n*sizeof(float));
return true;
}
void dcomplexEigenDecomposeSymmetric(int n, double *v, double *d,
double *a, bool eigenvectors) {
// SUBROUTINE ZHEEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK,
// $ INFO )
//* Purpose
//* =======
//*
//* ZHEEV computes all eigenvalues and, optionally, eigenvectors of a
//* complex Hermitian matrix A.
//*
//* Arguments
//* =========
//*
//* JOBZ (input) CHARACTER*1
//* = 'N': Compute eigenvalues only;
//* = 'V': Compute eigenvalues and eigenvectors.
//*
char JOBZ;
if (eigenvectors)
JOBZ = 'V';
else
JOBZ = 'N';
//* UPLO (input) CHARACTER*1
//* = 'U': Upper triangle of A is stored;
//* = 'L': Lower triangle of A is stored.
//*
char UPLO = 'U';
//* N (input) INTEGER
//* The order of the matrix A. N >= 0.
//*
int N = n;
//* A (input/output) COMPLEX array, dimension (LDA, N)
//* On entry, the Hermitian 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).
//*
int LDA = n;
//* W (output) REAL array, dimension (N)
//* If INFO = 0, the eigenvalues in ascending order.
//*
//* WORK (workspace/output) COMPLEX array, dimension (LWORK)
//* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
//*
//* LWORK (input) INTEGER
//* The length of the array WORK. LWORK >= max(1,2*N-1).
//* For optimal efficiency, LWORK >= (NB+1)*N,
//* where NB is the blocksize for CHETRD 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.
//*
//* RWORK (workspace) REAL array, dimension (max(1, 3*N-2))
//*
double *RWORK = (double*) Malloc(MAX(1,3*N-2)*sizeof(double));
//* 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.
//*
int LWORK;
int INFO;
double WORKSZE[2];
LWORK = -1;
zheev_(&JOBZ,&UPLO,&N,a,&LDA,d,WORKSZE,&LWORK,RWORK,&INFO);
LWORK = (int) WORKSZE[0];
double *WORK = (double*) Malloc(2*LWORK*sizeof(double));
zheev_(&JOBZ,&UPLO,&N,a,&LDA,d,WORK,&LWORK,RWORK,&INFO);
Free(WORK);
Free(RWORK);
if (eigenvectors)
memcpy(v,a,2*n*n*sizeof(double));
}
void dcomplexEigenDecompose(int n, double *v, double *d,
double *a, bool eigenvectors, bool balance) {
// SUBROUTINE ZGEEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, W, VL,
// $ LDVL, VR, LDVR, ILO, IHI, SCALE, ABNRM, RCONDE,
// $ RCONDV, WORK, LWORK, RWORK, INFO )
//* Purpose
//* =======
//*
//* ZGEEVX computes for an N-by-N complex nonsymmetric matrix A, the
//* eigenvalues and, optionally, the left and/or right eigenvectors.
//*
//* Optionally also, it computes a balancing transformation to improve
//* the conditioning of the eigenvalues and eigenvectors (ILO, IHI,
//* SCALE, and ABNRM), reciprocal condition numbers for the eigenvalues
//* (RCONDE), and reciprocal condition numbers for the right
//* eigenvectors (RCONDV).
//*
//* 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.
//*
//* Balancing a matrix means permuting the rows and columns to make it
//* more nearly upper triangular, and applying a diagonal similarity
//* transformation D * A * D**(-1), where D is a diagonal matrix, to
//* make its rows and columns closer in norm and the condition numbers
//* of its eigenvalues and eigenvectors smaller. The computed
//* reciprocal condition numbers correspond to the balanced matrix.
//* Permuting rows and columns will not change the condition numbers
//* (in exact arithmetic) but diagonal scaling will. For further
//* explanation of balancing, see section 4.10.2 of the LAPACK
//* Users' Guide.
//*
//* Arguments
//* =========
//*
//* BALANC (input) CHARACTER*1
//* Indicates how the input matrix should be diagonally scaled
//* and/or permuted to improve the conditioning of its
//* eigenvalues.
//* = 'N': Do not diagonally scale or permute;
//* = 'P': Perform permutations to make the matrix more nearly
//* upper triangular. Do not diagonally scale;
//* = 'S': Diagonally scale the matrix, ie. replace A by
//* D*A*D**(-1), where D is a diagonal matrix chosen
//* to make the rows and columns of A more equal in
//* norm. Do not permute;
//* = 'B': Both diagonally scale and permute A.
//*
//* Computed reciprocal condition numbers will be for the matrix
//* after balancing and/or permuting. Permuting does not change
//* condition numbers (in exact arithmetic), but balancing does.
char BALANC;
if (balance)
BALANC = 'B';
else
BALANC = 'N';
//* JOBVL (input) CHARACTER*1
//* = 'N': left eigenvectors of A are not computed;
//* = 'V': left eigenvectors of A are computed.
//* If SENSE = 'E' or 'B', JOBVL must = 'V'.
char JOBVL = 'N';
//* JOBVR (input) CHARACTER*1
//* = 'N': right eigenvectors of A are not computed;
//* = 'V': right eigenvectors of A are computed.
//* If SENSE = 'E' or 'B', JOBVR must = 'V'.
char JOBVR;
if (eigenvectors)
JOBVR = 'V';
else
JOBVR = 'N';
//* SENSE (input) CHARACTER*1
//* Determines which reciprocal condition numbers are computed.
//* = 'N': None are computed;
//* = 'E': Computed for eigenvalues only;
//* = 'V': Computed for right eigenvectors only;
//* = 'B': Computed for eigenvalues and right eigenvectors.
//*
//* If SENSE = 'E' or 'B', both left and right eigenvectors
//* must also be computed (JOBVL = 'V' and JOBVR = 'V').
char SENSE = 'N';
//* N (input) INTEGER
//* The order of the matrix A. N >= 0.
int N = n;
//* A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
//* On entry, the N-by-N matrix A.
//* On exit, A has been overwritten. If JOBVL = 'V' or
//* JOBVR = 'V', A contains the real Schur form of the balanced
//* version of the input matrix A.
double *Ain = a;
//* LDA (input) INTEGER
//* The leading dimension of the array A. LDA >= max(1,N).
int LDA = n;
//* W (output) COMPLEX*16 array, dimension (N)
//* W contains the computed eigenvalues.
double *W = d;
//* VL (output) COMPLEX*16 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.
//* u(j) = VL(:,j), the j-th column of VL.
//
double *VL = NULL;
//* LDVL (input) INTEGER
//* The leading dimension of the array VL. LDVL >= 1; if
//* JOBVL = 'V', LDVL >= N.
int LDVL = 1;
//* VR (output) COMPLEX*16 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.
//* v(j) = VR(:,j), the j-th column of VR.
double *VR = v;
//* LDVR (input) INTEGER
//* The leading dimension of the array VR. LDVR >= 1, and if
//* JOBVR = 'V', LDVR >= N.
int LDVR = n;
//* ILO,IHI (output) INTEGER
//* ILO and IHI are integer values determined when A was
//* balanced. The balanced A(i,j) = 0 if I > J and
//* J = 1,...,ILO-1 or I = IHI+1,...,N.
int ILO;
int IHI;
//* SCALE (output) DOUBLE PRECISION array, dimension (N)
//* Details of the permutations and scaling factors applied
//* when balancing A. If P(j) is the index of the row and column
//* interchanged with row and column j, and D(j) is the scaling
//* factor applied to row and column j, then
//* SCALE(J) = P(J), for J = 1,...,ILO-1
//* = D(J), for J = ILO,...,IHI
//* = P(J) for J = IHI+1,...,N.
//* The order in which the interchanges are made is N to IHI+1,
//* then 1 to ILO-1.
double *SCALE = (double*) Malloc(n*sizeof(double));
//* ABNRM (output) DOUBLE PRECISION
//* The one-norm of the balanced matrix (the maximum
//* of the sum of absolute values of elements of any column).
double ABNRM;
//* RCONDE (output) DOUBLE PRECISION array, dimension (N)
//* RCONDE(j) is the reciprocal condition number of the j-th
//* eigenvalue.
double *RCONDE = (double*) Malloc(n*sizeof(double));
//* RCONDV (output) DOUBLE PRECISION array, dimension (N)
//* RCONDV(j) is the reciprocal condition number of the j-th
//* right eigenvector.
double *RCONDV = (double*) Malloc(n*sizeof(double));
//* WORK (workspace/output) COMPLEX*16 array, dimension (LWORK)
//* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
//* LWORK (input) INTEGER
//* The dimension of the array WORK. If SENSE = 'N' or 'E',
//* LWORK >= max(1,2*N), and if SENSE = 'V' or 'B',
//* LWORK >= N*N+2*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.
int LWORK;
//* RWORK (workspace) DOUBLE PRECISION array, dimension (2*N)
double *RWORK = (double*) Malloc(2*n*sizeof(double));
//* 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 or condition numbers
//* have been computed; elements 1:ILO-1 and i+1:N of W
//* contain eigenvalues which have converged.
int INFO;
double WORKSZE[2];
LWORK = -1;
zgeevx_( &BALANC, &JOBVL, &JOBVR, &SENSE, &N, Ain, &LDA, W,
VL, &LDVL, VR, &LDVR, &ILO, &IHI, SCALE, &ABNRM,
RCONDE, RCONDV, WORKSZE, &LWORK, RWORK, &INFO );
LWORK = (int) WORKSZE[0];
double *WORK = (double*) Malloc(2*LWORK*sizeof(double));
zgeevx_( &BALANC, &JOBVL, &JOBVR, &SENSE, &N, Ain, &LDA, W,
VL, &LDVL, VR, &LDVR, &ILO, &IHI, SCALE, &ABNRM,
RCONDE, RCONDV, WORK, &LWORK, RWORK, &INFO );
Free(WORK);
Free(SCALE);
Free(RCONDE);
Free(RCONDV);
Free(RWORK);
}
void local_z_div(double *c, double *a, double *b) {
double ratio, den;
double abr, abi, cr;
if( (abr = b[0]) < 0.)
abr = - abr;
if( (abi = b[1]) < 0.)
abi = - abi;
if( abr <= abi ) {
if(abi == 0) {
if (a[1] != 0 || a[0] != 0)
abi = 1.;
c[1] = c[0] = abi / abr;
return;
}
ratio = b[0] / b[1] ;
den = b[1] * (1 + ratio*ratio);
cr = (a[0]*ratio + a[1]) / den;
c[1] = (a[1]*ratio - a[0]) / den;
} else {
ratio = b[1] / b[0] ;
den = b[0] * (1 + ratio*ratio);
cr = (a[0] + a[1]*ratio) / den;
c[1] = (a[1] - a[0]*ratio) / den;
}
c[0] = cr;
}
bool dcomplexGenEigenDecomposeSymmetric(int n, double *v, double *d,
double *a, double *b,
bool eigenvectors) {
int ITYPE = 1;
char JOBZ;
if (eigenvectors)
JOBZ = 'V';
else
JOBZ = 'N';
char UPLO = 'U';
int N = n;
double *A = a;
int LDA = n;
double *B = b;
int LDB = n;
double *W = d;
double *RWORK = (double*) Malloc(MAX(1,3*N-2)*sizeof(double));
int INFO;
int LWORK;
double WORKSZE[2];
LWORK = MAX(1,2*N-1);
double *WORK = (double*) Malloc(LWORK*sizeof(double)*2);
zhegv_(&ITYPE, &JOBZ, &UPLO, &N, A, &LDA, B, &LDB, d, WORK,
&LWORK, RWORK, &INFO );
Free(WORK);
Free(RWORK);
if (INFO>N) return false;
if (eigenvectors)
memcpy(v,a,2*n*n*sizeof(double));
return true;
}
void dcomplexGenEigenDecompose(int n, double *v, double *d,
double *a, double *b,
bool eigenvectors) {
char JOBVL = 'N';
char JOBVR;
if (eigenvectors)
JOBVR = 'V';
else
JOBVR = 'N';
int N = n;
double *A = a;
int LDA = n;
double *B = b;
int LDB = N;
double *ALPHA = (double*) Malloc(2*n*sizeof(double));
double *BETA = (double*) Malloc(2*n*sizeof(double));
double *VL = NULL;
int LDVL = n;
double *VR = v;
int LDVR = n;
double *RWORK = (double*) Malloc(8*n*sizeof(double));
double WORKSIZE[2];
int LWORK = -1;
int INFO;
zggev_( &JOBVL, &JOBVR, &N, A, &LDA, B, &LDB, ALPHA, BETA,
VL, &LDVL, VR, &LDVR, &WORKSIZE[0], &LWORK, RWORK, &INFO );
LWORK = (int) WORKSIZE[0];
double *WORK = (double*) Malloc(LWORK*2*sizeof(double));
zggev_( &JOBVL, &JOBVR, &N, A, &LDA, B, &LDB, ALPHA, BETA,
VL, &LDVL, VR, &LDVR, WORK, &LWORK, RWORK, &INFO );
int i;
for (i=0;i<n;i++)
local_z_div(d+2*i,ALPHA+2*i,BETA+2*i);
Free(ALPHA);
Free(BETA);
Free(RWORK);
Free(WORK);
}
syntax highlighted by Code2HTML, v. 0.9.1