htswanalysis/MACS/lib/gsl/gsl-1.11/linalg/balancemat.c

   1 /* linalg/balance.c
   2  *
   3  * Copyright (C) 2006 Patrick Alken
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 /* Balance a general matrix by scaling the rows and columns, so the
  21  * new row and column norms are the same order of magnitude.
  22  *
  23  * B =  D^-1 A D
  24  *
  25  * where D is a diagonal matrix
  26  *
  27  * This is necessary for the unsymmetric eigenvalue problem since the
  28  * calculation can become numerically unstable for unbalanced
  29  * matrices.
  30  *
  31  * See Golub & Van Loan, "Matrix Computations" (3rd ed), Section 7.5.7
  32  * and Wilkinson & Reinsch, "Handbook for Automatic Computation", II/11 p320.
  33  */
  34
  35 #include <config.h>
  36 #include <stdlib.h>
  37 #include <gsl/gsl_math.h>
  38 #include <gsl/gsl_vector.h>
  39 #include <gsl/gsl_matrix.h>
  40 #include <gsl/gsl_blas.h>
  41
  42 #include <gsl/gsl_linalg.h>
  43
  44 #define FLOAT_RADIX       2.0
  45 #define FLOAT_RADIX_SQ    (FLOAT_RADIX * FLOAT_RADIX)
  46
  47 int
  48 gsl_linalg_balance_matrix(gsl_matrix * A, gsl_vector * D)
  49 {
  50   const size_t N = A->size1;
  51
  52   if (N != D->size)
  53     {
  54       GSL_ERROR ("vector must match matrix size", GSL_EBADLEN);
  55     }
  56   else
  57     {
  58       double row_norm,
  59              col_norm;
  60       int not_converged;
  61       gsl_vector_view v;
  62
  63       /* initialize D to the identity matrix */
  64       gsl_vector_set_all(D, 1.0);
  65
  66       not_converged = 1;
  67
  68       while (not_converged)
  69         {
  70           size_t i, j;
  71           double g, f, s;
  72
  73           not_converged = 0;
  74
  75           for (i = 0; i < N; ++i)
  76             {
  77               row_norm = 0.0;
  78               col_norm = 0.0;
  79
  80               for (j = 0; j < N; ++j)
  81                 {
  82                   if (j != i)
  83                     {
  84                       col_norm += fabs(gsl_matrix_get(A, j, i));
  85                       row_norm += fabs(gsl_matrix_get(A, i, j));
  86                     }
  87                 }
  88
  89               if ((col_norm == 0.0) || (row_norm == 0.0))
  90                 {
  91                   continue;
  92                 }
  93
  94               g = row_norm / FLOAT_RADIX;
  95               f = 1.0;
  96               s = col_norm + row_norm;
  97
  98               /*
  99                * find the integer power of the machine radix which
 100                * comes closest to balancing the matrix
 101                */
 102               while (col_norm < g)
 103                 {
 104                   f *= FLOAT_RADIX;
 105                   col_norm *= FLOAT_RADIX_SQ;
 106                 }
 107
 108               g = row_norm * FLOAT_RADIX;
 109
 110               while (col_norm > g)
 111                 {
 112                   f /= FLOAT_RADIX;
 113                   col_norm /= FLOAT_RADIX_SQ;
 114                 }
 115
 116               if ((row_norm + col_norm) < 0.95 * s * f)
 117                 {
 118                   not_converged = 1;
 119
 120                   g = 1.0 / f;
 121
 122                   /*
 123                    * apply similarity transformation D, where
 124                    * D_{ij} = f_i * delta_{ij}
 125                    */
 126
 127                   /* multiply by D^{-1} on the left */
 128                   v = gsl_matrix_row(A, i);
 129                   gsl_blas_dscal(g, &v.vector);
 130
 131                   /* multiply by D on the right */
 132                   v = gsl_matrix_column(A, i);
 133                   gsl_blas_dscal(f, &v.vector);
 134
 135                   /* keep track of transformation */
 136                   gsl_vector_set(D, i, gsl_vector_get(D, i) * f);
 137                 }
 138             }
 139         }
 140
 141       return GSL_SUCCESS;
 142     }
 143 } /* gsl_linalg_balance_matrix() */
 144
 145 /*
 146 gsl_linalg_balance_accum()
 147   Accumulate a balancing transformation into a matrix.
 148 This is used during the computation of Schur vectors since the
 149 Schur vectors computed are the vectors for the balanced matrix.
 150 We must at some point accumulate the balancing transformation into
 151 the Schur vector matrix to get the vectors for the original matrix.
 152
 153 A -> D A
 154
 155 where D is the diagonal matrix
 156
 157 Inputs: A - matrix to transform
 158         D - vector containing diagonal elements of D
 159 */
 160
 161 int
 162 gsl_linalg_balance_accum(gsl_matrix *A, gsl_vector *D)
 163 {
 164   const size_t N = A->size1;
 165
 166   if (N != D->size)
 167     {
 168       GSL_ERROR ("vector must match matrix size", GSL_EBADLEN);
 169     }
 170   else
 171     {
 172       size_t i;
 173       double s;
 174       gsl_vector_view r;
 175
 176       for (i = 0; i < N; ++i)
 177         {
 178           s = gsl_vector_get(D, i);
 179           r = gsl_matrix_row(A, i);
 180
 181           gsl_blas_dscal(s, &r.vector);
 182         }
 183
 184       return GSL_SUCCESS;
 185     }
 186 } /* gsl_linalg_balance_accum() */