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

   1 /* eigen/nonsymm.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 #include <config.h>
  21 #include <stdlib.h>
  22 #include <math.h>
  23 #include <gsl/gsl_eigen.h>
  24 #include <gsl/gsl_linalg.h>
  25 #include <gsl/gsl_math.h>
  26 #include <gsl/gsl_blas.h>
  27 #include <gsl/gsl_vector.h>
  28 #include <gsl/gsl_vector_complex.h>
  29 #include <gsl/gsl_matrix.h>
  30
  31 /*
  32  * This module computes the eigenvalues of a real nonsymmetric
  33  * matrix, using the double shift Francis method.
  34  *
  35  * See the references in francis.c.
  36  *
  37  * This module gets the matrix ready by balancing it and
  38  * reducing it to Hessenberg form before passing it to the
  39  * francis module.
  40  */
  41
  42 /*
  43 gsl_eigen_nonsymm_alloc()
  44
  45 Allocate a workspace for solving the nonsymmetric eigenvalue problem.
  46 The size of this workspace is O(2n)
  47
  48 Inputs: n - size of matrix
  49
  50 Return: pointer to workspace
  51 */
  52
  53 gsl_eigen_nonsymm_workspace *
  54 gsl_eigen_nonsymm_alloc(const size_t n)
  55 {
  56   gsl_eigen_nonsymm_workspace *w;
  57
  58   if (n == 0)
  59     {
  60       GSL_ERROR_NULL ("matrix dimension must be positive integer",
  61                       GSL_EINVAL);
  62     }
  63
  64   w = (gsl_eigen_nonsymm_workspace *)
  65       calloc (1, sizeof (gsl_eigen_nonsymm_workspace));
  66
  67   if (w == 0)
  68     {
  69       GSL_ERROR_NULL ("failed to allocate space for workspace", GSL_ENOMEM);
  70     }
  71
  72   w->size = n;
  73   w->Z = NULL;
  74   w->do_balance = 0;
  75
  76   w->diag = gsl_vector_alloc(n);
  77
  78   if (w->diag == 0)
  79     {
  80       gsl_eigen_nonsymm_free(w);
  81       GSL_ERROR_NULL ("failed to allocate space for balancing vector", GSL_ENOMEM);
  82     }
  83
  84   w->tau = gsl_vector_alloc(n);
  85
  86   if (w->tau == 0)
  87     {
  88       gsl_eigen_nonsymm_free(w);
  89       GSL_ERROR_NULL ("failed to allocate space for hessenberg coefficients", GSL_ENOMEM);
  90     }
  91
  92   w->francis_workspace_p = gsl_eigen_francis_alloc();
  93
  94   if (w->francis_workspace_p == 0)
  95     {
  96       gsl_eigen_nonsymm_free(w);
  97       GSL_ERROR_NULL ("failed to allocate space for francis workspace", GSL_ENOMEM);
  98     }
  99
 100   return (w);
 101 } /* gsl_eigen_nonsymm_alloc() */
 102
 103 /*
 104 gsl_eigen_nonsymm_free()
 105   Free workspace w
 106 */
 107
 108 void
 109 gsl_eigen_nonsymm_free (gsl_eigen_nonsymm_workspace * w)
 110 {
 111   if (w->tau)
 112     gsl_vector_free(w->tau);
 113
 114   if (w->diag)
 115     gsl_vector_free(w->diag);
 116
 117   if (w->francis_workspace_p)
 118     gsl_eigen_francis_free(w->francis_workspace_p);
 119
 120   free(w);
 121 } /* gsl_eigen_nonsymm_free() */
 122
 123 /*
 124 gsl_eigen_nonsymm_params()
 125   Set some parameters which define how we solve the eigenvalue
 126 problem.
 127
 128 Inputs: compute_t - 1 if we want to compute T, 0 if not
 129         balance   - 1 if we want to balance the matrix, 0 if not
 130         w         - nonsymm workspace
 131 */
 132
 133 void
 134 gsl_eigen_nonsymm_params (const int compute_t, const int balance,
 135                           gsl_eigen_nonsymm_workspace *w)
 136 {
 137   gsl_eigen_francis_T(compute_t, w->francis_workspace_p);
 138   w->do_balance = balance;
 139 } /* gsl_eigen_nonsymm_params() */
 140
 141 /*
 142 gsl_eigen_nonsymm()
 143
 144 Solve the nonsymmetric eigenvalue problem
 145
 146 A x = \lambda x
 147
 148 for the eigenvalues \lambda using the Francis method.
 149
 150 Here we compute the real Schur form
 151
 152 T = Z^t A Z
 153
 154 with the diagonal blocks of T giving us the eigenvalues.
 155 Z is a matrix of Schur vectors which is not computed by
 156 this algorithm. See gsl_eigen_nonsymm_Z().
 157
 158 Inputs: A    - general real matrix
 159         eval - where to store eigenvalues
 160         w    - workspace
 161
 162 Return: success or error
 163
 164 Notes: If T is computed, it is stored in A on output. Otherwise
 165        the diagonal of A contains the 1-by-1 and 2-by-2 eigenvalue
 166        blocks.
 167 */
 168
 169 int
 170 gsl_eigen_nonsymm (gsl_matrix * A, gsl_vector_complex * eval,
 171                    gsl_eigen_nonsymm_workspace * w)
 172 {
 173   const size_t N = A->size1;
 174
 175   /* check matrix and vector sizes */
 176
 177   if (N != A->size2)
 178     {
 179       GSL_ERROR ("matrix must be square to compute eigenvalues", GSL_ENOTSQR);
 180     }
 181   else if (eval->size != N)
 182     {
 183       GSL_ERROR ("eigenvalue vector must match matrix size", GSL_EBADLEN);
 184     }
 185   else
 186     {
 187       int s;
 188
 189       if (w->do_balance)
 190         {
 191           /* balance the matrix */
 192           gsl_linalg_balance_matrix(A, w->diag);
 193         }
 194
 195       /* compute the Hessenberg reduction of A */
 196       gsl_linalg_hessenberg_decomp(A, w->tau);
 197
 198       if (w->Z)
 199         {
 200           /*
 201            * initialize the matrix Z to U, which is the matrix used
 202            * to construct the Hessenberg reduction.
 203            */
 204
 205           /* compute U and store it in Z */
 206           gsl_linalg_hessenberg_unpack(A, w->tau, w->Z);
 207
 208           /* find the eigenvalues and Schur vectors */
 209           s = gsl_eigen_francis_Z(A, eval, w->Z, w->francis_workspace_p);
 210
 211           if (w->do_balance)
 212             {
 213               /*
 214                * The Schur vectors in Z are the vectors for the balanced
 215                * matrix. We now must undo the balancing to get the
 216                * vectors for the original matrix A.
 217                */
 218               gsl_linalg_balance_accum(w->Z, w->diag);
 219             }
 220         }
 221       else
 222         {
 223           /* find the eigenvalues only */
 224           s = gsl_eigen_francis(A, eval, w->francis_workspace_p);
 225         }
 226
 227       w->n_evals = w->francis_workspace_p->n_evals;
 228
 229       return s;
 230     }
 231 } /* gsl_eigen_nonsymm() */
 232
 233 /*
 234 gsl_eigen_nonsymm_Z()
 235
 236 Solve the nonsymmetric eigenvalue problem
 237
 238 A x = \lambda x
 239
 240 for the eigenvalues \lambda.
 241
 242 Here we compute the real Schur form
 243
 244 T = Z^t A Z
 245
 246 with the diagonal blocks of T giving us the eigenvalues.
 247 Z is the matrix of Schur vectors.
 248
 249 Inputs: A    - general real matrix
 250         eval - where to store eigenvalues
 251         Z    - where to store Schur vectors
 252         w    - workspace
 253
 254 Return: success or error
 255
 256 Notes: If T is computed, it is stored in A on output. Otherwise
 257        the diagonal of A contains the 1-by-1 and 2-by-2 eigenvalue
 258        blocks.
 259 */
 260
 261 int
 262 gsl_eigen_nonsymm_Z (gsl_matrix * A, gsl_vector_complex * eval,
 263                      gsl_matrix * Z, gsl_eigen_nonsymm_workspace * w)
 264 {
 265   /* check matrix and vector sizes */
 266
 267   if (A->size1 != A->size2)
 268     {
 269       GSL_ERROR ("matrix must be square to compute eigenvalues", GSL_ENOTSQR);
 270     }
 271   else if (eval->size != A->size1)
 272     {
 273       GSL_ERROR ("eigenvalue vector must match matrix size", GSL_EBADLEN);
 274     }
 275   else if ((Z->size1 != Z->size2) || (Z->size1 != A->size1))
 276     {
 277       GSL_ERROR ("Z matrix has wrong dimensions", GSL_EBADLEN);
 278     }
 279   else
 280     {
 281       int s;
 282
 283       w->Z = Z;
 284
 285       s = gsl_eigen_nonsymm(A, eval, w);
 286
 287       w->Z = NULL;
 288
 289       return s;
 290     }
 291 } /* gsl_eigen_nonsymm_Z() */