Revision f584db857045d794fb5a587b18f1ca3ed40b8679 authored by Douglas Bates on 15 January 2006, 00:00:00 UTC, committed by Gabor Csardi on 15 January 2006, 00:00:00 UTC
1 parent a2d2fd0
lsCMatrix.c
#include "lsCMatrix.h"
/**
* Check the validity of the slots of an lsCMatrix object
*
* @param x Pointer to an lsCMatrix object
*
* @return an SEXP that is either TRUE or a character string
* describing the way in which the object failed the validity check
*/
SEXP lsCMatrix_validate(SEXP x)
{
SEXP val = symmetricMatrix_validate(x);
if(isString(val))
return(val);
else {
/* FIXME needed? ltC* inherits from lgC* which does this in validate*/
SEXP pslot = GET_SLOT(x, Matrix_pSym),
islot = GET_SLOT(x, Matrix_iSym);
int
ncol = length(pslot) - 1,
*xp = INTEGER(pslot),
*xi = INTEGER(islot);
if (csc_unsorted_columns(ncol, xp, xi))
csc_sort_columns(ncol, xp, xi, (double *) NULL);
return ScalarLogical(1);
}
}
/**
* Transpose an lsCMatrix
*
* @param x Pointer to an lsCMatrix object
*
* @return the transpose of x. It represents the same matrix but is
* stored in the opposite triangle.
*/
SEXP lsCMatrix_trans(SEXP x)
{
SEXP Xi = GET_SLOT(x, Matrix_iSym),
xDim = GET_SLOT(x, Matrix_DimSym),
ans = PROTECT(NEW_OBJECT(MAKE_CLASS("lsCMatrix")));
int n = INTEGER(xDim)[0], nz = length(Xi);
int *xj = expand_cmprPt(n, INTEGER(GET_SLOT(x, Matrix_pSym)),
Calloc(nz, int)),
*xi = Memcpy(Calloc(nz, int), Xi, nz);
int up = uplo_P(x)[0] == 'U';
SET_SLOT(ans, Matrix_DimSym, duplicate(xDim));
SET_SLOT(ans, Matrix_DimNamesSym,
duplicate(GET_SLOT(x, Matrix_DimNamesSym)));
SET_SLOT(ans, Matrix_uploSym, mkString(up ? "L" : "U"));
make_upper_triangular(up ? xj : xi, up ? xi : xj, nz);
triplet_to_col(n, n, nz, xi, xj, (double *) NULL,
INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, n + 1)),
INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)),
(double *) NULL);
Free(xj); Free(xi);
UNPROTECT(1);
return ans;
}
/**
* Create the symbolic Cholesky decomposition of an lsCMatrix object
*
* @param x Pointer to an lsCMatrix object
* @param pivot Pointer to a scalar logical indicating if a
* fill-reducing permutation should be determined
*
* @return an lCholCMatrix object
*/
SEXP lsCMatrix_chol(SEXP x, SEXP pivot)
{
int piv = asLogical(pivot);
SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("lCholCMatrix")));
int j, n = INTEGER(GET_SLOT(x, Matrix_DimSym))[0];
int *Xi = INTEGER(GET_SLOT(x, Matrix_iSym)),
*Xp = INTEGER(GET_SLOT(x, Matrix_pSym)),
*P, *Pinv = (int *) NULL, *Parent, *Lp;
double *D = Calloc(n, double), *Tx, *Xx = Calloc(Xp[n], double);
if (uplo_P(x)[0] != 'U')
error(_("Must have uplo == 'U' in x argument to lsCMatrix_chol"));
SET_SLOT(ans, Matrix_uploSym, mkString("L"));
SET_SLOT(ans, Matrix_diagSym, mkString("U"));
SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(x, Matrix_DimSym)));
SET_SLOT(ans, Matrix_DimNamesSym, duplicate(GET_SLOT(x, Matrix_DimNamesSym)));
P = INTEGER(ALLOC_SLOT(ans, Matrix_permSym, INTSXP, n));
if (piv) {
Pinv = Calloc(n, int);
ssc_metis_order(n, Xp, Xi, P, Pinv);
} else {
int i;
for (i = 0; i < n; i++) P[i] = i;
}
Lp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, n + 1));
Parent = INTEGER(ALLOC_SLOT(ans, Matrix_ParentSym, INTSXP, n));
R_ldl_symbolic(n, Xp, Xi, Lp, Parent,
(piv) ? P : (int *) NULL, (piv) ? Pinv : (int *) NULL);
/* Decompose the identity to get Li */
for (j = 0; j < n; j++) { /* Create an identity from Xp, Xi and Xx */
int ii, ii2 = Xp[j + 1];
for (ii = Xp[j]; ii < ii2; ii++)
Xx[ii] = (Xi[ii] == j) ? 1. : 0.;
}
Tx = Calloc(Lp[n], double);
R_ldl_numeric(n, Xp, Xi, Xx, Lp, Parent,
INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, Lp[n])),
Tx, D, (piv) ? P : (int *) NULL,
(piv) ? Pinv : (int *) NULL);
if (piv) Free(Pinv);
Free(Xx); Free(Tx); Free(D);
UNPROTECT(1);
return ans;
}
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...