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
#		All Rights Reserved
	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
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-