|
Packit |
67cb25 |
/* linalg/balance.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 2006 Patrick Alken
|
|
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 |
/* Balance a general matrix by scaling the rows and columns, so the
|
|
Packit |
67cb25 |
* new row and column norms are the same order of magnitude.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* B = D^-1 A D
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* where D is a diagonal matrix
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This is necessary for the unsymmetric eigenvalue problem since the
|
|
Packit |
67cb25 |
* calculation can become numerically unstable for unbalanced
|
|
Packit |
67cb25 |
* matrices.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* See Golub & Van Loan, "Matrix Computations" (3rd ed), Section 7.5.7
|
|
Packit |
67cb25 |
* and Wilkinson & Reinsch, "Handbook for Automatic Computation", II/11 p320.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <stdlib.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_vector.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_matrix.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_blas.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <gsl/gsl_linalg.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#define FLOAT_RADIX 2.0
|
|
Packit |
67cb25 |
#define FLOAT_RADIX_SQ (FLOAT_RADIX * FLOAT_RADIX)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_balance_matrix(gsl_matrix * A, gsl_vector * D)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t N = A->size1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (N != D->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("vector must match matrix size", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double row_norm,
|
|
Packit |
67cb25 |
col_norm;
|
|
Packit |
67cb25 |
int not_converged;
|
|
Packit |
67cb25 |
gsl_vector_view v;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* initialize D to the identity matrix */
|
|
Packit |
67cb25 |
gsl_vector_set_all(D, 1.0);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
not_converged = 1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
while (not_converged)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t i, j;
|
|
Packit |
67cb25 |
double g, f, s;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
not_converged = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < N; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
row_norm = 0.0;
|
|
Packit |
67cb25 |
col_norm = 0.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (j = 0; j < N; ++j)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (j != i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
col_norm += fabs(gsl_matrix_get(A, j, i));
|
|
Packit |
67cb25 |
row_norm += fabs(gsl_matrix_get(A, i, j));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if ((col_norm == 0.0) || (row_norm == 0.0))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
continue;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
g = row_norm / FLOAT_RADIX;
|
|
Packit |
67cb25 |
f = 1.0;
|
|
Packit |
67cb25 |
s = col_norm + row_norm;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* find the integer power of the machine radix which
|
|
Packit |
67cb25 |
* comes closest to balancing the matrix
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
while (col_norm < g)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
f *= FLOAT_RADIX;
|
|
Packit |
67cb25 |
col_norm *= FLOAT_RADIX_SQ;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
g = row_norm * FLOAT_RADIX;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
while (col_norm > g)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
f /= FLOAT_RADIX;
|
|
Packit |
67cb25 |
col_norm /= FLOAT_RADIX_SQ;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if ((row_norm + col_norm) < 0.95 * s * f)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
not_converged = 1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
g = 1.0 / f;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* apply similarity transformation D, where
|
|
Packit |
67cb25 |
* D_{ij} = f_i * delta_{ij}
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* multiply by D^{-1} on the left */
|
|
Packit |
67cb25 |
v = gsl_matrix_row(A, i);
|
|
Packit |
67cb25 |
gsl_blas_dscal(g, &v.vector);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* multiply by D on the right */
|
|
Packit |
67cb25 |
v = gsl_matrix_column(A, i);
|
|
Packit |
67cb25 |
gsl_blas_dscal(f, &v.vector);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* keep track of transformation */
|
|
Packit |
67cb25 |
gsl_vector_set(D, i, gsl_vector_get(D, i) * f);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
} /* gsl_linalg_balance_matrix() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_linalg_balance_accum()
|
|
Packit |
67cb25 |
Accumulate a balancing transformation into a matrix.
|
|
Packit |
67cb25 |
This is used during the computation of Schur vectors since the
|
|
Packit |
67cb25 |
Schur vectors computed are the vectors for the balanced matrix.
|
|
Packit |
67cb25 |
We must at some point accumulate the balancing transformation into
|
|
Packit |
67cb25 |
the Schur vector matrix to get the vectors for the original matrix.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A -> D A
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where D is the diagonal matrix
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: A - matrix to transform
|
|
Packit |
67cb25 |
D - vector containing diagonal elements of D
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_balance_accum(gsl_matrix *A, gsl_vector *D)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t N = A->size1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (N != D->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("vector must match matrix size", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
double s;
|
|
Packit |
67cb25 |
gsl_vector_view r;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < N; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
s = gsl_vector_get(D, i);
|
|
Packit |
67cb25 |
r = gsl_matrix_row(A, i);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_blas_dscal(s, &r.vector);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
} /* gsl_linalg_balance_accum() */
|