Hi all
I recently started to write a matrix exponentiation operator for R (by
adding a new operator definition to names.c, and adding the following code
to arrays.c). It is not finished yet, but I would like to solicit some
comments, as there are a few areas of R's internals that I am still feeling
my way around.
Firstly:
1) Would there be interest in adding a new operator %^% that performs the
matrix equivalent of the scalar ^ operator? I am implicitly assuming that
the benefits of a native exponentiation routine for Markov chain evaluation
or function generation would outstrip that of an R solution. Based on my
tests so far (comparing it to a couple of different pure R versions) it
does, but I still there is much room for optimization in my method.
2) Regarding the code below: Is there a better way to do the matrix
multiplication? I am creating quite a few copies for this implementation of
exponentiation by squaring. Is there a way to cut down on the number of
copies I am making here (I am assuming that the lhs and rhs of matprod()
must be different instances).
Any feedback appreciated !
Thanks
Rory
<snip>
/* Convenience function */
static void copyMatrixData(SEXP a, SEXP b, int nrows, int ncols, int mode) {
for (int i=0; i < ncols; ++i)
for (int j=0; j < nrows; ++j)
REAL(b)[i * nrows + j] = REAL(a)[i * nrows + j];
}
SEXP do_matexp(SEXP call, SEXP op, SEXP args, SEXP rho)
{
int nrows, ncols;
SEXP matrix, tmp, dims, dims2;
SEXP x, y, x_, x__;
int i,j,e,mode;
// Still need to fix full complex support
mode = isComplex(CAR(args)) ? CPLXSXP : REALSXP;
SETCAR(args, coerceVector(CAR(args), mode));
x = CAR(args);
y = CADR(args);
dims = getAttrib(x, R_DimSymbol);
nrows = INTEGER(dims)[0];
ncols = INTEGER(dims)[1];
if (nrows != ncols)
error(_("can only raise square matrix to power"));
if (!isNumeric(y))
error(_("exponent must be a scalar integer"));
e = asInteger(y);
if (e < -1)
error(_("exponent must be >= -1"));
else if (e == 1)
return x;
else if (e == -1) { /* return matrix inverse via solve() */
SEXP p1, p2, inv;
PROTECT(p1 = p2 = allocList(2));
SET_TYPEOF(p1, LANGSXP);
CAR(p2) = install("solve.default");
p2 = CDR(p2);
CAR(p2) = x;
inv = eval(p1, rho);
UNPROTECT(1);
return inv;
}
PROTECT(matrix = allocVector(mode, nrows * ncols));
PROTECT(tmp = allocVector(mode, nrows * ncols));
PROTECT(x_ = allocVector(mode, nrows * ncols));
PROTECT(x__ = allocVector(mode, nrows * ncols));
copyMatrixData(x, x_, nrows, ncols, mode);
// Initialize matrix to identity matrix
// Set x[i * ncols + i] = 1
for (i = 0; i < ncols*nrows; i++)
REAL(matrix)[i] = ((i % (ncols+1) == 0) ? 1 : 0);
if (e == 0) {
; // return identity matrix
}
else
while (e > 0) {
if (e & 1) {
if (mode == REALSXP)
matprod(REAL(matrix), nrows, ncols,
REAL(x_), nrows, ncols, REAL(tmp));
else
cmatprod(COMPLEX(tmp), nrows, ncols,
COMPLEX(x_), nrows, ncols, COMPLEX(matrix));
copyMatrixData(tmp, matrix, nrows, ncols, mode);
e--;
}
if (mode == REALSXP)
matprod(REAL(x_), nrows, ncols,
REAL(x_), nrows, ncols, REAL(x__));
else
cmatprod(COMPLEX(x_), nrows, ncols,
COMPLEX(x_), nrows, ncols, COMPLEX(x__));
copyMatrixData(x__, x_, nrows, ncols, mode);
e /= 2;
}
PROTECT(dims2 = allocVector(INTSXP, 2));
INTEGER(dims2)[0] = nrows;
INTEGER(dims2)[1] = ncols;
setAttrib(matrix, R_DimSymbol, dims2);
UNPROTECT(5);
return matrix;
}
[[alternative HTML version deleted]]
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
