# R-alpha: matrix exponential

Robert Gentleman (rgentlem@stat.auckland.ac.nz)
Thu, 21 Nov 1996 09:11:57 +1300 (NZDT)

```From: Robert Gentleman <rgentlem@stat.auckland.ac.nz>
Date: Thu, 21 Nov 1996 09:11:57 +1300 (NZDT)
Message-Id: <199611202011.JAA22700@stat13.stat.auckland.ac.nz>
To: r-testers@stat.math.ethz.ch
Subject: R-alpha: matrix exponential

If you don't have a symmetric matrix (as you don't with most of the
stochastic processes I have worked with) and you have linpack you can use
the following
eddcmp<-function(inmat)
{
#	Copyright 1993 by Robert Gentleman and Jane Lindsey
if(!is.matrix(inmat)) stop("eddcmp can only be used for matrices")
n <- nrow(inmat)
if(n != ncol(inmat))
stop("eddcmp can only be used for square matrices")
emat <- matrix(0, nrow = n, ncol = n)
eval1 <- vector(mode = "double", length = n)
eval2 <- vector(mode = "double", length = n)
tivec <- vector(mode = "integer", length = n)
trvec <- vector(mode = "double", length = n)
ierr <- 0
z <- .Fortran("eddcmp",
as.matrix(inmat),
as.matrix(emat),
as.vector(eval1),
as.vector(eval2),
as.integer(n),
as.vector(tivec),
as.vector(trvec),
as.integer(ierr))
if(z[[8]] != 0)
warning("maximum number of iterations exceeded in fortran")
if(max(abs(z[[4]])) != 0)
warning("there are complex eigenvalues")
if(sum(duplicated(z[[3]])) > 0)
warning("some eigenvalues are duplicated")
return(list(evectors = z[[2]], evalues = z[[3]], im.evalues = z[[4]]))

}

#######the fortran

SUBROUTINE EDDCMP(MMAT,EMAT,EVAL1,EVAL2,N,IV1,FV1,IERR)
INTEGER N,IS1,IS2,IV1(N),IERR
DOUBLE PRECISION MMAT(N,N),EMAT(N,N),EVAL1(N),EVAL2(N),FV1(N)

CALL BALANC(N,N,MMAT,IS1,IS2,FV1)
CALL ELMHES(N,N,IS1,IS2,MMAT,IV1)
CALL ELTRAN(N,N,IS1,IS2,MMAT,IV1,EMAT)
CALL HQR2(N,N,IS1,IS2,MMAT,EVAL1,EVAL2,EMAT,IERR)
CALL BALBAK(N,N,IS1,IS2,FV1,N,EMAT)
RETURN
END
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
r-testers mailing list -- For info or help, send "info" or "help",
To [un]subscribe, send "[un]subscribe"
(in the "body", not the subject !)  To: r-testers-request@stat.math.ethz.ch
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
```