1### ===== File part of R package expm =====
2###
3### Function to compute the matrix exponential
4###
5###    exp(M) = sum(n = 0:Inf; M^n / n!),
6###
7### where M is an (n x n) matrix.
8###
9
10expm.s.Pade.s <- function(x, order, n=nrow(x)) {
11    ## no checking here; this is not to be called by the user
12    ## try have this work with "mpfrMatrix" <==> solve(<mpfrArray>)
13
14    ## s := the number of [s]quarings
15    e <- ceiling(log2(max(rowSums(abs(x)))))
16    s <- max(e+1, 0)
17
18    ## preconditions x :
19    x <- x / (2^s)
20
21    ## Pade approximation for exp(x)
22    c <- .5
23    D <- diag(1, n)
24    E <- D + c*x
25    D <- D - c*x
26    X <- x
27    p <- TRUE
28    for(k in 2:order) {
29	c <- c * (order-k+1) / (k*(2*order-k+1))
30	X <- x %*% X # now  X = x ^ k
31	cX <- c*X
32	E <- E + cX
33	D <- if(p) D + cX  else	 D - cX
34	p <- !p
35    }
36    E <- solve(D, E)
37
38    ## Undo the scaling by repeated squaring :
39    for(k in seq_len(s))
40	E <- E %*% E
41    E
42}
43
44
45expm.methSparse <- c("Higham08", "R_Eigen", "R_Pade")
46## keep this list up-to-date - test by setting  R_EXPM_NO_DENSE_COERCION
47## but NOTE: It may make sense to *keep* the message() about coercion to dense (memory blow up!)
48
49expm <- function(x, method = c("Higham08.b", "Higham08",
50                    "AlMohy-Hi09",
51		    "Ward77", "PadeRBS", "Pade", "Taylor", "PadeO", "TaylorO",
52		    "R_Eigen", "R_Pade", "R_Ward77", "hybrid_Eigen_Ward"),
53		 order = 8,
54		 trySym = TRUE, tol = .Machine$double.eps,
55		 do.sparseMsg = TRUE,
56		 preconditioning = c("2bal", "1bal", "buggy"))
57{
58    ## some methods work for "matrix", "Matrix", or "mpfrMatrix", iff solve(.,.) worked:
59    stopifnot(is.numeric(x) || (isM <- inherits(x, "dMatrix")) || inherits(x, "mpfrMatrix"))
60    if(length(d <- dim(x)) != 2) stop("argument is not a matrix")
61    if (d[1] != d[2]) stop("matrix not square")
62    method <- match.arg(method)
63    checkSparse <- !nzchar(Sys.getenv("R_EXPM_NO_DENSE_COERCION"))
64    isM <- !is.numeric(x) && isM
65    if(isM && checkSparse) { # i.e., a "dMatrix"
66	if(!(method %in% expm.methSparse) && is(x, "sparseMatrix")) {
67	    if(do.sparseMsg)
68		message("coercing to dense matrix, as required by method ",
69			dQuote(method))
70	    x <- as(x, "denseMatrix")
71	}
72    }
73    switch(method,
74           "AlMohy-Hi09" = expm.AlMoHi09(x, p = order)
75          ,
76	   "Higham08.b" = expm.Higham08(x, balancing = TRUE)
77	  ,
78	   "Higham08"	= expm.Higham08(x, balancing = FALSE)
79          ,
80	   "Ward77" = {
81	       ## AUTHORS: Christophe Dutang, Vincent Goulet at act ulaval ca
82	       ##	 built on "Matrix" package, built on 'octave' code
83	       ##	 Martin Maechler, for the preconditioning etc
84               if(!is.numeric(x)) x <- as(x, "matrix")
85	       switch(match.arg(preconditioning),
86		      "2bal" = .Call(do_expm, x, "Ward77"),
87		      "1bal" = .Call(do_expm, x, "Ward77_1"),
88		      "buggy"= .Call(do_expm, x, "buggy_Ward77"),
89		      stop("invalid 'preconditioning'"))
90	   },
91	   "R_Eigen" = {
92	       ## matrix exponential using eigenvalues / spectral decomposition :
93	       ## ==  Dubious Way 'Method 14' : is
94	       ## good for 'symmetric' or 'orthogonal' (or other 'normal' : A'A = AA' ):
95
96	       ## MM: improved from mexp2() with 'trySym' and isSymmetric()
97	       isSym <- if(trySym) isSymmetric.matrix(x) else FALSE
98	       z <- eigen(x, symmetric = isSym)
99	       V <- z$vectors
100	       Vi <- if(isSym) t(V) else solve(V)
101	       Re(V %*% (    exp(z$values)   *	Vi)) ## ==
102	       ##(V %*% diag(exp(z$values)) %*% Vi)
103	   },
104	   "hybrid_Eigen_Ward" = {
105	       ## AUTHOR: Christophe Dutang
106	       ## matrix exponential using eigenvalues / spectral decomposition and
107	       ## Ward(1977) algorithm if x is numerically non diagonalisable
108               if(!is.numeric(x)) x <- as(x, "matrix")
109	       .Call(do_expm_eigen, x, tol)
110	   },
111	   "R_Pade"= { ## use scaling + Pade + squaring with R code:
112
113	       ## matrix exponential using a scaling and squaring algorithm
114	       ## with a Pade approximation
115	       ## source code translated from expmdemo1.m in Matlab
116	       ## by Stig Mortensen <sbm@imm.dtu.dk>,
117	       ## prettified by MM -- works for "matrix" or "Matrix" matrices !
118
119	       stopifnot(order >= 2)
120	       expm.s.Pade.s(x, order, n=d[1])
121	   },
122	   "R_Ward77" = { ## R implementation of "Ward(1977)"
123	       ## also works for "Matrix" matrices
124	       stopifnot(order >= 2)
125	       n <- d[1]
126	       ## Preconditioning  Step 1: shift diagonal by average diagonal
127	       trShift <- sum(d.x <- diag(x))
128	       if(trShift) {
129		   trShift <- trShift/n
130		   diag(x) <- d.x - trShift
131	       }
132
133	       ## Preconditioning  Step 2: balancing with balance.
134	       ##		   ------
135	       ## For now, do like the octave implementation
136	       ## TODO later:  use "B" (faster; better condition of result)
137	       baP <- balance(x,     "P")
138	       baS <- balance(baP$z, "S")
139
140	       x <- expm.s.Pade.s(baS$z, order)
141	       ##   -------------  scaling + Pade + squaring ------
142	       ##   i.e., entails  Preconditioning Step 3 (and its reverse)
143
144	       ## Reverse step 2: ------------------
145	       ##
146	       ## Step 2 b:  apply inverse scaling
147	       d <- baS$scale
148	       x <- x * (d * rep(1/d, each = n))
149	       ##
150	       ## Step 2 a:  apply inverse permutation (of rows and columns):
151	       pp <- as.integer(baP$scale)
152	       if(baP$i1 > 1) {		    ## The lower part
153		   for(i in (baP$i1-1):1) { # 'p1' in *reverse* order
154		       tt <- x[,i]; x[,i] <- x[,pp[i]]; x[,pp[i]] <- tt
155		       tt <- x[i,]; x[i,] <- x[pp[i],]; x[pp[i],] <- tt
156		   }
157	       }
158	       if(baP$i2 < n) {		    ## The upper part
159		   for(i in (baP$i2+1):n) { # 'p2' in *forward* order
160		       ## swap	 i <-> pp[i]   both rows and columns
161		       tt <- x[,i]; x[,i] <- x[,pp[i]]; x[,pp[i]] <- tt
162		       tt <- x[i,]; x[i,] <- x[pp[i],]; x[pp[i],] <- tt
163		   }
164	       }
165
166	       ## reverse step 1 (diagonal shift)
167	       if(trShift) {
168		   exp(trShift) * x
169	       }
170	       else x
171	   },
172	   "PadeRBS" =  { ## the "expofit" method by  Roger B. Sidje (U.Queensland, AU)
173               if(!is.numeric(x)) x <- as(x, "matrix")
174	       stopifnot((order <- as.integer(order)) >= 1)
175	       if(!is.double(x)) storage.mode(x) <- "double"
176	       Fobj <- .Fortran(matexpRBS,
177				order,		   # IDEG  1
178				as.integer(d[1]),  # M	   2
179				T = 1.,		   # T	   3
180				H = x,		   # H	   4
181				iflag = integer(1) # IFLAG 5
182				)[c("H","iflag")]
183	       if(Fobj[["iflag"]] < 0)
184		   stop("Unable to determine matrix exponential")
185	       Fobj[["H"]]
186           }
187	   , { ## the "mexp" methods by
188	       ## AUTHORS: Marina Shapira and David Firth --------------
189               if(!is.numeric(x)) x <- as(x, "matrix")
190	       if(!is.double(x)) storage.mode(x) <- "double"
191	       order <- as.integer(order)
192	       ## MM:	a "silly"  way to code the method / order
193	       ntaylor <- npade <- 0L
194	       if (substr(method,1,4) == "Pade")
195		   npade <- order else ntaylor <- order
196	       res <- if(identical(grep("O$", method), 1L))
197			  .Fortran(matrexpO,
198			       X = x,
199			       size = d[1],
200			       ntaylor,
201			       npade,
202			       accuracy = double(1))[c("X", "accuracy")]
203		      else
204			  .Fortran(matrexp,
205			       X = x,
206			       size = d[1],
207			       ntaylor,
208			       npade,
209			       accuracy = double(1))[c("X", "accuracy")]
210	       structure(res$X, accuracy = res$accuracy)
211	   })## end{switch}
212}
213
214