htswanalysis/MACS/lib/gsl/gsl-1.11/eigen/jacobi.c

   1 /* eigen/jacobi.c
   2  *
   3  * Copyright (C) 2004, 2007 Brian Gough, Gerard Jungman
   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 #include <config.h>
  21 #include <stdlib.h>
  22 #include <gsl/gsl_math.h>
  23 #include <gsl/gsl_vector.h>
  24 #include <gsl/gsl_matrix.h>
  25 #include <gsl/gsl_eigen.h>
  26
  27 /* Algorithm 8.4.3 - Cyclic Jacobi.  Golub & Van Loan, Matrix Computations */
  28
  29 static inline double
  30 symschur2 (gsl_matrix * A, size_t p, size_t q, double *c, double *s)
  31 {
  32   double Apq = gsl_matrix_get (A, p, q);
  33
  34   if (Apq != 0.0)
  35     {
  36       double App = gsl_matrix_get (A, p, p);
  37       double Aqq = gsl_matrix_get (A, q, q);
  38       double tau = (Aqq - App) / (2.0 * Apq);
  39       double t, c1;
  40
  41       if (tau >= 0.0)
  42         {
  43           t = 1.0 / (tau + hypot (1.0, tau));
  44         }
  45       else
  46         {
  47           t = -1.0 / (-tau + hypot (1.0, tau));
  48         }
  49
  50       c1 = 1.0 / hypot (1.0, t);
  51
  52       *c = c1;
  53       *s = t * c1;
  54     }
  55   else
  56     {
  57       *c = 1.0;
  58       *s = 0.0;
  59     }
  60
  61   /* reduction in off(A) is 2*(A_pq)^2 */
  62
  63   return fabs (Apq);
  64 }
  65
  66 inline static void
  67 apply_jacobi_L (gsl_matrix * A, size_t p, size_t q, double c, double s)
  68 {
  69   size_t j;
  70   const size_t N = A->size2;
  71
  72   /* Apply rotation to matrix A,  A' = J^T A */
  73
  74   for (j = 0; j < N; j++)
  75     {
  76       double Apj = gsl_matrix_get (A, p, j);
  77       double Aqj = gsl_matrix_get (A, q, j);
  78       gsl_matrix_set (A, p, j, Apj * c - Aqj * s);
  79       gsl_matrix_set (A, q, j, Apj * s + Aqj * c);
  80     }
  81 }
  82
  83 inline static void
  84 apply_jacobi_R (gsl_matrix * A, size_t p, size_t q, double c, double s)
  85 {
  86   size_t i;
  87   const size_t M = A->size1;
  88
  89   /* Apply rotation to matrix A,  A' = A J */
  90
  91   for (i = 0; i < M; i++)
  92     {
  93       double Aip = gsl_matrix_get (A, i, p);
  94       double Aiq = gsl_matrix_get (A, i, q);
  95       gsl_matrix_set (A, i, p, Aip * c - Aiq * s);
  96       gsl_matrix_set (A, i, q, Aip * s + Aiq * c);
  97     }
  98 }
  99
 100 inline static double
 101 norm (gsl_matrix * A)
 102 {
 103   size_t i, j, M = A->size1, N = A->size2;
 104   double sum = 0.0, scale = 0.0, ssq = 1.0;
 105
 106   for (i = 0; i < M; i++)
 107     {
 108       for (j = 0; j < N; j++)
 109         {
 110           double Aij = gsl_matrix_get (A, i, j);
 111
 112           if (Aij != 0.0)
 113             {
 114               double ax = fabs (Aij);
 115
 116               if (scale < ax)
 117                 {
 118                   ssq = 1.0 + ssq * (scale / ax) * (scale / ax);
 119                   scale = ax;
 120                 }
 121               else
 122                 {
 123                   ssq += (ax / scale) * (ax / scale);
 124                 }
 125             }
 126
 127         }
 128     }
 129
 130   sum = scale * sqrt (ssq);
 131
 132   return sum;
 133 }
 134
 135 int
 136 gsl_eigen_jacobi (gsl_matrix * a,
 137                   gsl_vector * eval,
 138                   gsl_matrix * evec, unsigned int max_rot, unsigned int *nrot)
 139 {
 140   size_t i, p, q;
 141   const size_t M = a->size1, N = a->size2;
 142   double red, redsum = 0.0;
 143
 144   if (M != N)
 145     {
 146       GSL_ERROR ("eigenproblem requires square matrix", GSL_ENOTSQR);
 147     }
 148   else if (M != evec->size1 || M != evec->size2)
 149     {
 150       GSL_ERROR ("eigenvector matrix must match input matrix", GSL_EBADLEN);
 151     }
 152   else if (M != eval->size)
 153     {
 154       GSL_ERROR ("eigenvalue vector must match input matrix", GSL_EBADLEN);
 155     }
 156
 157   gsl_vector_set_zero (eval);
 158   gsl_matrix_set_identity (evec);
 159
 160   for (i = 0; i < max_rot; i++)
 161     {
 162       double nrm = norm (a);
 163
 164       if (nrm == 0.0)
 165         break;
 166
 167       for (p = 0; p < N; p++)
 168         {
 169           for (q = p + 1; q < N; q++)
 170             {
 171               double c, s;
 172
 173               red = symschur2 (a, p, q, &c, &s);
 174               redsum += red;
 175
 176               /* Compute A <- J^T A J */
 177               apply_jacobi_L (a, p, q, c, s);
 178               apply_jacobi_R (a, p, q, c, s);
 179
 180               /* Compute V <- V J */
 181               apply_jacobi_R (evec, p, q, c, s);
 182             }
 183         }
 184     }
 185
 186   *nrot = i;
 187
 188   for (p = 0; p < N; p++)
 189     {
 190       double ep = gsl_matrix_get (a, p, p);
 191       gsl_vector_set (eval, p, ep);
 192     }
 193
 194   if (i == max_rot)
 195     {
 196       return GSL_EMAXITER;
 197     }
 198
 199   return GSL_SUCCESS;
 200 }
 201
 202 int
 203 gsl_eigen_invert_jacobi (const gsl_matrix * a,
 204                          gsl_matrix * ainv, unsigned int max_rot)
 205 {
 206   if (a->size1 != a->size2 || ainv->size1 != ainv->size2)
 207     {
 208       GSL_ERROR("jacobi method requires square matrix", GSL_ENOTSQR);
 209     }
 210   else if (a->size1 != ainv->size2)
 211     {
 212      GSL_ERROR ("inverse matrix must match input matrix", GSL_EBADLEN);
 213     }
 214
 215   {
 216     const size_t n = a->size2;
 217     size_t i,j,k;
 218     unsigned int nrot = 0;
 219     int status;
 220
 221     gsl_vector * eval = gsl_vector_alloc(n);
 222     gsl_matrix * evec = gsl_matrix_alloc(n, n);
 223     gsl_matrix * tmp = gsl_matrix_alloc(n, n);
 224
 225     gsl_matrix_memcpy (tmp, a);
 226
 227     status = gsl_eigen_jacobi(tmp, eval, evec, max_rot, &nrot);
 228
 229     for(i=0; i<n; i++)
 230       {
 231         for(j=0; j<n; j++)
 232           {
 233             double ainv_ij = 0.0;
 234
 235             for(k = 0; k<n; k++)
 236               {
 237                 double f = 1.0 / gsl_vector_get(eval, k);
 238                 double vik = gsl_matrix_get (evec, i, k);
 239                 double vjk = gsl_matrix_get (evec, j, k);
 240                 ainv_ij += vik * vjk * f;
 241               }
 242             gsl_matrix_set (ainv, i, j, ainv_ij);
 243           }
 244       }
 245
 246     gsl_vector_free(eval);
 247     gsl_matrix_free(evec);
 248     gsl_matrix_free(tmp);
 249
 250     if (status)
 251       {
 252         return status;
 253       }
 254     else
 255       {
 256         return GSL_SUCCESS;
 257       }
 258   }
 259 }