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