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

   1 /* linalg/exponential.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2007 Gerard Jungman, 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 /* Author:  G. Jungman */
  21
  22 /* Calculate the matrix exponential, following
  23  * Moler + Van Loan, SIAM Rev. 20, 801 (1978).
  24  */
  25
  26 #include <config.h>
  27 #include <stdlib.h>
  28 #include <gsl/gsl_math.h>
  29 #include <gsl/gsl_mode.h>
  30 #include <gsl/gsl_errno.h>
  31 #include <gsl/gsl_blas.h>
  32
  33 #include "gsl_linalg.h"
  34
  35
  36 /* store one of the suggested choices for the
  37  * Taylor series / square  method from Moler + VanLoan
  38  */
  39 struct moler_vanloan_optimal_suggestion
  40 {
  41   int k;
  42   int j;
  43 };
  44 typedef  struct moler_vanloan_optimal_suggestion  mvl_suggestion_t;
  45
  46
  47 /* table from Moler and Van Loan
  48  * mvl_tab[gsl_mode_t][matrix_norm_group]
  49  */
  50 static mvl_suggestion_t mvl_tab[3][6] =
  51 {
  52   /* double precision */
  53   {
  54     { 5, 1 }, { 5, 4 }, { 7, 5 }, { 9, 7 }, { 10, 10 }, { 8, 14 }
  55   },
  56
  57   /* single precision */
  58   {
  59     { 2, 1 }, { 4, 0 }, { 7, 1 }, { 6, 5 }, { 5, 9 }, { 7, 11 }
  60   },
  61
  62   /* approx precision */
  63   {
  64     { 1, 0 }, { 3, 0 }, { 5, 1 }, { 4, 5 }, { 4, 8 }, { 2, 11 }
  65   }
  66 };
  67
  68
  69 inline
  70 static double
  71 sup_norm(const gsl_matrix * A)
  72 {
  73   double min, max;
  74   gsl_matrix_minmax(A, &min, &max);
  75   return GSL_MAX_DBL(fabs(min), fabs(max));
  76 }
  77
  78
  79 static
  80 mvl_suggestion_t
  81 obtain_suggestion(const gsl_matrix * A, gsl_mode_t mode)
  82 {
  83   const unsigned int mode_prec = GSL_MODE_PREC(mode);
  84   const double norm_A = sup_norm(A);
  85   if(norm_A < 0.01) return mvl_tab[mode_prec][0];
  86   else if(norm_A < 0.1) return mvl_tab[mode_prec][1];
  87   else if(norm_A < 1.0) return mvl_tab[mode_prec][2];
  88   else if(norm_A < 10.0) return mvl_tab[mode_prec][3];
  89   else if(norm_A < 100.0) return mvl_tab[mode_prec][4];
  90   else if(norm_A < 1000.0) return mvl_tab[mode_prec][5];
  91   else
  92   {
  93     /* outside the table we simply increase the number
  94      * of squarings, bringing the reduced matrix into
  95      * the range of the table; this is obviously suboptimal,
  96      * but that is the price paid for not having those extra
  97      * table entries
  98      */
  99     const double extra = log(1.01*norm_A/1000.0) / M_LN2;
 100     const int extra_i = (unsigned int) ceil(extra);
 101     mvl_suggestion_t s = mvl_tab[mode][5];
 102     s.j += extra_i;
 103     return s;
 104   }
 105 }
 106
 107
 108 /* use series representation to calculate matrix exponential;
 109  * this is used for small matrices; we use the sup_norm
 110  * to measure the size of the terms in the expansion
 111  */
 112 static void
 113 matrix_exp_series(
 114   const gsl_matrix * B,
 115   gsl_matrix * eB,
 116   int number_of_terms
 117   )
 118 {
 119   int count;
 120   gsl_matrix * temp = gsl_matrix_calloc(B->size1, B->size2);
 121
 122   /* init the Horner polynomial evaluation,
 123    * eB = 1 + B/number_of_terms; we use
 124    * eB to collect the partial results
 125    */
 126   gsl_matrix_memcpy(eB, B);
 127   gsl_matrix_scale(eB, 1.0/number_of_terms);
 128   gsl_matrix_add_diagonal(eB, 1.0);
 129   for(count = number_of_terms-1; count >= 1; --count)
 130   {
 131     /*  mult_temp = 1 + B eB / count  */
 132     gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, B, eB, 0.0, temp);
 133     gsl_matrix_scale(temp, 1.0/count);
 134     gsl_matrix_add_diagonal(temp, 1.0);
 135
 136     /*  transfer partial result out of temp */
 137     gsl_matrix_memcpy(eB, temp);
 138   }
 139
 140   /* now eB holds the full result; we're done */
 141   gsl_matrix_free(temp);
 142 }
 143
 144
 145 int
 146 gsl_linalg_exponential_ss(
 147   const gsl_matrix * A,
 148   gsl_matrix * eA,
 149   gsl_mode_t mode
 150   )
 151 {
 152   if(A->size1 != A->size2)
 153   {
 154     GSL_ERROR("cannot exponentiate a non-square matrix", GSL_ENOTSQR);
 155   }
 156   else if(A->size1 != eA->size1 || A->size2 != eA->size2)
 157   {
 158     GSL_ERROR("exponential of matrix must have same dimension as matrix", GSL_EBADLEN);
 159   }
 160   else
 161   {
 162     int i;
 163     const mvl_suggestion_t sugg = obtain_suggestion(A, mode);
 164     const double divisor = exp(M_LN2 * sugg.j);
 165
 166     gsl_matrix * reduced_A = gsl_matrix_alloc(A->size1, A->size2);
 167
 168     /*  decrease A by the calculated divisor  */
 169     gsl_matrix_memcpy(reduced_A, A);
 170     gsl_matrix_scale(reduced_A, 1.0/divisor);
 171
 172     /*  calculate exp of reduced matrix; store in eA as temp  */
 173     matrix_exp_series(reduced_A, eA, sugg.k);
 174
 175     /*  square repeatedly; use reduced_A for scratch */
 176     for(i = 0; i < sugg.j; ++i)
 177     {
 178       gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, eA, eA, 0.0, reduced_A);
 179       gsl_matrix_memcpy(eA, reduced_A);
 180     }
 181
 182     gsl_matrix_free(reduced_A);
 183
 184     return GSL_SUCCESS;
 185   }
 186 }
 187