htswanalysis/MACS/lib/gsl/gsl-1.11/poly/balance.c

   1 /* poly/balance.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
   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 static void balance_companion_matrix (double *m, size_t n);
  21
  22 #define RADIX 2
  23 #define RADIX2 (RADIX*RADIX)
  24
  25 static void
  26 balance_companion_matrix (double *m, size_t nc)
  27 {
  28   int not_converged = 1;
  29
  30   double row_norm = 0;
  31   double col_norm = 0;
  32
  33   while (not_converged)
  34     {
  35       size_t i, j;
  36       double g, f, s;
  37
  38       not_converged = 0;
  39
  40       for (i = 0; i < nc; i++)
  41         {
  42           /* column norm, excluding the diagonal */
  43
  44           if (i != nc - 1)
  45             {
  46               col_norm = fabs (MAT (m, i + 1, i, nc));
  47             }
  48           else
  49             {
  50               col_norm = 0;
  51
  52               for (j = 0; j < nc - 1; j++)
  53                 {
  54                   col_norm += fabs (MAT (m, j, nc - 1, nc));
  55                 }
  56             }
  57
  58           /* row norm, excluding the diagonal */
  59
  60           if (i == 0)
  61             {
  62               row_norm = fabs (MAT (m, 0, nc - 1, nc));
  63             }
  64           else if (i == nc - 1)
  65             {
  66               row_norm = fabs (MAT (m, i, i - 1, nc));
  67             }
  68           else
  69             {
  70               row_norm = (fabs (MAT (m, i, i - 1, nc))
  71                           + fabs (MAT (m, i, nc - 1, nc)));
  72             }
  73
  74           if (col_norm == 0 || row_norm == 0)
  75             {
  76               continue;
  77             }
  78
  79           g = row_norm / RADIX;
  80           f = 1;
  81           s = col_norm + row_norm;
  82
  83           while (col_norm < g)
  84             {
  85               f *= RADIX;
  86               col_norm *= RADIX2;
  87             }
  88
  89           g = row_norm * RADIX;
  90
  91           while (col_norm > g)
  92             {
  93               f /= RADIX;
  94               col_norm /= RADIX2;
  95             }
  96
  97           if ((row_norm + col_norm) < 0.95 * s * f)
  98             {
  99               not_converged = 1;
 100
 101               g = 1 / f;
 102
 103               if (i == 0)
 104                 {
 105                   MAT (m, 0, nc - 1, nc) *= g;
 106                 }
 107               else
 108                 {
 109                   MAT (m, i, i - 1, nc) *= g;
 110                   MAT (m, i, nc - 1, nc) *= g;
 111                 }
 112
 113               if (i == nc - 1)
 114                 {
 115                   for (j = 0; j < nc; j++)
 116                     {
 117                       MAT (m, j, i, nc) *= f;
 118                     }
 119                 }
 120               else
 121                 {
 122                   MAT (m, i + 1, i, nc) *= f;
 123                 }
 124             }
 125         }
 126     }
 127 }