Hi Martin

I have changed the implementation slightly so that it now handles complex
matrices as well. Take a look and see what you think. I realise that a lot
of the if/else mode checking could be merged.

Cheers
Rory

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;

    // necessary?
    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));

    if (mode == REALSXP)
        Memcpy(REAL(x_), REAL(x), (size_t)nrows*ncols);
    else
        Memcpy(COMPLEX(x_), COMPLEX(x), (size_t)nrows*ncols);

    // Initialize matrix to identity matrix
    // Set x[i * ncols + i] = 1
    if (mode == REALSXP)
        for (i = 0; i < ncols*nrows; i++)
            REAL(matrix)[i] = ((i % (ncols+1) == 0) ? 1 : 0);
    else
        for (i = 0; i < ncols*nrows;i++) {
            COMPLEX(matrix)[i].i = 0.0;
            COMPLEX(matrix)[i].r = ((i % (ncols+1) == 0) ? 1.0 : 0.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(matrix), nrows, ncols,
                        COMPLEX(x_), nrows, ncols, COMPLEX(tmp));

            //copyMatrixData(tmp, matrix, nrows, ncols, mode);
            if (mode == REALSXP)
                Memcpy(REAL(matrix), REAL(tmp), (size_t)nrows*ncols);
            else
                Memcpy(COMPLEX(matrix), COMPLEX(tmp), (size_t)nrows*ncols);
            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);
        if (mode == REALSXP)
            Memcpy(REAL(x_), REAL(x__), (size_t)nrows*ncols);
        else
            Memcpy(COMPLEX(x_), COMPLEX(x__), (size_t)nrows*ncols);

        e >>= 1;
    }

    PROTECT(dims2 = allocVector(INTSXP, 2));
    INTEGER(dims2)[0] = nrows;
    INTEGER(dims2)[1] = ncols;
    setAttrib(matrix, R_DimSymbol, dims2);

    UNPROTECT(5);
    return matrix;
}



On Sun, Apr 6, 2008 at 12:01 PM, Rory Winston <rory.winston@gmail.com>
wrote:

> Hi Martin
>
> Thanks for the detailed reply. I had a look at the matrix power
> implementation in the actuar package and the modified version in the expm
> package. I have a couple of questions/comments:
>
> 1. Firstly, I seem to have trouble loading expm.
>
> > install.packages("expm",repos="http://R-Forge.R-project.org")
> trying URL '
> http://R-Forge.R-project.org/bin/macosx/universal/contrib/2.6/expm_0.9-1.tgz
> '
> Content type 'application/x-gzip' length 149679 bytes (146 Kb)
> opened URL
> ==================================================
> downloaded 146 Kb
> ...
> > library("expm")
> Error in namespaceExport(ns, exports) : undefined exports :matpow
> Error: package/namespace load failed for 'expm'
>
> Possibly a namespace file issue? My version is:
> platform       i386-apple-darwin8.10.1
> arch           i386
> os             darwin8.10.1
> system         i386, darwin8.10.1
> status
> major          2
> minor          6.1
> year           2007
> month          11
> day            26
> svn rev        43537
> language       R
> version.string R version 2.6.1 (2007-11-26)
> >
>
> 2. On to the package implementation, I see we actually have very similar
> implementations. The main differences are:
>
>  i) For an exponent equal to -1, I call back into R for the solve()
> function using eval() and CAR/CDR'ing the arguments into place. The actuar
> package calls dgesv() directly. I suspect that the direct route is more
> efficient and thus the more preferable one. I see that your implementation
> doesnt calculate the inverse for an exponent of -1,is there a specific
> reason for doing that?
> ii) Regarding complex matrices: I guess we should have support for this,
> as its not unreasonable that someone may do this, and it should be pretty
> easy to implement. My function doesnt have full support yet.
> iii) A philosophical question - where the the "right" place for the %^%
> operator? Is it in the operator list at a C level along with %*% and the
> like? Or is it in an R file as a function definition? I dont really have a
> preference either way...have you an opinion on this?
>
> Thanks
> Rory
>
>
>
> On Sat, Apr 5, 2008 at 6:52 PM, Martin Maechler <
> maechler@stat.math.ethz.ch> wrote:
>
> > >>>>> "RW" == Rory Winston <rory.winston@gmail.com>
> > >>>>>     on Sat, 5 Apr 2008 14:44:44 +0100 writes:
> >
> >    RW> Hi all I recently started to write a matrix
> >    RW> exponentiation operator for R (by adding a new operator
> >    RW> definition to names.c, and adding the following code to
> >    RW> arrays.c). It is not finished yet, but I would like to
> >    RW> solicit some comments, as there are a few areas of R's
> >    RW> internals that I am still feeling my way around.
> >
> >    RW> Firstly:
> >
> >    RW> 1) Would there be interest in adding a new operator %^%
> >    RW> that performs the matrix equivalent of the scalar ^
> >    RW> operator?
> >
> > Yes. A few weeks ago, I had investigated a bit about this and
> > found several R-help/R-devel postings with code proposals
> > and then about half dozen CRAN packages with diverse
> > implementations of the matrix power (I say "power" very much on
> > purpose, in order to not confuse it with the matrix exponential
> > which is another much more interesting topic, also recently (+-
> > two months?) talked about.
> >
> > Consequently I made a few timing tests and found that indeed,
> > the "smart matrix power" {computing m^2, m^4, ... and only those
> > multiplications needed} as you find it in many good books about
> > algorithms and e.g. also in *the* standard Golub & Van Loan
> > "Matrix Computation" is better than "the eigen" method even for
> > large powers.
> >
> > matPower <- function(X,n)
> > ## Function to calculate the n-th power of a matrix X
> > {
> >    if(n != round(n)) {
> >        n <- round(n)
> >        warning("rounding exponent `n' to", n)
> >    }
> >    if(n == 0)
> >        return(diag(nrow = nrow(X)))
> >    n <- n - 1
> >    phi <- X
> >    ## pot <- X # the first power of the matrix.
> >    while (n > 0)
> >    {
> >        if (n %% 2)
> >            phi <- phi %*% X
> >        if (n == 1) break
> >        n <- n %/% 2
> >        X <- X %*% X
> >    }
> >    return(phi)
> > }
> >
> > "Simultaneously" people where looking at the matrix exponential
> > expm() in the Matrix package,
> > and some of us had consequently started the 'expm' project on
> > R-forge.
> > The main goal there has been to investigate several algorithms
> > for the matrix exponential, notably the one buggy implementation
> > (in the 'Matrix' package until a couple of weeks ago, the bug
> >  stemming from 'octave' implementation).
> > The authors of 'actuar', Vincent and Christophe, notably also
> > had code for the matrix *power* in a C (building on BLAS) and I
> > had added an R-interface '%^%'  there as well.
> >
> > Yes, with the goal to move that (not the matrix exponential yet)
> > into standard R.
> > Even though it's not used so often (in percentage of all uses of
> > R), it's simple to *right*, and I have seen very many versions
> > of the matrix power that were much slower / inaccurate / ...
> > such that a reference implementation seems to be called for.
> >
> > So, please consider
> >
> >    install.packages("expm",repos="http://R-Forge.R-project.org")
> >
> > -- but only from tomorrow for Windows (which installs a
> >   pre-compiled package), since I found that we had accidentally
> >   broken the package trivially by small changes two weeks ago.
> >
> > and then
> >
> >    library(expm)
> >    ?%^%
> >
> >
> > Best regards,
> > Martin Maechler, ETH Zurich
> >
> >
> >
> >
> >    RW> operator? I am implicitly assuming that the benefits of
> >    RW> a native exponentiation routine for Markov chain
> >    RW> evaluation or function generation would outstrip that of
> >    RW> an R solution. Based on my tests so far (comparing it to
> >    RW> a couple of different pure R versions) it does, but I
> >    RW> still there is much room for optimization in my method.
> >    RW> 2) Regarding the code below: Is there a better way to do
> >    RW> the matrix multiplication? I am creating quite a few
> >    RW> copies for this implementation of exponentiation by
> >    RW> squaring. Is there a way to cut down on the number of
> >    RW> copies I am making here (I am assuming that the lhs and
> >    RW> rhs of matprod() must be different instances).
> >
> >    RW> Any feedback appreciated !  Thanks Rory
> >
> >    RW> <snip>
> >
> >    RW> /* Convenience function */ static void
> >    RW> copyMatrixData(SEXP a, SEXP b, int nrows, int ncols, int
> >    RW> mode) { for (int i=0; i < ncols; ++i) for (int j=0; j <
> >    RW> nrows; ++j) REAL(b)[i * nrows + j] = REAL(a)[i * nrows +
> >    RW> j]; }
> >
> >    RW> SEXP do_matexp(SEXP call, SEXP op, SEXP args, SEXP rho)
> >    RW> { int nrows, ncols; SEXP matrix, tmp, dims, dims2; SEXP
> >    RW> x, y, x_, x__; int i,j,e,mode;
> >
> >    RW>     // Still need to fix full complex support mode =
> >    RW> isComplex(CAR(args)) ? CPLXSXP : REALSXP;
> >
> >    RW>     SETCAR(args, coerceVector(CAR(args), mode)); x =
> >    RW> CAR(args); y = CADR(args);
> >
> >    RW>     dims = getAttrib(x, R_DimSymbol); nrows =
> >    RW> INTEGER(dims)[0]; ncols = INTEGER(dims)[1];
> >
> >
> >    RW>     if (nrows != ncols) error(_("can only raise square
> >    RW> matrix to power"));
> >
> >    RW>     if (!isNumeric(y)) error(_("exponent must be a
> >    RW> scalar integer"));
> >
> >    RW>     e = asInteger(y);
> >
> >    RW>     if (e < -1) error(_("exponent must be >= -1")); else
> >    RW> if (e == 1) return x;
> >
> >    RW>     else if (e == -1) { /* return matrix inverse via
> >    RW> solve() */ SEXP p1, p2, inv; PROTECT(p1 = p2 =
> >    RW> allocList(2)); SET_TYPEOF(p1, LANGSXP); CAR(p2) =
> >    RW> install("solve.default"); p2 = CDR(p2); CAR(p2) = x; inv
> >    RW> = eval(p1, rho); UNPROTECT(1); return inv; }
> >
> >    RW>     PROTECT(matrix = allocVector(mode, nrows * ncols));
> >    RW> PROTECT(tmp = allocVector(mode, nrows * ncols));
> >    RW> PROTECT(x_ = allocVector(mode, nrows * ncols));
> >    RW> PROTECT(x__ = allocVector(mode, nrows * ncols));
> >
> >    RW>     copyMatrixData(x, x_, nrows, ncols, mode);
> >
> >    RW>     // Initialize matrix to identity matrix // Set x[i *
> >    RW> ncols + i] = 1 for (i = 0; i < ncols*nrows; i++)
> >    RW> REAL(matrix)[i] = ((i % (ncols+1) == 0) ? 1 : 0);
> >
> >    RW>     if (e == 0) { ; // return identity matrix } else
> >    RW> while (e > 0) { if (e & 1) { if (mode == REALSXP)
> >    RW> matprod(REAL(matrix), nrows, ncols, REAL(x_), nrows,
> >    RW> ncols, REAL(tmp)); else cmatprod(COMPLEX(tmp), nrows,
> >    RW> ncols, COMPLEX(x_), nrows, ncols, COMPLEX(matrix));
> >
> >    RW>             copyMatrixData(tmp, matrix, nrows, ncols,
> >    RW> mode); e--; }
> >
> >    RW>         if (mode == REALSXP) matprod(REAL(x_), nrows,
> >    RW> ncols, REAL(x_), nrows, ncols, REAL(x__)); else
> >    RW> cmatprod(COMPLEX(x_), nrows, ncols, COMPLEX(x_), nrows,
> >    RW> ncols, COMPLEX(x__));
> >
> >    RW>         copyMatrixData(x__, x_, nrows, ncols, mode); e
> >    RW> /= 2; }
> >
> >    RW>     PROTECT(dims2 = allocVector(INTSXP, 2));
> >    RW> INTEGER(dims2)[0] = nrows; INTEGER(dims2)[1] = ncols;
> >    RW> setAttrib(matrix, R_DimSymbol, dims2);
> >
> >    RW>     UNPROTECT(5); return matrix; }
> >
> >    RW>         [[alternative HTML version deleted]]
> >
> >    RW> ______________________________________________
> >    RW> R-devel@r-project.org mailing list
> >    RW> https://stat.ethz.ch/mailman/listinfo/r-devel
> >
>
>

	[[alternative HTML version deleted]]

