The R Project SVN R-packages

#include "dgeMatrix.h"

SEXP dgeMatrix_validate(SEXP obj)
{
    SEXP x = GET_SLOT(obj, Matrix_xSym),
    Dim = GET_SLOT(obj, Matrix_DimSym),
    fact = GET_SLOT(obj, Matrix_factorSym),
    rc = GET_SLOT(obj, Matrix_rcondSym);
    int m, n;

    if (length(Dim) != 2)
    return mkString(_("Dim slot must have length 2"));
    m = INTEGER(Dim)[0]; n = INTEGER(Dim)[1];
    if (m < 0 || n < 0)
    return mkString(_("Negative value(s) in Dim"));
    if (length(x) != m * n)
    return mkString(_("length of x slot != prod(Dim)"));
    if (!isReal(x))
    return mkString(_("x slot must be numeric \"double\""));
    if (length(fact) > 0 && getAttrib(fact, R_NamesSymbol) == R_NilValue)
    return mkString(_("factors slot must be named list"));
    if (length(rc) > 0 && getAttrib(rc, R_NamesSymbol) == R_NilValue)
    return mkString(_("rcond slot must be named numeric vector"));
    return ScalarLogical(1);
}

static
double get_norm(SEXP obj, char *typstr)
{
    char typnm[] = {'\0', '\0'};
    int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym));
    double *work = (double *) NULL;

    typnm[0] = norm_type(typstr);
    if (*typnm == 'I') {
        work = (double *) R_alloc(dims[0], sizeof(double));
    }
    return F77_CALL(dlange)(typstr, dims, dims+1,
                REAL(GET_SLOT(obj, Matrix_xSym)),
                dims, work);
}

SEXP dgeMatrix_norm(SEXP obj, SEXP type)
{
    return ScalarReal(get_norm(obj, CHAR(asChar(type))));
}

static
double set_rcond(SEXP obj, char *typstr)
{
    char typnm[] = {'\0', '\0'};
    SEXP rcv = GET_SLOT(obj, Matrix_rcondSym);
    double rcond = get_double_by_name(rcv, typnm);

    typnm[0] = rcond_type(typstr);
    if (R_IsNA(rcond)) {
        SEXP LU = dgeMatrix_LU(obj);
    int *dims = INTEGER(GET_SLOT(LU, Matrix_DimSym)), info;
    double anorm = get_norm(obj, typstr);

    if (dims[0] != dims[1] || dims[0] < 1)
            error(_("rcond requires a square, non-empty matrix"));
    F77_CALL(dgecon)(typnm,
             dims, REAL(GET_SLOT(LU, Matrix_xSym)),
             dims, &anorm, &rcond,
             (double *) R_alloc(4*dims[0], sizeof(double)),
             (int *) R_alloc(dims[0], sizeof(int)), &info);
    SET_SLOT(obj, Matrix_rcondSym,
         set_double_by_name(rcv, rcond, typnm));
    }
    return rcond;
}

SEXP dgeMatrix_rcond(SEXP obj, SEXP type)
{
  return ScalarReal(set_rcond(obj, CHAR(asChar(type))));
}

SEXP dgeMatrix_crossprod(SEXP x)
{
    SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dpoMatrix")));
    int *Dims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
    *vDims;
    int i, n = Dims[1];
    int nsqr = n * n;
    double one = 1.0, *xvals, zero = 0.0;

    SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));
    SET_SLOT(val, Matrix_rcondSym, allocVector(REALSXP, 0));
    SET_SLOT(val, Matrix_uploSym, mkString("U"));
    SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));
    vDims = INTEGER(GET_SLOT(val, Matrix_DimSym));
    vDims[0] = vDims[1] = n;
    SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, nsqr));
    xvals = REAL(GET_SLOT(val, Matrix_xSym));
    for (i = 0; i < nsqr; i++) xvals[i] = 0.; /* keep valgrind happy */
    F77_CALL(dsyrk)("U", "T", vDims, Dims,
            &one, REAL(GET_SLOT(x, Matrix_xSym)), Dims,
            &zero, xvals, vDims);
    UNPROTECT(1);
    return val;
}

SEXP dgeMatrix_dgeMatrix_crossprod(SEXP x, SEXP y)
{
    SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
    int *xDims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
    *yDims = INTEGER(GET_SLOT(y, Matrix_DimSym)),
    *vDims;
    int m = xDims[1], n = yDims[1];
    double one = 1.0, zero = 0.0;

    SET_SLOT(val, Matrix_rcondSym, allocVector(REALSXP, 0));
    SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));
    SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));
    vDims = INTEGER(GET_SLOT(val, Matrix_DimSym));
    if ((*xDims) > 0 && (*yDims) > 0 && n > 0 && m > 0) {
    if (*xDims != *yDims)
        error(_("Dimensions of x and y are not compatible for crossprod"));
    vDims[0] = m; vDims[1] = n;
    SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, m * n));
    F77_CALL(dgemm)("T", "N", xDims + 1, yDims + 1, xDims, &one,
            REAL(GET_SLOT(x, Matrix_xSym)), xDims,
            REAL(GET_SLOT(y, Matrix_xSym)), yDims,
            &zero, REAL(GET_SLOT(val, Matrix_xSym)),
            xDims + 1);
    }
    UNPROTECT(1);
    return val;
}

SEXP dgeMatrix_matrix_crossprod(SEXP x, SEXP y)
{
    SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
    int *xDims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
    *yDims = INTEGER(getAttrib(y, R_DimSymbol)),
    *vDims;
    int m = xDims[1], n = yDims[1];
    double one = 1.0, zero = 0.0;

    if (!(isMatrix(y) && isReal(y)))
    error(_("Argument y must be a numeric (real) matrix"));
    SET_SLOT(val, Matrix_rcondSym, allocVector(REALSXP, 0));
    SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));
    SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));
    vDims = INTEGER(GET_SLOT(val, Matrix_DimSym));
    if ((*xDims) > 0 && (*yDims) > 0 && n > 0 && m > 0) {
    if (*xDims != *yDims)
        error(_("Dimensions of x and y are not compatible for crossprod"));
    vDims[0] = m; vDims[1] = n;
    SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, m * n));
    F77_CALL(dgemm)("T", "N", xDims + 1, yDims + 1, xDims, &one,
            REAL(GET_SLOT(x, Matrix_xSym)), xDims,
            REAL(y), yDims,
            &zero, REAL(GET_SLOT(val, Matrix_xSym)),
            xDims + 1);
    }
    UNPROTECT(1);
    return val;
}

SEXP dgeMatrix_getDiag(SEXP x)
{
    int *dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
    int i, m = dims[0], nret = (m < dims[1]) ? m : dims[1];
    SEXP ret = PROTECT(allocVector(REALSXP, nret)),
    xv = GET_SLOT(x, Matrix_xSym);

    for (i = 0; i < nret; i++) {
    REAL(ret)[i] = REAL(xv)[i * (m + 1)];
    }
    UNPROTECT(1);
    return ret;
}

SEXP dgeMatrix_LU(SEXP x)
{
    SEXP val = get_factors(x, "LU");
    int *dims, npiv, info;

    if (val != R_NilValue) return val;
    dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
    if (dims[0] < 1 || dims[1] < 1)
    error(_("Cannot factor a matrix with zero extents"));
    npiv = (dims[0] <dims[1]) ? dims[0] : dims[1];
    val = PROTECT(NEW_OBJECT(MAKE_CLASS("LU")));
    SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(x, Matrix_xSym)));
    SET_SLOT(val, Matrix_DimSym, duplicate(GET_SLOT(x, Matrix_DimSym)));
    SET_SLOT(val, install("pivot"), allocVector(INTSXP, npiv));
    F77_CALL(dgetrf)(dims, dims + 1, REAL(GET_SLOT(val, Matrix_xSym)),
             dims, INTEGER(GET_SLOT(val, install("pivot"))),
             &info);
    if (info) error(_("Lapack routine dgetrf returned error code %d"), info);
    UNPROTECT(1);
    return set_factors(x, val, "LU");
}

SEXP dgeMatrix_determinant(SEXP x, SEXP logarithm)
{
    int lg = asLogical(logarithm);
    SEXP lu = dgeMatrix_LU(x);
    int *dims = INTEGER(GET_SLOT(lu, Matrix_DimSym)),
    *jpvt = INTEGER(GET_SLOT(lu, install("pivot"))),
    i, n = dims[0], sign = 1;
    double *luvals = REAL(GET_SLOT(lu, Matrix_xSym)), modulus;

    if (n != dims[1])
    error(_("Determinant requires a square matrix"));
    for (i = 0; i < n; i++) if (jpvt[i] != (i + 1)) sign = -sign;
    if (lg) {
    modulus = 0.0;
    for (i = 0; i < n; i++) {
        double dii = luvals[i*(n + 1)]; /* ith diagonal element */
        modulus += log(dii < 0 ? -dii : dii);
        if (dii < 0) sign = -sign;
    }
    } else {
    modulus = 1.0;
    for (i = 0; i < n; i++)
        modulus *= luvals[i*(n + 1)];
    if (modulus < 0) {
        modulus = -modulus;
        sign = -sign;
    }
    }
    return as_det_obj(modulus, lg, sign);
}

SEXP dgeMatrix_solve(SEXP a)
{
    SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix"))),
    lu = dgeMatrix_LU(a);
    int *dims = INTEGER(GET_SLOT(lu, Matrix_DimSym)),
    *pivot = INTEGER(GET_SLOT(lu, install("pivot")));
    double *x, tmp;
    int info, lwork = -1;


    if (dims[0] != dims[1]) error(_("Solve requires a square matrix"));
    SET_SLOT(val, Matrix_rcondSym, allocVector(REALSXP, 0));
    SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));
    SET_SLOT(val, Matrix_xSym, duplicate(GET_SLOT(lu, Matrix_xSym)));
    x = REAL(GET_SLOT(val, Matrix_xSym));
    SET_SLOT(val, Matrix_DimSym, duplicate(GET_SLOT(lu, Matrix_DimSym)));
    F77_CALL(dgetri)(dims, x, dims, pivot, &tmp, &lwork, &info);
    lwork = (int) tmp;
    F77_CALL(dgetri)(dims, x, dims, pivot,
             (double *) R_alloc((size_t) lwork, sizeof(double)),
             &lwork, &info);
    UNPROTECT(1);
    return val;
}

SEXP dgeMatrix_dgeMatrix_mm(SEXP a, SEXP b)
{
    int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
    *bdims = INTEGER(GET_SLOT(b, Matrix_DimSym)),
    *cdims,
    m = adims[0], n = bdims[1], k = adims[1];
    SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
    char *trans = "N";
    double one = 1., zero = 0.;

    if (bdims[0] != k)
    error(_("Matrices are not conformable for multiplication"));
    if (m < 1 || n < 1 || k < 1)
    error(_("Matrices with zero extents cannot be multiplied"));
    SET_SLOT(val, Matrix_rcondSym, allocVector(REALSXP, 0));
    SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));
    SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, m * n));
    SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));
    cdims = INTEGER(GET_SLOT(val, Matrix_DimSym));
    cdims[0] = m; cdims[1] = n;
    F77_CALL(dgemm)(trans, trans, adims, bdims+1, bdims, &one,
            REAL(GET_SLOT(a, Matrix_xSym)), adims,
            REAL(GET_SLOT(b, Matrix_xSym)), bdims,
            &zero, REAL(GET_SLOT(val, Matrix_xSym)), adims);
    UNPROTECT(1);
    return val;
}

SEXP dgeMatrix_svd(SEXP x, SEXP nnu, SEXP nnv)
{
    int /* nu = asInteger(nnu),
       nv = asInteger(nnv), */
    *dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
    double *xx = REAL(GET_SLOT(x, Matrix_xSym));
    SEXP val = PROTECT(allocVector(VECSXP, 3));

    if (dims[0] && dims[1]) {
    int m = dims[0], n = dims[1], mm = (m < n)?m:n,
        lwork = -1, info;
    int *iwork = Calloc(8 * mm, int);
    double tmp, *work;
/*  int bdspac = 3*m*m + 4*m, */
/*      wrkbl, maxwrk, minwrk, itmp, */
/*      ione = 1, iminus1 = -1; */
/*  int i1, i2, i3; */

    SET_VECTOR_ELT(val, 0, allocVector(REALSXP, mm));
    SET_VECTOR_ELT(val, 1, allocMatrix(REALSXP, m, mm));
    SET_VECTOR_ELT(val, 2, allocMatrix(REALSXP, mm, n));
    F77_CALL(dgesdd)("S", &m, &n, xx, &m,
             REAL(VECTOR_ELT(val, 0)),
             REAL(VECTOR_ELT(val, 1)), &m,
             REAL(VECTOR_ELT(val, 2)), &mm,
             &tmp, &lwork, iwork, &info);
    lwork = (int) tmp;
/*  F77_CALL(foo)(&i1, &i2, &i3); */
/*  wrkbl = 3*m+(m+n)*i1; */
/*  if (wrkbl < (itmp = 3*m + m*i2)) wrkbl = itmp; */
/*  if (wrkbl < (itmp = 3*m + m*i3)) wrkbl = itmp; */
/*  itmp = bdspac+3*m; */
/*  maxwrk = (wrkbl > itmp) ? wrkbl : itmp; */
/*  minwrk = 3*m + ((bdspac > n) ?  bdspac : n); */
    work = Calloc(lwork, double);
    F77_CALL(dgesdd)("S", &m, &n, xx, &m,
             REAL(VECTOR_ELT(val, 0)),
             REAL(VECTOR_ELT(val, 1)), &m,
             REAL(VECTOR_ELT(val, 2)), &mm,
             work, &lwork, iwork, &info);
    Free(iwork); Free(work);
    }
    UNPROTECT(1);
    return val;
}

static double padec [] =   /*  Constants for matrix exponential calculation. */
{
  5.0000000000000000e-1,
  1.1666666666666667e-1,
  1.6666666666666667e-2,
  1.6025641025641026e-3,
  1.0683760683760684e-4,
  4.8562548562548563e-6,
  1.3875013875013875e-7,
  1.9270852604185938e-9,
};

/**
 * Matrix exponential - based on the code for Octave's expm function.
 *
 * @param x real square matrix to exponentiate
 *
 * @return matrix exponential of x
 */
SEXP dgeMatrix_exp(SEXP x)
{
    SEXP val = PROTECT(duplicate(x));
    int *Dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
    int i, ilo, ilos, ihi, ihis, j, nc = Dims[1], sqpow;
    int ncp1 = Dims[1] + 1, ncsqr = nc * nc;
    int *pivot = Calloc(nc, int);
    int *iperm = Calloc(nc, int);
    double *dpp = Calloc(ncsqr, double), /* denominator power Pade' */
    *npp = Calloc(ncsqr, double), /* numerator power Pade' */
    *perm = Calloc(nc, double),
    *scale = Calloc(nc, double),
    *v = REAL(GET_SLOT(val, Matrix_xSym)),
    *work = Calloc(ncsqr, double), inf_norm, m1_j, /* (-1)^j */
    one = 1., trshift, zero = 0.;

    if (nc < 1 || Dims[0] != nc)
    error(_("Matrix exponential requires square, non-null matrix"));

    /* FIXME: Add special treatment for nc == 1 */

    /* Preconditioning 1.  Shift diagonal by average diagonal if positive. */
    trshift = 0;        /* determine average diagonal element */
    for (i = 0; i < nc; i++) trshift += v[i * ncp1];
    trshift /= nc;
    if (trshift > 0.) {     /* shift diagonal by -trshift */
    for (i = 0; i < nc; i++) v[i * ncp1] -= trshift;
    }

    /* Preconditioning 2. Balancing with dgebal. */
    F77_CALL(dgebal)("P", &nc, v, &nc, &ilo, &ihi, perm, &j);
    if (j) error(_("dgeMatrix_exp: LAPACK routine dgebal returned %d"), j);
    F77_CALL(dgebal)("S", &nc, v, &nc, &ilos, &ihis, scale, &j);
    if (j) error(_("dgeMatrix_exp: LAPACK routine dgebal returned %d"), j);

    /* Preconditioning 3. Scaling according to infinity norm */
    inf_norm = F77_CALL(dlange)("I", &nc, &nc, v, &nc, work);
    sqpow = (inf_norm > 0) ? (int) (1 + log(inf_norm)/log(2.)) : 0;
    if (sqpow < 0) sqpow = 0;
    if (sqpow > 0) {
    double scale_factor = 1.0;
    for (i = 0; i < sqpow; i++) scale_factor *= 2.;
    for (i = 0; i < ncsqr; i++) v[i] /= scale_factor;
    }

    /* Pade' approximation. Powers v^8, v^7, ..., v^1 */
    AZERO(npp, ncsqr);
    AZERO(dpp, ncsqr);
    m1_j = -1;
    for (j = 7; j >=0; j--) {
    double mult = padec[j];
    /* npp = m * npp + padec[j] *m */
    F77_CALL(dgemm)("N", "N", &nc, &nc, &nc, &one, v, &nc, npp, &nc,
            &zero, work, &nc);
    for (i = 0; i < ncsqr; i++) npp[i] = work[i] + mult * v[i];
    /* dpp = m * dpp * (m1_j * padec[j]) * m */
    mult *= m1_j;
    F77_CALL(dgemm)("N", "N", &nc, &nc, &nc, &one, v, &nc, dpp, &nc,
            &zero, work, &nc);
    for (i = 0; i < ncsqr; i++) dpp[i] = work[i] + mult * v[i];
    m1_j *= -1;
    }
    /* Zero power */
    for (i = 0; i < ncsqr; i++) dpp[i] *= -1.;
    for (j = 0; j < nc; j++) {
    npp[j * ncp1] += 1.;
    dpp[j * ncp1] += 1.;
    }

    /* Pade' approximation is solve(dpp, npp) */
    F77_CALL(dgetrf)(&nc, &nc, dpp, &nc, pivot, &j);
    if (j) error(_("dgeMatrix_exp: dgetrf returned error code %d"), j);
    F77_CALL(dgetrs)("N", &nc, &nc, dpp, &nc, pivot, npp, &nc, &j);
    if (j) error(_("dgeMatrix_exp: dgetrs returned error code %d"), j);
    Memcpy(v, npp, ncsqr);

    /* Now undo all of the preconditioning */
    /* Preconditioning 3: square the result for every power of 2 */
    while (sqpow--) {
    F77_CALL(dgemm)("N", "N", &nc, &nc, &nc, &one, v, &nc, v, &nc,
            &zero, work, &nc);
    Memcpy(v, work, ncsqr);
    }
    /* Preconditioning 2: apply inverse scaling */
    for (j = 0; j < nc; j++)
    for (i = 0; i < nc; i++)
        v[i + j * nc] *= scale[i]/scale[j];
    /* Construct balancing permutation vector */
    for (i = 0; i < nc; i++) iperm[i] = i; /* identity permutation */
    /* Leading permutations applied in forward order */
    for (i = 0; i < (ilo - 1); i++) {
    int swapidx = (int) (perm[i]) - 1;
    int tmp = iperm[i];
    iperm[i] = iperm[swapidx];
    iperm[swapidx] = tmp;
    }
    /* Trailing permutations applied in reverse order */
    for (i = nc - 1; i >= ihi; i--) {
    int swapidx = (int) (perm[i]) - 1;
    int tmp = iperm[i];
    iperm[i] = iperm[swapidx];
    iperm[swapidx] = tmp;
    }
    /* Construct inverse balancing permutation vector */
    Memcpy(pivot, iperm, nc);
    for (i = 0; i < nc; i++) iperm[pivot[i]] = i;
    /* Apply inverse permutation */
    Memcpy(work, v, ncsqr);
    for (j = 0; j < nc; j++)
    for (i = 0; i < nc; i++)
        v[i + j * nc] = work[iperm[i] + iperm[j] * nc];

    /* Preconditioning 1: Trace normalization */
    if (trshift > 0.) {
    double mult = exp(trshift);
    for (i = 0; i < ncsqr; i++) v[i] *= mult;
    }

    /* Clean up */
    Free(dpp); Free(npp); Free(perm); Free(iperm); Free(pivot); Free(scale); Free(work);
    UNPROTECT(1);
    return val;
}

SEXP dgeMatrix_Schur(SEXP x, SEXP vectors)
{
    int *dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
    int vecs = asLogical(vectors), info, izero = 0, lwork = -1, n = dims[0];
    double *work, tmp;
    char *nms[] = {"WR", "WI", "T", "Z", ""};
    SEXP val = PROTECT(Matrix_make_named(VECSXP, nms));

    if (n != dims[1] || n < 1)
    error(_("dgeMatrix_Schur: argument x must be a non-null square matrix"));
    SET_VECTOR_ELT(val, 0, allocVector(REALSXP, n));
    SET_VECTOR_ELT(val, 1, allocVector(REALSXP, n));
    SET_VECTOR_ELT(val, 2, allocMatrix(REALSXP, n, n));
    Memcpy(REAL(VECTOR_ELT(val, 2)), REAL(GET_SLOT(x, Matrix_xSym)), n * n);
    SET_VECTOR_ELT(val, 3, allocMatrix(REALSXP, vecs ? n : 0, vecs ? n : 0));
    F77_CALL(dgees)(vecs ? "V" : "N", "N", NULL, dims, (double *) NULL, dims, &izero,
            (double *) NULL, (double *) NULL, (double *) NULL, dims,
            &tmp, &lwork, (int *) NULL, &info);
    if (info) error(_("dgeMatrix_Schur: first call to dgees failed"));
    lwork = (int) tmp;
    work = Calloc(lwork, double);
    F77_CALL(dgees)(vecs ? "V" : "N", "N", NULL, dims, REAL(VECTOR_ELT(val, 2)), dims,
            &izero, REAL(VECTOR_ELT(val, 0)), REAL(VECTOR_ELT(val, 1)),
            REAL(VECTOR_ELT(val, 3)), dims, work, &lwork,
            (int *) NULL, &info);
    if (info) error(_("dgeMatrix_Schur: dgees returned code %d"), info);
    Free(work);
    UNPROTECT(1);
    return val;
}