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]]

