The R Project SVN R-packages

#include "HBMM.h"
#include "iohb.h"
#include "mmio.h"

#if 0
SEXP Matrix_readHarwellBoeing(SEXP filename)
{
    char *fnm = CHAR(asChar(filename)), *Type = Calloc(4, char);
    int M, N, nz, Nrhs;
    SEXP ans = R_NilValue;

    readHB_info(fnm, &M, &N, &nz, &Type, &Nrhs);
    Rprintf("Filename: %s, M=%d, N=%d, nz=%d, Type=\"%s\", Nrhs=%d\n",
        fnm, M, N, nz, Type, Nrhs);

    if (toupper(Type[0]) == 'R') { /* Real (double precision) matrix */
    char upT1 = toupper(Type[1]);
    int *dims;

    if (upT1 == 'S') {
        ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dsCMatrix")));
        SET_SLOT(ans, Matrix_uploSym, mkString("L"));
    }
    if (upT1 == 'R' || upT1 == 'U')
        ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgCMatrix")));
    if (ans == R_NilValue) {Free(Type); return ans;}
    SET_SLOT(ans, Matrix_pSym, allocVector(INTSXP, N + 1));
    SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nz));
    SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nz));
    SET_SLOT(ans, Matrix_DimSym, allocVector(INTSXP, 2));
    dims = INTEGER(GET_SLOT(ans, Matrix_DimSym));
    dims[0] = M; dims[1] = N;
    readHB_mat_double(fnm, INTEGER(GET_SLOT(ans, Matrix_pSym)),
              INTEGER(GET_SLOT(ans, Matrix_iSym)),
              REAL(GET_SLOT(ans, Matrix_xSym)));
    }
    Free(Type);
    UNPROTECT(1);
    return ans;
}

SEXP Matrix_readMatrixMarket(SEXP filename)
{
    FILE *conn = (FILE *) NULL; /* -Wall */
    MM_typecode code;
    int *dims, M, N, i, nz;
    SEXP ans = R_NilValue;

    if (isString(filename)) {
    conn = fopen(CHAR(asChar(filename)), "r");
    if (conn == NULL) {
        error("Unable to open file: %s", CHAR(asChar(filename)));
        return R_NilValue;
    }
    } else {
    error("non-string values not presently accepted");
    }

    if ((nz = mm_read_banner(conn, &code))) {
    fclose(conn);
    error("mm_read_banner returned code %d", nz);
    }
    if (!mm_is_valid(code)) {
    fclose(conn);
    error("Invalid code: %s", mm_typecode_to_str(code));
    }

    if (mm_is_sparse(code)) mm_read_mtx_crd_size(conn, &M, &N, &nz);
    if (mm_is_dense(code)) mm_read_mtx_array_size(conn, &M, &N);

    if (mm_is_real(code)) {
    if (mm_is_sparse(code)) {
        int *ii, *jj;
        double *vv;

        if (mm_is_general(code))
        ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgTMatrix")));
        if (mm_is_symmetric(code)) {
        ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dsTMatrix")));
        SET_SLOT(ans, Matrix_uploSym, mkString("L"));
        }
        if (ans == R_NilValue)
        error("Unrecognized matrix type: %s",
               mm_typecode_to_str(code));
        SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nz));
        ii = INTEGER(GET_SLOT(ans, Matrix_iSym));
        SET_SLOT(ans, Matrix_jSym, allocVector(INTSXP, nz));
        jj = INTEGER(GET_SLOT(ans, Matrix_jSym));
        SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nz));
        vv = REAL(GET_SLOT(ans, Matrix_xSym));
        for (i = 0; i < nz; i++) {
        if (fscanf(conn, "%d %d %lg", &ii[i], &jj[i], &vv[i]) != 3)
            error("Premature end of file");
        ii[i] -= 1; /* change indices to zero-based */
        jj[i] -= 1;
        }
    }
    if (mm_is_dense(code)) {
        int j;
        double *vv;

        if (mm_is_general(code))
        ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
        /* Need to adjust reading loop for symmetric matrices */
/*      if (mm_is_symmetric(code)) { */
/*      ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dsyMatrix"))); */
/*      SET_SLOT(ans, Matrix_uploSym, mkString("L")); */
/*      } */
        if (ans == R_NilValue)
        error("Unrecognized matrix type: %s",
               mm_typecode_to_str(code));
        SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, M * N));
        vv = REAL(GET_SLOT(ans, Matrix_xSym));
        for (j = 0; j < N; j++) {
        for (i = 0; i < M; i++) {
            if (fscanf(conn, "%lg", &vv[i + j * M]) != 1)
            error("Premature end of file");
        }
        }
    }
    } else {
    error("Only real matrices handled at this point");
    }
    SET_SLOT(ans, Matrix_DimSym, allocVector(INTSXP, 2));
    dims = INTEGER(GET_SLOT(ans, Matrix_DimSym));
    dims[0] = M; dims[1] = N;

    fclose(conn);
    UNPROTECT(1);
    return ans;
}
#endif

SEXP Matrix_writeHarwellBoeing(SEXP obj, SEXP file, SEXP typep)
{
    char *type = CHAR(asChar(typep)), Type[4] = "RUA";
    int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym)),
    *ii = (int *) NULL, *pp = (int *) NULL;
    int M = dims[0], N = dims[1], nz = -1;
    double *xx = (double *) NULL;

    if (type[2] == 'C' || type[2] == 'T') {
    SEXP islot = GET_SLOT(obj, Matrix_iSym);
    nz = LENGTH(islot);
    ii = INTEGER(islot);
    if (type[2] == 'T') {   /* create column pointers */
        int *i1 = Calloc(nz, int);
        double *x1 = Calloc(nz, double);

        pp = Calloc(N + 1, int);
        triplet_to_col(M, N, nz, ii,
               INTEGER(GET_SLOT(obj, Matrix_jSym)), xx,
               pp, i1, x1);
        nz = pp[N];
        xx = x1;
        ii = i1;
    } else pp = INTEGER(GET_SLOT(obj, Matrix_pSym));
    } else error("Only types 'C' and 'T' allowed");

    if (type[0] == 'D') {
    xx = REAL(GET_SLOT(obj, Matrix_xSym));
    } else error("Only real matrices allowed");

    if (!isString(file))
    error("non-string values for file not presently accepted");

    if (type[1] == 'S') {
    if (*uplo_P(obj) != 'L')
        error("Symmetric matrices must be stored in lower triangle");
    Type[1] = 'S';
    }

    writeHB_mat_double(CHAR(asChar(file)), M, N, nz, pp, ii, xx, 0,
               (double *)NULL, (double *)NULL, (double *)NULL,
               "", "", Type, (char*)NULL, (char*)NULL,
               (char*)NULL, (char*)NULL, "RUA");

    if (type[2] == 'T') {Free(ii); Free(pp); Free(xx);}
    return R_NilValue;
}

SEXP Matrix_writeMatrixMarket(SEXP obj, SEXP file, SEXP typep)
{
    char *type = CHAR(asChar(typep));
    int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym)),
    *ii = (int *) NULL, *jj = (int *) NULL;
    int M = dims[0], N = dims[1], nz = -1;
    MM_typecode matcode;
    double *xx = (double *) NULL;

    mm_set_matrix(&matcode);
    if (type[2] == 'C' || type[2] == 'T') {
    SEXP islot = GET_SLOT(obj, Matrix_iSym);
    nz = LENGTH(islot);
    ii = INTEGER(islot);
    mm_set_coordinate(&matcode);
    } else error("Only types 'C' and 'T' allowed");

    if (type[0] == 'D') {
    xx = REAL(GET_SLOT(obj, Matrix_xSym));
    mm_set_real(&matcode);
    } else error("Only real matrices allowed");

    if (!isString(file))
    error("non-string values for file not currently allowed");

    if (type[1] == 'S') {
    if (*uplo_P(obj) != 'L')
        error("Symmetric matrices must be stored in lower triangle");
    mm_set_symmetric(&matcode);
    }
    if (type[1] == 'G') mm_set_general(&matcode);

    if (type[2] == 'C')
    jj = expand_cmprPt(N, INTEGER(GET_SLOT(obj, Matrix_pSym)),
               Calloc(nz, int));
    if (type[2] == 'T')
    jj = INTEGER(GET_SLOT(obj, Matrix_jSym));
    if (!jj) error("storage mode must be T or C");

    mm_write_mtx_crd(CHAR(STRING_ELT(file, 0)), M, N, nz, ii, jj, xx,
             matcode);

    if (type[2] == 'C') Free(jj);
    return R_NilValue;

}