The R Project SVN R-packages

Rev

Rev 4565 | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 4565 Rev 4573
Line 1... Line 1...
1 
				/* Sparse triangular numeric matrices */
 1 
				/* Sparse triangular numeric matrices */
2 
#include "dtCMatrix.h"
 2 
#include "dtCMatrix.h"
3 
#include "cs_utils.h"
 3 
#include "cs_utils.h"
4 
 
 4 
 
- 
 
 5 
#define RETURN(_CH_)   UNPROTECT(1); return (_CH_);
- 
 
 6 
 
5 
/* This should be use for *BOTH* triangular and symmetric Csparse: */
 7 
/* This is used for *BOTH* triangular and symmetric Csparse: */
6 
SEXP tCMatrix_validate(SEXP x)
 8 
SEXP tCMatrix_validate(SEXP x)
7 
{
 9 
{
8 
    SEXP val = xCMatrix_validate(x);/* checks x slot */
 10 
    SEXP val = xCMatrix_validate(x);/* checks x slot */
9 
    if(isString(val))
 11 
    if(isString(val))
10 
	return(val);
 12 
	return(val);
Line 15... Line 17...
15 
	int uploT = (*uplo_P(x) == 'U'),
 17 
	int uploT = (*uplo_P(x) == 'U'),
16 
	    k, nnz = length(islot),
 18 
	    k, nnz = length(islot),
17 
	    *xi = INTEGER(islot),
 19 
	    *xi = INTEGER(islot),
18 
	    *xj = INTEGER(PROTECT(allocVector(INTSXP, nnz)));
 20 
	    *xj = INTEGER(PROTECT(allocVector(INTSXP, nnz)));
19 
 
 21 
 
20 
#define RETURN(_CH_)   UNPROTECT(1); return (_CH_);
 - 
 
21 
 
 - 
 
22 
	expand_cmprPt(length(pslot) - 1, INTEGER(pslot), xj);
 22 
	expand_cmprPt(length(pslot) - 1, INTEGER(pslot), xj);
23 
 
 23 
 
24 
	/* Maybe FIXME: ">" should be ">="	for diag = 'U' (uplo = 'U') */
 24 
	/* Maybe FIXME: ">" should be ">="	for diag = 'U' (uplo = 'U') */
25 
	if(uploT) {
 25 
	if(uploT) {
26 
	    for (k = 0; k < nnz; k++)
 26 
	    for (k = 0; k < nnz; k++)
Line 36... Line 36...
36 
	}
 36 
	}
37 
 
 37 
 
38 
	RETURN(ScalarLogical(1));
 38 
	RETURN(ScalarLogical(1));
39 
    }
 39 
    }
40 
}
 40 
}
- 
 
 41 
 
- 
 
 42 
/* This is used for *BOTH* triangular and symmetric Rsparse: */
- 
 
 43 
SEXP tRMatrix_validate(SEXP x)
- 
 
 44 
{
- 
 
 45 
    SEXP val = xRMatrix_validate(x);/* checks x slot */
- 
 
 46 
    if(isString(val))
- 
 
 47 
	return(val);
- 
 
 48 
    else {
- 
 
 49 
	SEXP
- 
 
 50 
	    jslot = GET_SLOT(x, Matrix_jSym),
- 
 
 51 
	    pslot = GET_SLOT(x, Matrix_pSym);
- 
 
 52 
	int uploT = (*uplo_P(x) == 'U'),
- 
 
 53 
	    k, nnz = length(jslot),
- 
 
 54 
	    *xj = INTEGER(jslot),
- 
 
 55 
	    *xi = INTEGER(PROTECT(allocVector(INTSXP, nnz)));
- 
 
 56 
 
- 
 
 57 
	expand_cmprPt(length(pslot) - 1, INTEGER(pslot), xi);
- 
 
 58 
 
- 
 
 59 
	/* Maybe FIXME: ">" should be ">="	for diag = 'U' (uplo = 'U') */
- 
 
 60 
	if(uploT) {
- 
 
 61 
	    for (k = 0; k < nnz; k++)
- 
 
 62 
		if(xi[k] > xj[k]) {
- 
 
 63 
		    RETURN(mkString(_("uplo='U' must not have sparse entries in lower diagonal")));
- 
 
 64 
		}
- 
 
 65 
	}
- 
 
 66 
	else {
- 
 
 67 
	    for (k = 0; k < nnz; k++)
- 
 
 68 
		if(xi[k] < xj[k]) {
- 
 
 69 
		    RETURN(mkString(_("uplo='L' must not have sparse entries in upper diagonal")));
- 
 
 70 
		}
- 
 
 71 
	}
- 
 
 72 
 
- 
 
 73 
	RETURN(ScalarLogical(1));
- 
 
 74 
    }
- 
 
 75 
}
- 
 
 76 
 
- 
 
 77 
 
41 
#undef RETURN
 78 
#undef RETURN
42 
 
 79 
 
43 
/**
 80 
/**
44 
 * Derive the column pointer vector for the inverse of L from the parent array
 81 
 * Derive the column pointer vector for the inverse of L from the parent array
45 
 *
 82 
 *
Line 122... Line 159...
122 
    cs *u = cs_spalloc(A->n, 1,1,1,0);	/* Sparse unit vector */
 159 
    cs *u = cs_spalloc(A->n, 1,1,1,0);	/* Sparse unit vector */
123 
    double  *wrk = Alloca(A->n, double);
 160 
    double  *wrk = Alloca(A->n, double);
124 
    int *xi = Alloca(2*A->n, int);	/* for cs_reach */
 161 
    int *xi = Alloca(2*A->n, int);	/* for cs_reach */
125 
    R_CheckStack();
 162 
    R_CheckStack();
126 
 
 163 
 
127 
    SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(a, Matrix_DimSym)));
 164 
    slot_dup(ans, a, Matrix_DimSym);
128 
    SET_DimNames(ans, a);
 165 
    SET_DimNames(ans, a);
129 
    SET_SLOT(ans, Matrix_uploSym, duplicate(GET_SLOT(a, Matrix_uploSym)));
 166 
    slot_dup(ans, a, Matrix_uploSym);
130 
    SET_SLOT(ans, Matrix_diagSym, duplicate(GET_SLOT(a, Matrix_diagSym)));
 167 
    slot_dup(ans, a, Matrix_diagSym);
131 
    /* initialize the "constant part" of the sparse unit vector */
 168 
    /* initialize the "constant part" of the sparse unit vector */
132 
    u->x[0] = 1.;
 169 
    u->x[0] = 1.;
133 
    u->p[0] = 0; u->p[1] = 1;
 170 
    u->p[0] = 0; u->p[1] = 1;
134 
    bp[0] = 0;
 171 
    bp[0] = 0;
135 
    for (j = 0; j < A->n; j++) {
 172 
    for (j = 0; j < A->n; j++) {
Line 223... Line 260...
223 
	for (i = 0; i < n; i++) if (tmp[i]) {ti[i1] = i; tx[i1] = tmp[i]; i1++;}
 260 
	for (i = 0; i < n; i++) if (tmp[i]) {ti[i1] = i; tx[i1] = tmp[i]; i1++;}
224 
    }
 261 
    }
225 
    nz = bp[n];
 262 
    nz = bp[n];
226 
    Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)), ti, nz);
 263 
    Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)), ti, nz);
227 
    Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);
 264 
    Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);
228 
    SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(a, Matrix_DimSym)));
 265 
    slot_dup(ans, a, Matrix_DimSym);
229 
    SET_DimNames(ans, a);
 266 
    SET_DimNames(ans, a);
230 
    SET_SLOT(ans, Matrix_uploSym, duplicate(GET_SLOT(a, Matrix_uploSym)));
 267 
    slot_dup(ans, a, Matrix_uploSym);
231 
    SET_SLOT(ans, Matrix_diagSym, duplicate(GET_SLOT(a, Matrix_diagSym)));
 268 
    slot_dup(ans, a, Matrix_diagSym);
232 
    UNPROTECT(1);
 269 
    UNPROTECT(1);
233 
    return ans;
 270 
    return ans;
234 
}
 271 
}