src/mat_expo.c

   1 /* matexpo.c       2007-10-08 */
   2
   3 /* Copyright 2007 Emmanuel Paradis
   4
   5 /* This file is part of the R-package `ape'. */
   6 /* See the file ../COPYING for licensing issues. */
   7
   8 #include <R.h>
   9 #include <R_ext/Lapack.h>
  10
  11 void mat_expo(double *P, int *nr)
  12 /* This function computes the exponential of a nr x nr matrix */
  13 {
  14         double *U, *vl, *WR, *Uinv, *WI, *work;
  15         int i, j, k, l, info, *ipiv, n = *nr, nc = n*n, lw = nc << 1, *ord;
  16         char yes = 'V', no = 'N';
  17
  18         U = (double *)R_alloc(nc, sizeof(double));
  19         vl = (double *)R_alloc(n, sizeof(double));
  20         WR = (double *)R_alloc(n, sizeof(double));
  21         Uinv = (double *)R_alloc(nc, sizeof(double));
  22         WI = (double *)R_alloc(n, sizeof(double));
  23         work = (double *)R_alloc(lw, sizeof(double));
  24
  25         ipiv = (int *)R_alloc(nc, sizeof(int));
  26         ord = (int *)R_alloc(n, sizeof(int));
  27
  28 /* The matrix is not symmetric, so we use 'dgeev'.
  29    We take the real part of the eigenvalues -> WR
  30    and the right eigenvectors (vr) -> U */
  31         F77_CALL(dgeev)(&no, &yes, &n, P, &n, WR, WI, vl, &n,
  32                         U, &n, work, &lw, &info);
  33
  34 /* It is not necessary to sort the eigenvalues...
  35    Copy U -> P */
  36         memcpy(P, U, nc*sizeof(double));
  37
  38 /* For the inversion, we first make Uinv an identity matrix */
  39         memset(Uinv, 0, nc*sizeof(double));
  40         for (i = 0; i < nc; i += n + 1) Uinv[i] = 1;
  41
  42 /* The matrix is not symmetric, so we use 'dgesv'.
  43    This subroutine puts the result in Uinv (B)
  44    (P [= U] is erased) */
  45         F77_CALL(dgesv)(&n, &n, P, &n, ipiv, Uinv, &n, &info);
  46
  47 /* The matrix product of U with the eigenvalues diagonal matrix: */
  48         for (i = 0; i < n; i++)
  49                 for (j = 0; j < n; j++)
  50                         U[j + i*n] *= exp(WR[i]);
  51
  52 /* The second matrix product with U^-1 */
  53         memset(P, 0, nc*sizeof(double));
  54
  55         for (k = 0; k < n; k++) {
  56                 for (l = 0; l < n; l++) {
  57                         lw = l + k*n;
  58                         for (i = 0 + n*k, j = l; j < nc; i++, j += n)
  59                                 P[lw] += U[j]*Uinv[i];
  60                 }
  61         }
  62 }