balance.c

00001 #include "solve_poly_root.h"
00002 
00003 
00004 #define RADIX 2
00005 #define RADIX2 (RADIX*RADIX)
00006 
00007 void
00008 balance_companion_matrix (double *m, size_t nc)
00009 {
00010   int not_converged = 1;
00011 
00012   double row_norm = 0;
00013   double col_norm = 0;
00014 
00015   while (not_converged)
00016     {
00017       size_t i, j;
00018       double g, f, s;
00019 
00020       not_converged = 0;
00021 
00022       for (i = 0; i < nc; i++)
00023     {
00024       /* column norm, excluding the diagonal */
00025 
00026       if (i != nc - 1)
00027         {
00028           col_norm = fabs (MAT (m, i + 1, i, nc));
00029         }
00030       else
00031         {
00032           col_norm = 0;
00033 
00034           for (j = 0; j < nc - 1; j++)
00035         {
00036           col_norm += fabs (MAT (m, j, nc - 1, nc));
00037         }
00038         }
00039 
00040       /* row norm, excluding the diagonal */
00041 
00042       if (i == 0)
00043         {
00044           row_norm = fabs (MAT (m, 0, nc - 1, nc));
00045         }
00046       else if (i == nc - 1)
00047         {
00048           row_norm = fabs (MAT (m, i, i - 1, nc));
00049         }
00050       else
00051         {
00052           row_norm = (fabs (MAT (m, i, i - 1, nc)) 
00053                           + fabs (MAT (m, i, nc - 1, nc)));
00054         }
00055 
00056       if (col_norm == 0 || row_norm == 0)
00057         {
00058           continue;
00059         }
00060 
00061       g = row_norm / RADIX;
00062       f = 1;
00063       s = col_norm + row_norm;
00064 
00065       while (col_norm < g)
00066         {
00067           f *= RADIX;
00068           col_norm *= RADIX2;
00069         }
00070 
00071       g = row_norm * RADIX;
00072 
00073       while (col_norm > g)
00074         {
00075           f /= RADIX;
00076           col_norm /= RADIX2;
00077         }
00078 
00079       if ((row_norm + col_norm) < 0.95 * s * f)
00080         {
00081           not_converged = 1;
00082 
00083           g = 1 / f;
00084 
00085           if (i == 0)
00086         {
00087           MAT (m, 0, nc - 1, nc) *= g;
00088         }
00089           else
00090         {
00091           MAT (m, i, i - 1, nc) *= g;
00092           MAT (m, i, nc - 1, nc) *= g;
00093         }
00094 
00095           if (i == nc - 1)
00096         {
00097           for (j = 0; j < nc; j++)
00098             {
00099               MAT (m, j, i, nc) *= f;
00100             }
00101         }
00102           else
00103         {
00104           MAT (m, i + 1, i, nc) *= f;
00105         }
00106         }
00107     }
00108     }
00109 }

Generated on Wed Oct 26 13:08:51 2005 for SINFONI Pipeline Reference Manual by doxygen1.2.13.1 written by Dimitri van Heesch, © 1997-2001