|
Packit |
67cb25 |
/* eigen/test.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 1996, 1997, 1998, 1999, 2000, 2006, 2007 Gerard Jungman, Patrick Alken, Brian Gough
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
67cb25 |
* it under the terms of the GNU General Public License as published by
|
|
Packit |
67cb25 |
* the Free Software Foundation; either version 3 of the License, or (at
|
|
Packit |
67cb25 |
* your option) any later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is distributed in the hope that it will be useful, but
|
|
Packit |
67cb25 |
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
67cb25 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
67cb25 |
* General Public License for more details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* You should have received a copy of the GNU General Public License
|
|
Packit |
67cb25 |
* along with this program; if not, write to the Free Software
|
|
Packit |
67cb25 |
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Author: G. Jungman
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <stdlib.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_test.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_blas.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_ieee_utils.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_complex_math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_eigen.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_linalg.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_rng.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_randist.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_errno.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_sort.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_sort_vector.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_matrix.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_vector.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/******************************************
|
|
Packit |
67cb25 |
* common test code *
|
|
Packit |
67cb25 |
******************************************/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double
|
|
Packit |
67cb25 |
chop_subnormals (double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* Chop any subnormal values */
|
|
Packit |
67cb25 |
return fabs(x) < GSL_DBL_MIN ? 0 : x;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
create_random_symm_matrix(gsl_matrix *m, gsl_rng *r, int lower, int upper)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t i, j;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < m->size1; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for (j = i; j < m->size2; ++j)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double x = gsl_rng_uniform(r) * (upper - lower) + lower;
|
|
Packit |
67cb25 |
gsl_matrix_set(m, i, j, x);
|
|
Packit |
67cb25 |
gsl_matrix_set(m, j, i, x);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
} /* create_random_symm_matrix() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
create_random_herm_matrix(gsl_matrix_complex *m, gsl_rng *r, int lower,
|
|
Packit |
67cb25 |
int upper)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t i, j;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < m->size1; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for (j = i; j < m->size2; ++j)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_complex z;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
GSL_REAL(z) = gsl_rng_uniform(r) * (upper - lower) + lower;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (i == j)
|
|
Packit |
67cb25 |
GSL_IMAG(z) = 0.0;
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
GSL_IMAG(z) = gsl_rng_uniform(r) * (upper - lower) + lower;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_complex_set(m, i, j, z);
|
|
Packit |
67cb25 |
gsl_matrix_complex_set(m, j, i, gsl_complex_conjugate(z));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
} /* create_random_herm_matrix() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* with r \in (0,1) if m_{ij} = r^{|i - j|} then m is positive definite */
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
create_random_posdef_matrix(gsl_matrix *m, gsl_rng *r)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t i, j;
|
|
Packit |
67cb25 |
double x = gsl_rng_uniform(r);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < m->size1; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for (j = i; j < m->size2; ++j)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double a = pow(x, (double) (j - i));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_set(m, i, j, a);
|
|
Packit |
67cb25 |
gsl_matrix_set(m, j, i, a);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
} /* create_random_posdef_matrix() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
create_random_complex_posdef_matrix(gsl_matrix_complex *m, gsl_rng *r,
|
|
Packit |
67cb25 |
gsl_vector_complex *work)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t N = m->size1;
|
|
Packit |
67cb25 |
size_t i, j;
|
|
Packit |
67cb25 |
double x, y;
|
|
Packit |
67cb25 |
gsl_complex z;
|
|
Packit |
67cb25 |
gsl_complex tau;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
GSL_SET_IMAG(&z, 0.0);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* make a positive diagonal matrix */
|
|
Packit |
67cb25 |
gsl_matrix_complex_set_zero(m);
|
|
Packit |
67cb25 |
for (i = 0; i < N; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
x = gsl_rng_uniform(r);
|
|
Packit |
67cb25 |
GSL_SET_REAL(&z, x);
|
|
Packit |
67cb25 |
gsl_matrix_complex_set(m, i, i, z);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* now generate random householder reflections and form P D P^H */
|
|
Packit |
67cb25 |
for (i = 0; i < N; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* form complex vector */
|
|
Packit |
67cb25 |
for (j = 0; j < N; ++j)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
x = 2.0 * gsl_rng_uniform(r) - 1.0;
|
|
Packit |
67cb25 |
y = 2.0 * gsl_rng_uniform(r) - 1.0;
|
|
Packit |
67cb25 |
GSL_SET_COMPLEX(&z, x, y);
|
|
Packit |
67cb25 |
gsl_vector_complex_set(work, j, z);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
tau = gsl_linalg_complex_householder_transform(work);
|
|
Packit |
67cb25 |
gsl_linalg_complex_householder_hm(tau, work, m);
|
|
Packit |
67cb25 |
gsl_linalg_complex_householder_mh(gsl_complex_conjugate(tau), work, m);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
} /* create_random_complex_posdef_matrix() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
create_random_nonsymm_matrix(gsl_matrix *m, gsl_rng *r, int lower,
|
|
Packit |
67cb25 |
int upper)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t i, j;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < m->size1; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for (j = 0; j < m->size2; ++j)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_set(m,
|
|
Packit |
67cb25 |
i,
|
|
Packit |
67cb25 |
j,
|
|
Packit |
67cb25 |
gsl_rng_uniform(r) * (upper - lower) + lower);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
} /* create_random_nonsymm_matrix() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* test if A Z = Q S */
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
test_eigen_schur(const gsl_matrix * A, const gsl_matrix * S,
|
|
Packit |
67cb25 |
const gsl_matrix * Q, const gsl_matrix * Z,
|
|
Packit |
67cb25 |
size_t count, const char * desc,
|
|
Packit |
67cb25 |
const char * desc2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t N = A->size1;
|
|
Packit |
67cb25 |
size_t i, j;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix * T1 = gsl_matrix_alloc(N, N);
|
|
Packit |
67cb25 |
gsl_matrix * T2 = gsl_matrix_alloc(N, N);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute T1 = A Z */
|
|
Packit |
67cb25 |
gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, A, Z, 0.0, T1);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute T2 = Q S */
|
|
Packit |
67cb25 |
gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, Q, S, 0.0, T2);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < N; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for (j = 0; j < N; ++j)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double x = gsl_matrix_get(T1, i, j);
|
|
Packit |
67cb25 |
double y = gsl_matrix_get(T2, i, j);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_test_abs(x, y, 1.0e8 * GSL_DBL_EPSILON,
|
|
Packit |
67cb25 |
"%s(N=%u,cnt=%u), %s, schur(%d,%d)", desc, N, count, desc2, i, j);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_free (T1);
|
|
Packit |
67cb25 |
gsl_matrix_free (T2);
|
|
Packit |
67cb25 |
} /* test_eigen_schur() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* test if A is orthogonal */
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
test_eigen_orthog(gsl_matrix *A, size_t count, const char *desc,
|
|
Packit |
67cb25 |
const char *desc2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t N = A->size1;
|
|
Packit |
67cb25 |
gsl_matrix *M = gsl_matrix_alloc(A->size1, A->size2);
|
|
Packit |
67cb25 |
size_t i, j;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1.0, A, A, 0.0, M);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < A->size1; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for (j = 0; j < A->size2; ++j)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double val;
|
|
Packit |
67cb25 |
double mij = gsl_matrix_get(M, i, j);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (i == j)
|
|
Packit |
67cb25 |
val = 1.0;
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
val = 0.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_test_abs(mij, val, 1.0e8 * GSL_DBL_EPSILON,
|
|
Packit |
67cb25 |
"%s(N=%u,cnt=%u), %s, orthog(%d,%d)", desc, N, count, desc2, i, j);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_free(M);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return 1;
|
|
Packit |
67cb25 |
} /* test_eigen_orthog() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
test_eigenvalues_real (const gsl_vector *eval, const gsl_vector * eval2,
|
|
Packit |
67cb25 |
const char * desc, const char * desc2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t N = eval->size;
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double emax = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* check eigenvalues */
|
|
Packit |
67cb25 |
for (i = 0; i < N; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double ei = gsl_vector_get (eval, i);
|
|
Packit |
67cb25 |
if (fabs(ei) > emax) emax = fabs(ei);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < N; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double ei = gsl_vector_get (eval, i);
|
|
Packit |
67cb25 |
double e2i = gsl_vector_get (eval2, i);
|
|
Packit |
67cb25 |
e2i = chop_subnormals(e2i);
|
|
Packit |
67cb25 |
gsl_test_abs(ei, e2i, emax * 1e8 * GSL_DBL_EPSILON,
|
|
Packit |
67cb25 |
"%s, direct eigenvalue(%d), %s",
|
|
Packit |
67cb25 |
desc, i, desc2);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
test_eigenvalues_complex (const gsl_vector_complex * eval,
|
|
Packit |
67cb25 |
const gsl_vector_complex * eval2,
|
|
Packit |
67cb25 |
const char * desc, const char * desc2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t N = eval->size;
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < N; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_complex ei = gsl_vector_complex_get (eval, i);
|
|
Packit |
67cb25 |
gsl_complex e2i = gsl_vector_complex_get (eval2, i);
|
|
Packit |
67cb25 |
gsl_test_rel(GSL_REAL(ei), GSL_REAL(e2i), 10*N*GSL_DBL_EPSILON,
|
|
Packit |
67cb25 |
"%s, direct eigenvalue(%d) real, %s",
|
|
Packit |
67cb25 |
desc, i, desc2);
|
|
Packit |
67cb25 |
gsl_test_rel(GSL_IMAG(ei), GSL_IMAG(e2i), 10*N*GSL_DBL_EPSILON,
|
|
Packit |
67cb25 |
"%s, direct eigenvalue(%d) imag, %s",
|
|
Packit |
67cb25 |
desc, i, desc2);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/******************************************
|
|
Packit |
67cb25 |
* symm test code *
|
|
Packit |
67cb25 |
******************************************/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
test_eigen_symm_results (const gsl_matrix * A,
|
|
Packit |
67cb25 |
const gsl_vector * eval,
|
|
Packit |
67cb25 |
const gsl_matrix * evec,
|
|
Packit |
67cb25 |
size_t count,
|
|
Packit |
67cb25 |
const char * desc,
|
|
Packit |
67cb25 |
const char * desc2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t N = A->size1;
|
|
Packit |
67cb25 |
size_t i, j;
|
|
Packit |
67cb25 |
double emax = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector * x = gsl_vector_alloc(N);
|
|
Packit |
67cb25 |
gsl_vector * y = gsl_vector_alloc(N);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* check eigenvalues */
|
|
Packit |
67cb25 |
for (i = 0; i < N; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double ei = gsl_vector_get (eval, i);
|
|
Packit |
67cb25 |
if (fabs(ei) > emax) emax = fabs(ei);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < N; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double ei = gsl_vector_get (eval, i);
|
|
Packit |
67cb25 |
gsl_vector_const_view vi = gsl_matrix_const_column(evec, i);
|
|
Packit |
67cb25 |
gsl_vector_memcpy(x, &vi.vector);
|
|
Packit |
67cb25 |
/* compute y = A x (should = lambda v) */
|
|
Packit |
67cb25 |
gsl_blas_dgemv (CblasNoTrans, 1.0, A, x, 0.0, y);
|
|
Packit |
67cb25 |
for (j = 0; j < N; j++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double xj = gsl_vector_get (x, j);
|
|
Packit |
67cb25 |
double yj = gsl_vector_get (y, j);
|
|
Packit |
67cb25 |
double eixj = chop_subnormals(ei * xj);
|
|
Packit |
67cb25 |
gsl_test_abs(yj, eixj, emax * 1e8 * GSL_DBL_EPSILON,
|
|
Packit |
67cb25 |
"%s, eigenvalue(%d,%d), %s", desc, i, j, desc2);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* check eigenvectors are orthonormal */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < N; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_const_view vi = gsl_matrix_const_column(evec, i);
|
|
Packit |
67cb25 |
double nrm_v = gsl_blas_dnrm2(&vi.vector);
|
|
Packit |
67cb25 |
gsl_test_rel (nrm_v, 1.0, N * GSL_DBL_EPSILON, "%s, normalized(%d), %s",
|
|
Packit |
67cb25 |
desc, i, desc2);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < N; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_const_view vi = gsl_matrix_const_column(evec, i);
|
|
Packit |
67cb25 |
for (j = i + 1; j < N; j++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_const_view vj = gsl_matrix_const_column(evec, j);
|
|
Packit |
67cb25 |
double vivj;
|
|
Packit |
67cb25 |
gsl_blas_ddot (&vi.vector, &vj.vector, &vivj);
|
|
Packit |
67cb25 |
gsl_test_abs (vivj, 0.0, N * GSL_DBL_EPSILON,
|
|
Packit |
67cb25 |
"%s, orthogonal(%d,%d), %s", desc, i, j, desc2);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_free(x);
|
|
Packit |
67cb25 |
gsl_vector_free(y);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
test_eigen_symm_matrix(const gsl_matrix * m, size_t count,
|
|
Packit |
67cb25 |
const char * desc)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t N = m->size1;
|
|
Packit |
67cb25 |
gsl_matrix * A = gsl_matrix_alloc(N, N);
|
|
Packit |
67cb25 |
gsl_vector * eval = gsl_vector_alloc(N);
|
|
Packit |
67cb25 |
gsl_vector * evalv = gsl_vector_alloc(N);
|
|
Packit |
67cb25 |
gsl_vector * x = gsl_vector_alloc(N);
|
|
Packit |
67cb25 |
gsl_vector * y = gsl_vector_alloc(N);
|
|
Packit |
67cb25 |
gsl_matrix * evec = gsl_matrix_alloc(N, N);
|
|
Packit |
67cb25 |
gsl_eigen_symm_workspace * w = gsl_eigen_symm_alloc(N);
|
|
Packit |
67cb25 |
gsl_eigen_symmv_workspace * wv = gsl_eigen_symmv_alloc(N);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_memcpy(A, m);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_eigen_symmv(A, evalv, evec, wv);
|
|
Packit |
67cb25 |
test_eigen_symm_results(m, evalv, evec, count, desc, "unsorted");
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_memcpy(A, m);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_eigen_symm(A, eval, w);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* sort eval and evalv */
|
|
Packit |
67cb25 |
gsl_vector_memcpy(x, eval);
|
|
Packit |
67cb25 |
gsl_vector_memcpy(y, evalv);
|
|
Packit |
67cb25 |
gsl_sort_vector(x);
|
|
Packit |
67cb25 |
gsl_sort_vector(y);
|
|
Packit |
67cb25 |
test_eigenvalues_real(y, x, desc, "unsorted");
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_eigen_symmv_sort(evalv, evec, GSL_EIGEN_SORT_VAL_ASC);
|
|
Packit |
67cb25 |
test_eigen_symm_results(m, evalv, evec, count, desc, "val/asc");
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_eigen_symmv_sort(evalv, evec, GSL_EIGEN_SORT_VAL_DESC);
|
|
Packit |
67cb25 |
test_eigen_symm_results(m, evalv, evec, count, desc, "val/desc");
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_eigen_symmv_sort(evalv, evec, GSL_EIGEN_SORT_ABS_ASC);
|
|
Packit |
67cb25 |
test_eigen_symm_results(m, evalv, evec, count, desc, "abs/asc");
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_eigen_symmv_sort(evalv, evec, GSL_EIGEN_SORT_ABS_DESC);
|
|
Packit |
67cb25 |
test_eigen_symm_results(m, evalv, evec, count, desc, "abs/desc");
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_free(A);
|
|
Packit |
67cb25 |
gsl_vector_free(eval);
|
|
Packit |
67cb25 |
gsl_vector_free(evalv);
|
|
Packit |
67cb25 |
gsl_vector_free(x);
|
|
Packit |
67cb25 |
gsl_vector_free(y);
|
|
Packit |
67cb25 |
gsl_matrix_free(evec);
|
|
Packit |
67cb25 |
gsl_eigen_symm_free(w);
|
|
Packit |
67cb25 |
gsl_eigen_symmv_free(wv);
|
|
Packit |
67cb25 |
} /* test_eigen_symm_matrix() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
test_eigen_symm(void)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t N_max = 20;
|
|
Packit |
67cb25 |
size_t n, i;
|
|
Packit |
67cb25 |
gsl_rng *r = gsl_rng_alloc(gsl_rng_default);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (n = 1; n <= N_max; ++n)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix * A = gsl_matrix_alloc(n, n);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < 5; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
create_random_symm_matrix(A, r, -10, 10);
|
|
Packit |
67cb25 |
test_eigen_symm_matrix(A, i, "symm random");
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_free(A);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_rng_free(r);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double dat1[] = { 0, 0, -1, 0,
|
|
Packit |
67cb25 |
0, 1, 0, 1,
|
|
Packit |
67cb25 |
-1, 0, 0, 0,
|
|
Packit |
67cb25 |
0, 1, 0, 0 };
|
|
Packit |
67cb25 |
double dat2[] = { 1, 0, 0, 0,
|
|
Packit |
67cb25 |
0, 2, 0, 0,
|
|
Packit |
67cb25 |
0, 0, 3, 0,
|
|
Packit |
67cb25 |
0, 0, 0, 4 };
|
|
Packit |
67cb25 |
gsl_matrix_view m;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
m = gsl_matrix_view_array (dat1, 4, 4);
|
|
Packit |
67cb25 |
test_eigen_symm_matrix(&m.matrix, 0, "symm(4)");
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
m = gsl_matrix_view_array (dat2, 4, 4);
|
|
Packit |
67cb25 |
test_eigen_symm_matrix(&m.matrix, 0, "symm(4) diag");
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double dat[27*27] = {
|
|
Packit |
67cb25 |
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
|
|
Packit |
67cb25 |
0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,
|
|
Packit |
67cb25 |
0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,
|
|
Packit |
67cb25 |
0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
|
|
Packit |
67cb25 |
0,0,0,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,
|
|
Packit |
67cb25 |
0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
|
|
Packit |
67cb25 |
0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
|
|
Packit |
67cb25 |
0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
|
|
Packit |
67cb25 |
0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
|
|
Packit |
67cb25 |
0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
|
|
Packit |
67cb25 |
0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
|
|
Packit |
67cb25 |
0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
|
|
Packit |
67cb25 |
0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
|
|
Packit |
67cb25 |
0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
|
|
Packit |
67cb25 |
0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
|
|
Packit |
67cb25 |
0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
|
|
Packit |
67cb25 |
0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
|
|
Packit |
67cb25 |
0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
|
|
Packit |
67cb25 |
0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
|
|
Packit |
67cb25 |
0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
|
|
Packit |
67cb25 |
0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
|
|
Packit |
67cb25 |
0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
|
|
Packit |
67cb25 |
0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
|
|
Packit |
67cb25 |
0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
|
|
Packit |
67cb25 |
0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
|
|
Packit |
67cb25 |
0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
|
|
Packit |
67cb25 |
0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_view m;
|
|
Packit |
67cb25 |
m = gsl_matrix_view_array (dat, 27, 27);
|
|
Packit |
67cb25 |
test_eigen_symm_matrix(&m.matrix, 0, "symm(27)");
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
} /* test_eigen_symm() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/******************************************
|
|
Packit |
67cb25 |
* herm test code *
|
|
Packit |
67cb25 |
******************************************/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
test_eigen_herm_results (const gsl_matrix_complex * A,
|
|
Packit |
67cb25 |
const gsl_vector * eval,
|
|
Packit |
67cb25 |
const gsl_matrix_complex * evec,
|
|
Packit |
67cb25 |
size_t count,
|
|
Packit |
67cb25 |
const char * desc,
|
|
Packit |
67cb25 |
const char * desc2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t N = A->size1;
|
|
Packit |
67cb25 |
size_t i, j;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_complex * x = gsl_vector_complex_alloc(N);
|
|
Packit |
67cb25 |
gsl_vector_complex * y = gsl_vector_complex_alloc(N);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* check eigenvalues */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < N; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double ei = gsl_vector_get (eval, i);
|
|
Packit |
67cb25 |
gsl_vector_complex_const_view vi =
|
|
Packit |
67cb25 |
gsl_matrix_complex_const_column(evec, i);
|
|
Packit |
67cb25 |
gsl_vector_complex_memcpy(x, &vi.vector);
|
|
Packit |
67cb25 |
/* compute y = m x (should = lambda v) */
|
|
Packit |
67cb25 |
gsl_blas_zgemv (CblasNoTrans, GSL_COMPLEX_ONE, A, x,
|
|
Packit |
67cb25 |
GSL_COMPLEX_ZERO, y);
|
|
Packit |
67cb25 |
for (j = 0; j < N; j++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_complex xj = gsl_vector_complex_get (x, j);
|
|
Packit |
67cb25 |
gsl_complex yj = gsl_vector_complex_get (y, j);
|
|
Packit |
67cb25 |
gsl_test_rel(GSL_REAL(yj), ei * GSL_REAL(xj), 1e8*GSL_DBL_EPSILON,
|
|
Packit |
67cb25 |
"%s, eigenvalue(%d,%d), real, %s", desc, i, j, desc2);
|
|
Packit |
67cb25 |
gsl_test_rel(GSL_IMAG(yj), ei * GSL_IMAG(xj), 1e8*GSL_DBL_EPSILON,
|
|
Packit |
67cb25 |
"%s, eigenvalue(%d,%d), imag, %s", desc, i, j, desc2);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* check eigenvectors are orthonormal */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < N; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_complex_const_view vi = gsl_matrix_complex_const_column(evec, i);
|
|
Packit |
67cb25 |
double nrm_v = gsl_blas_dznrm2(&vi.vector);
|
|
Packit |
67cb25 |
gsl_test_rel (nrm_v, 1.0, N * GSL_DBL_EPSILON, "%s, normalized(%d), %s",
|
|
Packit |
67cb25 |
desc, i, desc2);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < N; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_complex_const_view vi = gsl_matrix_complex_const_column(evec, i);
|
|
Packit |
67cb25 |
for (j = i + 1; j < N; j++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_complex_const_view vj
|
|
Packit |
67cb25 |
= gsl_matrix_complex_const_column(evec, j);
|
|
Packit |
67cb25 |
gsl_complex vivj;
|
|
Packit |
67cb25 |
gsl_blas_zdotc (&vi.vector, &vj.vector, &vivj);
|
|
Packit |
67cb25 |
gsl_test_abs (gsl_complex_abs(vivj), 0.0, 10.0 * N * GSL_DBL_EPSILON,
|
|
Packit |
67cb25 |
"%s, orthogonal(%d,%d), %s", desc, i, j, desc2);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_complex_free(x);
|
|
Packit |
67cb25 |
gsl_vector_complex_free(y);
|
|
Packit |
67cb25 |
} /* test_eigen_herm_results() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
test_eigen_herm_matrix(const gsl_matrix_complex * m, size_t count,
|
|
Packit |
67cb25 |
const char * desc)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t N = m->size1;
|
|
Packit |
67cb25 |
gsl_matrix_complex * A = gsl_matrix_complex_alloc(N, N);
|
|
Packit |
67cb25 |
gsl_vector * eval = gsl_vector_alloc(N);
|
|
Packit |
67cb25 |
gsl_vector * evalv = gsl_vector_alloc(N);
|
|
Packit |
67cb25 |
gsl_vector * x = gsl_vector_alloc(N);
|
|
Packit |
67cb25 |
gsl_vector * y = gsl_vector_alloc(N);
|
|
Packit |
67cb25 |
gsl_matrix_complex * evec = gsl_matrix_complex_alloc(N, N);
|
|
Packit |
67cb25 |
gsl_eigen_herm_workspace * w = gsl_eigen_herm_alloc(N);
|
|
Packit |
67cb25 |
gsl_eigen_hermv_workspace * wv = gsl_eigen_hermv_alloc(N);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_complex_memcpy(A, m);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_eigen_hermv(A, evalv, evec, wv);
|
|
Packit |
67cb25 |
test_eigen_herm_results(m, evalv, evec, count, desc, "unsorted");
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_complex_memcpy(A, m);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_eigen_herm(A, eval, w);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* sort eval and evalv */
|
|
Packit |
67cb25 |
gsl_vector_memcpy(x, eval);
|
|
Packit |
67cb25 |
gsl_vector_memcpy(y, evalv);
|
|
Packit |
67cb25 |
gsl_sort_vector(x);
|
|
Packit |
67cb25 |
gsl_sort_vector(y);
|
|
Packit |
67cb25 |
test_eigenvalues_real(y, x, desc, "unsorted");
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_eigen_hermv_sort(evalv, evec, GSL_EIGEN_SORT_VAL_ASC);
|
|
Packit |
67cb25 |
test_eigen_herm_results(m, evalv, evec, count, desc, "val/asc");
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_eigen_hermv_sort(evalv, evec, GSL_EIGEN_SORT_VAL_DESC);
|
|
Packit |
67cb25 |
test_eigen_herm_results(m, evalv, evec, count, desc, "val/desc");
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_eigen_hermv_sort(evalv, evec, GSL_EIGEN_SORT_ABS_ASC);
|
|
Packit |
67cb25 |
test_eigen_herm_results(m, evalv, evec, count, desc, "abs/asc");
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_eigen_hermv_sort(evalv, evec, GSL_EIGEN_SORT_ABS_DESC);
|
|
Packit |
67cb25 |
test_eigen_herm_results(m, evalv, evec, count, desc, "abs/desc");
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_complex_free(A);
|
|
Packit |
67cb25 |
gsl_vector_free(eval);
|
|
Packit |
67cb25 |
gsl_vector_free(evalv);
|
|
Packit |
67cb25 |
gsl_vector_free(x);
|
|
Packit |
67cb25 |
gsl_vector_free(y);
|
|
Packit |
67cb25 |
gsl_matrix_complex_free(evec);
|
|
Packit |
67cb25 |
gsl_eigen_herm_free(w);
|
|
Packit |
67cb25 |
gsl_eigen_hermv_free(wv);
|
|
Packit |
67cb25 |
} /* test_eigen_herm_matrix() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
test_eigen_herm(void)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t N_max = 20;
|
|
Packit |
67cb25 |
size_t n, i;
|
|
Packit |
67cb25 |
gsl_rng *r = gsl_rng_alloc(gsl_rng_default);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (n = 1; n <= N_max; ++n)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_complex * A = gsl_matrix_complex_alloc(n, n);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < 5; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
create_random_herm_matrix(A, r, -10, 10);
|
|
Packit |
67cb25 |
test_eigen_herm_matrix(A, i, "herm random");
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_complex_free(A);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_rng_free(r);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double dat1[] = { 0,0, 0,0, -1,0, 0,0,
|
|
Packit |
67cb25 |
0,0, 1,0, 0,0, 1,0,
|
|
Packit |
67cb25 |
-1,0, 0,0, 0,0, 0,0,
|
|
Packit |
67cb25 |
0,0, 1,0, 0,0, 0,0 };
|
|
Packit |
67cb25 |
double dat2[] = { 1,0, 0,0, 0,0, 0,0,
|
|
Packit |
67cb25 |
0,0, 2,0, 0,0, 0,0,
|
|
Packit |
67cb25 |
0,0, 0,0, 3,0, 0,0,
|
|
Packit |
67cb25 |
0,0, 0,0, 0,0, 4,0 };
|
|
Packit |
67cb25 |
gsl_matrix_complex_view m;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
m = gsl_matrix_complex_view_array (dat1, 4, 4);
|
|
Packit |
67cb25 |
test_eigen_herm_matrix(&m.matrix, 0, "herm(4)");
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
m = gsl_matrix_complex_view_array (dat2, 4, 4);
|
|
Packit |
67cb25 |
test_eigen_herm_matrix(&m.matrix, 1, "herm(4) diag");
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
} /* test_eigen_herm() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/******************************************
|
|
Packit |
67cb25 |
* nonsymm test code *
|
|
Packit |
67cb25 |
******************************************/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
test_eigen_nonsymm_results (const gsl_matrix * m,
|
|
Packit |
67cb25 |
const gsl_vector_complex * eval,
|
|
Packit |
67cb25 |
const gsl_matrix_complex * evec,
|
|
Packit |
67cb25 |
size_t count,
|
|
Packit |
67cb25 |
const char * desc,
|
|
Packit |
67cb25 |
const char * desc2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t i,j;
|
|
Packit |
67cb25 |
size_t N = m->size1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_complex * x = gsl_vector_complex_alloc(N);
|
|
Packit |
67cb25 |
gsl_vector_complex * y = gsl_vector_complex_alloc(N);
|
|
Packit |
67cb25 |
gsl_matrix_complex * A = gsl_matrix_complex_alloc(N, N);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* we need a complex matrix for the blas routines, so copy m into A */
|
|
Packit |
67cb25 |
for (i = 0; i < N; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for (j = 0; j < N; ++j)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_complex z;
|
|
Packit |
67cb25 |
GSL_SET_COMPLEX(&z, gsl_matrix_get(m, i, j), 0.0);
|
|
Packit |
67cb25 |
gsl_matrix_complex_set(A, i, j, z);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < N; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_complex ei = gsl_vector_complex_get (eval, i);
|
|
Packit |
67cb25 |
gsl_vector_complex_const_view vi = gsl_matrix_complex_const_column(evec, i);
|
|
Packit |
67cb25 |
double norm = gsl_blas_dznrm2(&vi.vector);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* check that eigenvector is normalized */
|
|
Packit |
67cb25 |
gsl_test_rel(norm, 1.0, N * GSL_DBL_EPSILON,
|
|
Packit |
67cb25 |
"nonsymm(N=%u,cnt=%u), %s, normalized(%d), %s", N, count, desc, i, desc2);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_complex_memcpy(x, &vi.vector);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute y = m x (should = lambda v) */
|
|
Packit |
67cb25 |
gsl_blas_zgemv (CblasNoTrans, GSL_COMPLEX_ONE, A, x,
|
|
Packit |
67cb25 |
GSL_COMPLEX_ZERO, y);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute x = lambda v */
|
|
Packit |
67cb25 |
gsl_blas_zscal(ei, x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* now test if y = x */
|
|
Packit |
67cb25 |
for (j = 0; j < N; j++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_complex xj = gsl_vector_complex_get (x, j);
|
|
Packit |
67cb25 |
gsl_complex yj = gsl_vector_complex_get (y, j);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* use abs here in case the values are close to 0 */
|
|
Packit |
67cb25 |
gsl_test_abs(GSL_REAL(yj), GSL_REAL(xj), 1e8*GSL_DBL_EPSILON,
|
|
Packit |
67cb25 |
"nonsymm(N=%u,cnt=%u), %s, eigenvalue(%d,%d), real, %s", N, count, desc, i, j, desc2);
|
|
Packit |
67cb25 |
gsl_test_abs(GSL_IMAG(yj), GSL_IMAG(xj), 1e8*GSL_DBL_EPSILON,
|
|
Packit |
67cb25 |
"nonsymm(N=%u,cnt=%u), %s, eigenvalue(%d,%d), imag, %s", N, count, desc, i, j, desc2);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_complex_free(A);
|
|
Packit |
67cb25 |
gsl_vector_complex_free(x);
|
|
Packit |
67cb25 |
gsl_vector_complex_free(y);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
test_eigen_nonsymm_matrix(const gsl_matrix * m, size_t count,
|
|
Packit |
67cb25 |
const char * desc,
|
|
Packit |
67cb25 |
gsl_eigen_nonsymmv_workspace *w)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t N = m->size1;
|
|
Packit |
67cb25 |
gsl_matrix * A = gsl_matrix_alloc(N, N);
|
|
Packit |
67cb25 |
gsl_matrix * Z = gsl_matrix_alloc(N, N);
|
|
Packit |
67cb25 |
gsl_matrix_complex * evec = gsl_matrix_complex_alloc(N, N);
|
|
Packit |
67cb25 |
gsl_vector_complex * eval = gsl_vector_complex_alloc(N);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* calculate eigenvalues and eigenvectors - it is sufficient to
|
|
Packit |
67cb25 |
* test gsl_eigen_nonsymmv() since that function calls
|
|
Packit |
67cb25 |
* gsl_eigen_nonsymm() for the eigenvalues
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
gsl_matrix_memcpy(A, m);
|
|
Packit |
67cb25 |
gsl_eigen_nonsymmv(A, eval, evec, w);
|
|
Packit |
67cb25 |
test_eigen_nonsymm_results (m, eval, evec, count, desc, "unsorted");
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* test sort routines */
|
|
Packit |
67cb25 |
gsl_eigen_nonsymmv_sort (eval, evec, GSL_EIGEN_SORT_ABS_ASC);
|
|
Packit |
67cb25 |
test_eigen_nonsymm_results (m, eval, evec, count, desc, "abs/asc");
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_eigen_nonsymmv_sort (eval, evec, GSL_EIGEN_SORT_ABS_DESC);
|
|
Packit |
67cb25 |
test_eigen_nonsymm_results (m, eval, evec, count, desc, "abs/desc");
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* test Schur vectors */
|
|
Packit |
67cb25 |
gsl_matrix_memcpy(A, m);
|
|
Packit |
67cb25 |
gsl_eigen_nonsymmv_Z(A, eval, evec, Z, w);
|
|
Packit |
67cb25 |
gsl_linalg_hessenberg_set_zero(A);
|
|
Packit |
67cb25 |
test_eigen_schur(m, A, Z, Z, count, "nonsymm", desc);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* test if Z is an orthogonal matrix */
|
|
Packit |
67cb25 |
if (w->nonsymm_workspace_p->do_balance == 0)
|
|
Packit |
67cb25 |
test_eigen_orthog(Z, count, "nonsymm", desc);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_free(A);
|
|
Packit |
67cb25 |
gsl_matrix_free(Z);
|
|
Packit |
67cb25 |
gsl_matrix_complex_free(evec);
|
|
Packit |
67cb25 |
gsl_vector_complex_free(eval);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
test_eigen_nonsymm(void)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t N_max = 20;
|
|
Packit |
67cb25 |
size_t n, i;
|
|
Packit |
67cb25 |
gsl_rng *r = gsl_rng_alloc(gsl_rng_default);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (n = 1; n <= N_max; ++n)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix * m = gsl_matrix_alloc(n, n);
|
|
Packit |
67cb25 |
gsl_eigen_nonsymmv_workspace * w = gsl_eigen_nonsymmv_alloc(n);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < 5; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
create_random_nonsymm_matrix(m, r, -10, 10);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_eigen_nonsymmv_params(0, w);
|
|
Packit |
67cb25 |
test_eigen_nonsymm_matrix(m, i, "random, unbalanced", w);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_eigen_nonsymmv_params(1, w);
|
|
Packit |
67cb25 |
test_eigen_nonsymm_matrix(m, i, "random, balanced", w);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_free(m);
|
|
Packit |
67cb25 |
gsl_eigen_nonsymmv_free(w);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_rng_free(r);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double dat1[] = { 0, 1, 1, 1,
|
|
Packit |
67cb25 |
1, 1, 1, 1,
|
|
Packit |
67cb25 |
0, 0, 0, 0,
|
|
Packit |
67cb25 |
0, 0, 0, 0 };
|
|
Packit |
67cb25 |
double dat2[] = { 1, 1, 0, 1,
|
|
Packit |
67cb25 |
1, 1, 1, 1,
|
|
Packit |
67cb25 |
1, 1, 1, 1,
|
|
Packit |
67cb25 |
0, 1, 0, 0 };
|
|
Packit |
67cb25 |
gsl_matrix_view v;
|
|
Packit |
67cb25 |
gsl_eigen_nonsymmv_workspace * w = gsl_eigen_nonsymmv_alloc(4);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
v = gsl_matrix_view_array (dat1, 4, 4);
|
|
Packit |
67cb25 |
test_eigen_nonsymm_matrix(&v.matrix, 0, "integer", w);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
v = gsl_matrix_view_array (dat2, 4, 4);
|
|
Packit |
67cb25 |
test_eigen_nonsymm_matrix(&v.matrix, 1, "integer", w);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_eigen_nonsymmv_free(w);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
} /* test_eigen_nonsymm() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/******************************************
|
|
Packit |
67cb25 |
* gensymm test code *
|
|
Packit |
67cb25 |
******************************************/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
test_eigen_gensymm_results (const gsl_matrix * A,
|
|
Packit |
67cb25 |
const gsl_matrix * B,
|
|
Packit |
67cb25 |
const gsl_vector * eval,
|
|
Packit |
67cb25 |
const gsl_matrix * evec,
|
|
Packit |
67cb25 |
size_t count,
|
|
Packit |
67cb25 |
const char * desc,
|
|
Packit |
67cb25 |
const char * desc2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t N = A->size1;
|
|
Packit |
67cb25 |
size_t i, j;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector * x = gsl_vector_alloc(N);
|
|
Packit |
67cb25 |
gsl_vector * y = gsl_vector_alloc(N);
|
|
Packit |
67cb25 |
gsl_vector * z = gsl_vector_alloc(N);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* check A v = lambda B v */
|
|
Packit |
67cb25 |
for (i = 0; i < N; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double ei = gsl_vector_get (eval, i);
|
|
Packit |
67cb25 |
gsl_vector_const_view vi = gsl_matrix_const_column(evec, i);
|
|
Packit |
67cb25 |
double norm = gsl_blas_dnrm2(&vi.vector);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* check that eigenvector is normalized */
|
|
Packit |
67cb25 |
gsl_test_rel(norm, 1.0, N * GSL_DBL_EPSILON,
|
|
Packit |
67cb25 |
"gensymm(N=%u,cnt=%u), %s, normalized(%d), %s", N, count,
|
|
Packit |
67cb25 |
desc, i, desc2);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_memcpy(z, &vi.vector);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute y = A z */
|
|
Packit |
67cb25 |
gsl_blas_dgemv (CblasNoTrans, 1.0, A, z, 0.0, y);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute x = B z */
|
|
Packit |
67cb25 |
gsl_blas_dgemv (CblasNoTrans, 1.0, B, z, 0.0, x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute x = lambda B z */
|
|
Packit |
67cb25 |
gsl_blas_dscal(ei, x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* now test if y = x */
|
|
Packit |
67cb25 |
for (j = 0; j < N; j++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double xj = gsl_vector_get (x, j);
|
|
Packit |
67cb25 |
double yj = gsl_vector_get (y, j);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_test_rel(yj, xj, 1e9 * GSL_DBL_EPSILON,
|
|
Packit |
67cb25 |
"gensymm(N=%u,cnt=%u), %s, eigenvalue(%d,%d), real, %s", N, count, desc, i, j, desc2);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_free(x);
|
|
Packit |
67cb25 |
gsl_vector_free(y);
|
|
Packit |
67cb25 |
gsl_vector_free(z);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
test_eigen_gensymm(void)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t N_max = 20;
|
|
Packit |
67cb25 |
size_t n, i;
|
|
Packit |
67cb25 |
gsl_rng *r = gsl_rng_alloc(gsl_rng_default);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (n = 1; n <= N_max; ++n)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix * A = gsl_matrix_alloc(n, n);
|
|
Packit |
67cb25 |
gsl_matrix * B = gsl_matrix_alloc(n, n);
|
|
Packit |
67cb25 |
gsl_matrix * ma = gsl_matrix_alloc(n, n);
|
|
Packit |
67cb25 |
gsl_matrix * mb = gsl_matrix_alloc(n, n);
|
|
Packit |
67cb25 |
gsl_vector * eval = gsl_vector_alloc(n);
|
|
Packit |
67cb25 |
gsl_vector * evalv = gsl_vector_alloc(n);
|
|
Packit |
67cb25 |
gsl_vector * x = gsl_vector_alloc(n);
|
|
Packit |
67cb25 |
gsl_vector * y = gsl_vector_alloc(n);
|
|
Packit |
67cb25 |
gsl_matrix * evec = gsl_matrix_alloc(n, n);
|
|
Packit |
67cb25 |
gsl_eigen_gensymm_workspace * w = gsl_eigen_gensymm_alloc(n);
|
|
Packit |
67cb25 |
gsl_eigen_gensymmv_workspace * wv = gsl_eigen_gensymmv_alloc(n);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < 5; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
create_random_symm_matrix(A, r, -10, 10);
|
|
Packit |
67cb25 |
create_random_posdef_matrix(B, r);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_memcpy(ma, A);
|
|
Packit |
67cb25 |
gsl_matrix_memcpy(mb, B);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_eigen_gensymmv(ma, mb, evalv, evec, wv);
|
|
Packit |
67cb25 |
test_eigen_gensymm_results(A, B, evalv, evec, i, "random", "unsorted");
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_memcpy(ma, A);
|
|
Packit |
67cb25 |
gsl_matrix_memcpy(mb, B);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_eigen_gensymm(ma, mb, eval, w);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* eval and evalv have to be sorted? not sure why */
|
|
Packit |
67cb25 |
gsl_vector_memcpy(x, eval);
|
|
Packit |
67cb25 |
gsl_vector_memcpy(y, evalv);
|
|
Packit |
67cb25 |
gsl_sort_vector(x);
|
|
Packit |
67cb25 |
gsl_sort_vector(y);
|
|
Packit |
67cb25 |
test_eigenvalues_real(y, x, "gensymm, random", "unsorted");
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_eigen_gensymmv_sort(evalv, evec, GSL_EIGEN_SORT_VAL_ASC);
|
|
Packit |
67cb25 |
test_eigen_gensymm_results(A, B, evalv, evec, i, "random", "val/asc");
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_eigen_gensymmv_sort(evalv, evec, GSL_EIGEN_SORT_VAL_DESC);
|
|
Packit |
67cb25 |
test_eigen_gensymm_results(A, B, evalv, evec, i, "random", "val/desc");
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_eigen_gensymmv_sort(evalv, evec, GSL_EIGEN_SORT_ABS_ASC);
|
|
Packit |
67cb25 |
test_eigen_gensymm_results(A, B, evalv, evec, i, "random", "abs/asc");
|
|
Packit |
67cb25 |
gsl_eigen_gensymmv_sort(evalv, evec, GSL_EIGEN_SORT_ABS_DESC);
|
|
Packit |
67cb25 |
test_eigen_gensymm_results(A, B, evalv, evec, i, "random", "abs/desc");
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_free(A);
|
|
Packit |
67cb25 |
gsl_matrix_free(B);
|
|
Packit |
67cb25 |
gsl_matrix_free(ma);
|
|
Packit |
67cb25 |
gsl_matrix_free(mb);
|
|
Packit |
67cb25 |
gsl_vector_free(eval);
|
|
Packit |
67cb25 |
gsl_vector_free(evalv);
|
|
Packit |
67cb25 |
gsl_vector_free(x);
|
|
Packit |
67cb25 |
gsl_vector_free(y);
|
|
Packit |
67cb25 |
gsl_matrix_free(evec);
|
|
Packit |
67cb25 |
gsl_eigen_gensymm_free(w);
|
|
Packit |
67cb25 |
gsl_eigen_gensymmv_free(wv);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_rng_free(r);
|
|
Packit |
67cb25 |
} /* test_eigen_gensymm() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/******************************************
|
|
Packit |
67cb25 |
* genherm test code *
|
|
Packit |
67cb25 |
******************************************/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
test_eigen_genherm_results (const gsl_matrix_complex * A,
|
|
Packit |
67cb25 |
const gsl_matrix_complex * B,
|
|
Packit |
67cb25 |
const gsl_vector * eval,
|
|
Packit |
67cb25 |
const gsl_matrix_complex * evec,
|
|
Packit |
67cb25 |
size_t count,
|
|
Packit |
67cb25 |
const char * desc,
|
|
Packit |
67cb25 |
const char * desc2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t N = A->size1;
|
|
Packit |
67cb25 |
size_t i, j;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_complex * x = gsl_vector_complex_alloc(N);
|
|
Packit |
67cb25 |
gsl_vector_complex * y = gsl_vector_complex_alloc(N);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* check A v = lambda B v */
|
|
Packit |
67cb25 |
for (i = 0; i < N; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double ei = gsl_vector_get (eval, i);
|
|
Packit |
67cb25 |
gsl_vector_complex_const_view vi =
|
|
Packit |
67cb25 |
gsl_matrix_complex_const_column(evec, i);
|
|
Packit |
67cb25 |
double norm = gsl_blas_dznrm2(&vi.vector);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* check that eigenvector is normalized */
|
|
Packit |
67cb25 |
gsl_test_rel(norm, 1.0, N * GSL_DBL_EPSILON,
|
|
Packit |
67cb25 |
"genherm(N=%u,cnt=%u), %s, normalized(%d), %s", N, count,
|
|
Packit |
67cb25 |
desc, i, desc2);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute y = A z */
|
|
Packit |
67cb25 |
gsl_blas_zgemv (CblasNoTrans, GSL_COMPLEX_ONE, A, &vi.vector, GSL_COMPLEX_ZERO, y);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute x = B z */
|
|
Packit |
67cb25 |
gsl_blas_zgemv (CblasNoTrans, GSL_COMPLEX_ONE, B, &vi.vector, GSL_COMPLEX_ZERO, x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute x = lambda B z */
|
|
Packit |
67cb25 |
gsl_blas_zdscal(ei, x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* now test if y = x */
|
|
Packit |
67cb25 |
for (j = 0; j < N; j++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_complex xj = gsl_vector_complex_get (x, j);
|
|
Packit |
67cb25 |
gsl_complex yj = gsl_vector_complex_get (y, j);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_test_rel(GSL_REAL(yj), GSL_REAL(xj), 1e9 * GSL_DBL_EPSILON,
|
|
Packit |
67cb25 |
"genherm(N=%u,cnt=%u), %s, eigenvalue(%d,%d), real, %s", N, count, desc, i, j, desc2);
|
|
Packit |
67cb25 |
gsl_test_abs(GSL_IMAG(yj), GSL_IMAG(xj), 1e9 * GSL_DBL_EPSILON,
|
|
Packit |
67cb25 |
"genherm(N=%u,cnt=%u), %s, eigenvalue(%d,%d), imag, %s", N, count, desc, i, j, desc2);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_complex_free(x);
|
|
Packit |
67cb25 |
gsl_vector_complex_free(y);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
test_eigen_genherm(void)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t N_max = 20;
|
|
Packit |
67cb25 |
size_t n, i;
|
|
Packit |
67cb25 |
gsl_rng *r = gsl_rng_alloc(gsl_rng_default);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (n = 1; n <= N_max; ++n)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_complex * A = gsl_matrix_complex_alloc(n, n);
|
|
Packit |
67cb25 |
gsl_matrix_complex * B = gsl_matrix_complex_alloc(n, n);
|
|
Packit |
67cb25 |
gsl_matrix_complex * ma = gsl_matrix_complex_alloc(n, n);
|
|
Packit |
67cb25 |
gsl_matrix_complex * mb = gsl_matrix_complex_alloc(n, n);
|
|
Packit |
67cb25 |
gsl_vector * eval = gsl_vector_alloc(n);
|
|
Packit |
67cb25 |
gsl_vector * evalv = gsl_vector_alloc(n);
|
|
Packit |
67cb25 |
gsl_vector * x = gsl_vector_alloc(n);
|
|
Packit |
67cb25 |
gsl_vector * y = gsl_vector_alloc(n);
|
|
Packit |
67cb25 |
gsl_vector_complex * work = gsl_vector_complex_alloc(n);
|
|
Packit |
67cb25 |
gsl_matrix_complex * evec = gsl_matrix_complex_alloc(n, n);
|
|
Packit |
67cb25 |
gsl_eigen_genherm_workspace * w = gsl_eigen_genherm_alloc(n);
|
|
Packit |
67cb25 |
gsl_eigen_genhermv_workspace * wv = gsl_eigen_genhermv_alloc(n);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < 5; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
create_random_herm_matrix(A, r, -10, 10);
|
|
Packit |
67cb25 |
create_random_complex_posdef_matrix(B, r, work);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_complex_memcpy(ma, A);
|
|
Packit |
67cb25 |
gsl_matrix_complex_memcpy(mb, B);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_eigen_genhermv(ma, mb, evalv, evec, wv);
|
|
Packit |
67cb25 |
test_eigen_genherm_results(A, B, evalv, evec, i, "random", "unsorted");
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_complex_memcpy(ma, A);
|
|
Packit |
67cb25 |
gsl_matrix_complex_memcpy(mb, B);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_eigen_genherm(ma, mb, eval, w);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* eval and evalv have to be sorted? not sure why */
|
|
Packit |
67cb25 |
gsl_vector_memcpy(x, eval);
|
|
Packit |
67cb25 |
gsl_vector_memcpy(y, evalv);
|
|
Packit |
67cb25 |
gsl_sort_vector(x);
|
|
Packit |
67cb25 |
gsl_sort_vector(y);
|
|
Packit |
67cb25 |
test_eigenvalues_real(y, x, "genherm, random", "unsorted");
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_eigen_genhermv_sort(evalv, evec, GSL_EIGEN_SORT_VAL_ASC);
|
|
Packit |
67cb25 |
test_eigen_genherm_results(A, B, evalv, evec, i, "random", "val/asc");
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_eigen_genhermv_sort(evalv, evec, GSL_EIGEN_SORT_VAL_DESC);
|
|
Packit |
67cb25 |
test_eigen_genherm_results(A, B, evalv, evec, i, "random", "val/desc");
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_eigen_genhermv_sort(evalv, evec, GSL_EIGEN_SORT_ABS_ASC);
|
|
Packit |
67cb25 |
test_eigen_genherm_results(A, B, evalv, evec, i, "random", "abs/asc");
|
|
Packit |
67cb25 |
gsl_eigen_genhermv_sort(evalv, evec, GSL_EIGEN_SORT_ABS_DESC);
|
|
Packit |
67cb25 |
test_eigen_genherm_results(A, B, evalv, evec, i, "random", "abs/desc");
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_complex_free(A);
|
|
Packit |
67cb25 |
gsl_matrix_complex_free(B);
|
|
Packit |
67cb25 |
gsl_matrix_complex_free(ma);
|
|
Packit |
67cb25 |
gsl_matrix_complex_free(mb);
|
|
Packit |
67cb25 |
gsl_vector_free(eval);
|
|
Packit |
67cb25 |
gsl_vector_free(evalv);
|
|
Packit |
67cb25 |
gsl_vector_free(x);
|
|
Packit |
67cb25 |
gsl_vector_free(y);
|
|
Packit |
67cb25 |
gsl_vector_complex_free(work);
|
|
Packit |
67cb25 |
gsl_matrix_complex_free(evec);
|
|
Packit |
67cb25 |
gsl_eigen_genherm_free(w);
|
|
Packit |
67cb25 |
gsl_eigen_genhermv_free(wv);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_rng_free(r);
|
|
Packit |
67cb25 |
} /* test_eigen_genherm() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/******************************************
|
|
Packit |
67cb25 |
* gen test code *
|
|
Packit |
67cb25 |
******************************************/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
typedef struct
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix *A;
|
|
Packit |
67cb25 |
gsl_matrix *B;
|
|
Packit |
67cb25 |
gsl_vector_complex *alpha;
|
|
Packit |
67cb25 |
gsl_vector *beta;
|
|
Packit |
67cb25 |
gsl_vector_complex *alphav;
|
|
Packit |
67cb25 |
gsl_vector *betav;
|
|
Packit |
67cb25 |
gsl_vector_complex *eval;
|
|
Packit |
67cb25 |
gsl_vector_complex *evalv;
|
|
Packit |
67cb25 |
gsl_vector *x;
|
|
Packit |
67cb25 |
gsl_vector *y;
|
|
Packit |
67cb25 |
gsl_matrix *Q;
|
|
Packit |
67cb25 |
gsl_matrix *Z;
|
|
Packit |
67cb25 |
gsl_matrix_complex *evec;
|
|
Packit |
67cb25 |
gsl_eigen_gen_workspace *gen_p;
|
|
Packit |
67cb25 |
gsl_eigen_genv_workspace *genv_p;
|
|
Packit |
67cb25 |
} test_eigen_gen_workspace;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
test_eigen_gen_workspace *
|
|
Packit |
67cb25 |
test_eigen_gen_alloc(const size_t n)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
test_eigen_gen_workspace *w;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
w = (test_eigen_gen_workspace *) calloc(1, sizeof(test_eigen_gen_workspace));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
w->A = gsl_matrix_alloc(n, n);
|
|
Packit |
67cb25 |
w->B = gsl_matrix_alloc(n, n);
|
|
Packit |
67cb25 |
w->alpha = gsl_vector_complex_alloc(n);
|
|
Packit |
67cb25 |
w->beta = gsl_vector_alloc(n);
|
|
Packit |
67cb25 |
w->alphav = gsl_vector_complex_alloc(n);
|
|
Packit |
67cb25 |
w->betav = gsl_vector_alloc(n);
|
|
Packit |
67cb25 |
w->eval = gsl_vector_complex_alloc(n);
|
|
Packit |
67cb25 |
w->evalv = gsl_vector_complex_alloc(n);
|
|
Packit |
67cb25 |
w->x = gsl_vector_alloc(n);
|
|
Packit |
67cb25 |
w->y = gsl_vector_alloc(n);
|
|
Packit |
67cb25 |
w->Q = gsl_matrix_alloc(n, n);
|
|
Packit |
67cb25 |
w->Z = gsl_matrix_alloc(n, n);
|
|
Packit |
67cb25 |
w->evec = gsl_matrix_complex_alloc(n, n);
|
|
Packit |
67cb25 |
w->gen_p = gsl_eigen_gen_alloc(n);
|
|
Packit |
67cb25 |
w->genv_p = gsl_eigen_genv_alloc(n);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return (w);
|
|
Packit |
67cb25 |
} /* test_eigen_gen_alloc() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
test_eigen_gen_free(test_eigen_gen_workspace *w)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_free(w->A);
|
|
Packit |
67cb25 |
gsl_matrix_free(w->B);
|
|
Packit |
67cb25 |
gsl_vector_complex_free(w->alpha);
|
|
Packit |
67cb25 |
gsl_vector_free(w->beta);
|
|
Packit |
67cb25 |
gsl_vector_complex_free(w->alphav);
|
|
Packit |
67cb25 |
gsl_vector_free(w->betav);
|
|
Packit |
67cb25 |
gsl_vector_complex_free(w->eval);
|
|
Packit |
67cb25 |
gsl_vector_complex_free(w->evalv);
|
|
Packit |
67cb25 |
gsl_vector_free(w->x);
|
|
Packit |
67cb25 |
gsl_vector_free(w->y);
|
|
Packit |
67cb25 |
gsl_matrix_free(w->Q);
|
|
Packit |
67cb25 |
gsl_matrix_free(w->Z);
|
|
Packit |
67cb25 |
gsl_matrix_complex_free(w->evec);
|
|
Packit |
67cb25 |
gsl_eigen_gen_free(w->gen_p);
|
|
Packit |
67cb25 |
gsl_eigen_genv_free(w->genv_p);
|
|
Packit |
67cb25 |
free(w);
|
|
Packit |
67cb25 |
} /* test_eigen_gen_free() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
test_eigen_gen_results (const gsl_matrix * A, const gsl_matrix * B,
|
|
Packit |
67cb25 |
const gsl_vector_complex * alpha,
|
|
Packit |
67cb25 |
const gsl_vector * beta,
|
|
Packit |
67cb25 |
const gsl_matrix_complex * evec,
|
|
Packit |
67cb25 |
size_t count, const char * desc,
|
|
Packit |
67cb25 |
const char * desc2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t N = A->size1;
|
|
Packit |
67cb25 |
size_t i, j;
|
|
Packit |
67cb25 |
gsl_matrix_complex *ma, *mb;
|
|
Packit |
67cb25 |
gsl_vector_complex *x, *y;
|
|
Packit |
67cb25 |
gsl_complex z_one, z_zero;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
ma = gsl_matrix_complex_alloc(N, N);
|
|
Packit |
67cb25 |
mb = gsl_matrix_complex_alloc(N, N);
|
|
Packit |
67cb25 |
y = gsl_vector_complex_alloc(N);
|
|
Packit |
67cb25 |
x = gsl_vector_complex_alloc(N);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* ma <- A, mb <- B */
|
|
Packit |
67cb25 |
for (i = 0; i < N; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for (j = 0; j < N; ++j)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_complex z;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
GSL_SET_COMPLEX(&z, gsl_matrix_get(A, i, j), 0.0);
|
|
Packit |
67cb25 |
gsl_matrix_complex_set(ma, i, j, z);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
GSL_SET_COMPLEX(&z, gsl_matrix_get(B, i, j), 0.0);
|
|
Packit |
67cb25 |
gsl_matrix_complex_set(mb, i, j, z);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
GSL_SET_COMPLEX(&z_one, 1.0, 0.0);
|
|
Packit |
67cb25 |
GSL_SET_COMPLEX(&z_zero, 0.0, 0.0);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* check eigenvalues */
|
|
Packit |
67cb25 |
for (i = 0; i < N; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_complex_const_view vi =
|
|
Packit |
67cb25 |
gsl_matrix_complex_const_column(evec, i);
|
|
Packit |
67cb25 |
gsl_complex ai = gsl_vector_complex_get(alpha, i);
|
|
Packit |
67cb25 |
double bi = gsl_vector_get(beta, i);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute x = alpha * B * v */
|
|
Packit |
67cb25 |
gsl_blas_zgemv(CblasNoTrans, z_one, mb, &vi.vector, z_zero, x);
|
|
Packit |
67cb25 |
gsl_blas_zscal(ai, x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute y = beta * A v */
|
|
Packit |
67cb25 |
gsl_blas_zgemv(CblasNoTrans, z_one, ma, &vi.vector, z_zero, y);
|
|
Packit |
67cb25 |
gsl_blas_zdscal(bi, y);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* now test if y = x */
|
|
Packit |
67cb25 |
for (j = 0; j < N; ++j)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_complex xj = gsl_vector_complex_get(x, j);
|
|
Packit |
67cb25 |
gsl_complex yj = gsl_vector_complex_get(y, j);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_test_abs(GSL_REAL(yj), GSL_REAL(xj), 1e8*GSL_DBL_EPSILON,
|
|
Packit |
67cb25 |
"gen(N=%u,cnt=%u), %s, eigenvalue(%d,%d), real, %s",
|
|
Packit |
67cb25 |
N, count, desc, i, j, desc2);
|
|
Packit |
67cb25 |
gsl_test_abs(GSL_IMAG(yj), GSL_IMAG(xj), 1e8*GSL_DBL_EPSILON,
|
|
Packit |
67cb25 |
"gen(N=%u,cnt=%u), %s, eigenvalue(%d,%d), real, %s",
|
|
Packit |
67cb25 |
N, count, desc, i, j, desc2);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_complex_free(ma);
|
|
Packit |
67cb25 |
gsl_matrix_complex_free(mb);
|
|
Packit |
67cb25 |
gsl_vector_complex_free(y);
|
|
Packit |
67cb25 |
gsl_vector_complex_free(x);
|
|
Packit |
67cb25 |
} /* test_eigen_gen_results() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
test_eigen_gen_pencil(const gsl_matrix * A, const gsl_matrix * B,
|
|
Packit |
67cb25 |
size_t count, const char * desc, int test_schur,
|
|
Packit |
67cb25 |
test_eigen_gen_workspace *w)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t N = A->size1;
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_memcpy(w->A, A);
|
|
Packit |
67cb25 |
gsl_matrix_memcpy(w->B, B);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (test_schur)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_eigen_genv_QZ(w->A, w->B, w->alphav, w->betav, w->evec, w->Q, w->Z, w->genv_p);
|
|
Packit |
67cb25 |
test_eigen_schur(A, w->A, w->Q, w->Z, count, "genv/A", desc);
|
|
Packit |
67cb25 |
test_eigen_schur(B, w->B, w->Q, w->Z, count, "genv/B", desc);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
gsl_eigen_genv(w->A, w->B, w->alphav, w->betav, w->evec, w->genv_p);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
test_eigen_gen_results(A, B, w->alphav, w->betav, w->evec, count, desc, "unsorted");
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_memcpy(w->A, A);
|
|
Packit |
67cb25 |
gsl_matrix_memcpy(w->B, B);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (test_schur)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_eigen_gen_params(1, 1, 0, w->gen_p);
|
|
Packit |
67cb25 |
gsl_eigen_gen_QZ(w->A, w->B, w->alpha, w->beta, w->Q, w->Z, w->gen_p);
|
|
Packit |
67cb25 |
test_eigen_schur(A, w->A, w->Q, w->Z, count, "gen/A", desc);
|
|
Packit |
67cb25 |
test_eigen_schur(B, w->B, w->Q, w->Z, count, "gen/B", desc);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_eigen_gen_params(0, 0, 0, w->gen_p);
|
|
Packit |
67cb25 |
gsl_eigen_gen(w->A, w->B, w->alpha, w->beta, w->gen_p);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute eval = alpha / beta values */
|
|
Packit |
67cb25 |
for (i = 0; i < N; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_complex z, ai;
|
|
Packit |
67cb25 |
double bi;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
ai = gsl_vector_complex_get(w->alpha, i);
|
|
Packit |
67cb25 |
bi = gsl_vector_get(w->beta, i);
|
|
Packit |
67cb25 |
GSL_SET_COMPLEX(&z, GSL_REAL(ai) / bi, GSL_IMAG(ai) / bi);
|
|
Packit |
67cb25 |
gsl_vector_complex_set(w->eval, i, z);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
ai = gsl_vector_complex_get(w->alphav, i);
|
|
Packit |
67cb25 |
bi = gsl_vector_get(w->betav, i);
|
|
Packit |
67cb25 |
GSL_SET_COMPLEX(&z, GSL_REAL(ai) / bi, GSL_IMAG(ai) / bi);
|
|
Packit |
67cb25 |
gsl_vector_complex_set(w->evalv, i, z);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* sort eval and evalv and test them */
|
|
Packit |
67cb25 |
gsl_eigen_nonsymmv_sort(w->eval, NULL, GSL_EIGEN_SORT_VAL_ASC);
|
|
Packit |
67cb25 |
gsl_eigen_nonsymmv_sort(w->evalv, NULL, GSL_EIGEN_SORT_VAL_ASC);
|
|
Packit |
67cb25 |
test_eigenvalues_complex(w->evalv, w->eval, "gen", desc);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_eigen_genv_sort(w->alphav, w->betav, w->evec, GSL_EIGEN_SORT_ABS_ASC);
|
|
Packit |
67cb25 |
test_eigen_gen_results(A, B, w->alphav, w->betav, w->evec, count, desc, "abs/asc");
|
|
Packit |
67cb25 |
gsl_eigen_genv_sort(w->alphav, w->betav, w->evec, GSL_EIGEN_SORT_ABS_DESC);
|
|
Packit |
67cb25 |
test_eigen_gen_results(A, B, w->alphav, w->betav, w->evec, count, desc, "abs/desc");
|
|
Packit |
67cb25 |
} /* test_eigen_gen_pencil() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
test_eigen_gen(void)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t N_max = 20;
|
|
Packit |
67cb25 |
size_t n, i;
|
|
Packit |
67cb25 |
gsl_rng *r = gsl_rng_alloc(gsl_rng_default);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (n = 1; n <= N_max; ++n)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix * A = gsl_matrix_alloc(n, n);
|
|
Packit |
67cb25 |
gsl_matrix * B = gsl_matrix_alloc(n, n);
|
|
Packit |
67cb25 |
test_eigen_gen_workspace * w = test_eigen_gen_alloc(n);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < 5; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
create_random_nonsymm_matrix(A, r, -10, 10);
|
|
Packit |
67cb25 |
create_random_nonsymm_matrix(B, r, -10, 10);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
test_eigen_gen_pencil(A, B, i, "random", 0, w);
|
|
Packit |
67cb25 |
test_eigen_gen_pencil(A, B, i, "random", 1, w);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_matrix_free(A);
|
|
Packit |
67cb25 |
gsl_matrix_free(B);
|
|
Packit |
67cb25 |
test_eigen_gen_free(w);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_rng_free(r);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* this system will test the exceptional shift code */
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double datA[] = { 1, 1, 0,
|
|
Packit |
67cb25 |
0, 0, -1,
|
|
Packit |
67cb25 |
1, 0, 0 };
|
|
Packit |
67cb25 |
double datB[] = { -1, 0, -1,
|
|
Packit |
67cb25 |
0, -1, 0,
|
|
Packit |
67cb25 |
0, 0, -1 };
|
|
Packit |
67cb25 |
gsl_matrix_view va = gsl_matrix_view_array (datA, 3, 3);
|
|
Packit |
67cb25 |
gsl_matrix_view vb = gsl_matrix_view_array (datB, 3, 3);
|
|
Packit |
67cb25 |
test_eigen_gen_workspace * w = test_eigen_gen_alloc(3);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
test_eigen_gen_pencil(&va.matrix, &vb.matrix, 0, "integer", 0, w);
|
|
Packit |
67cb25 |
test_eigen_gen_pencil(&va.matrix, &vb.matrix, 0, "integer", 1, w);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
test_eigen_gen_free(w);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
} /* test_eigen_gen() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
main()
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_ieee_env_setup ();
|
|
Packit |
67cb25 |
gsl_rng_env_setup ();
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
test_eigen_symm();
|
|
Packit |
67cb25 |
test_eigen_herm();
|
|
Packit |
67cb25 |
test_eigen_nonsymm();
|
|
Packit |
67cb25 |
test_eigen_gensymm();
|
|
Packit |
67cb25 |
test_eigen_genherm();
|
|
Packit |
67cb25 |
test_eigen_gen();
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
exit (gsl_test_summary());
|
|
Packit |
67cb25 |
}
|