The R Project SVN R-packages

Rev

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

Rev 4543 Rev 4560
Line 50... Line 50...
50 
 *
 50 
 *
51 
 * @return the number of non-zero entries (ap[n])
 51 
 * @return the number of non-zero entries (ap[n])
52 
 */
 52 
 */
53 
int parent_inv_ap(int n, int countDiag, const int pr[], int ap[])
 53 
int parent_inv_ap(int n, int countDiag, const int pr[], int ap[])
54 
{
 54 
{
55 
    int *sz = Calloc(n, int), j;
 55 
    int *sz = Alloca(n, int), j;
- 
 
 56 
    R_CheckStack();
56 
 
 57 
 
57 
    for (j = n - 1; j >= 0; j--) {
 58 
    for (j = n - 1; j >= 0; j--) {
58 
	int parent = pr[j];
 59 
	int parent = pr[j];
59 
	sz[j] = (parent < 0) ?  countDiag : (1 + sz[parent]);
 60 
	sz[j] = (parent < 0) ?  countDiag : (1 + sz[parent]);
60 
    }
 61 
    }
61 
    ap[0] = 0;
 62 
    ap[0] = 0;
62 
    for (j = 0; j < n; j++)
 63 
    for (j = 0; j < n; j++)
63 
	ap[j+1] = ap[j] + sz[j];
 64 
	ap[j+1] = ap[j] + sz[j];
64 
    Free(sz);
 - 
 
65 
    return ap[n];
 65 
    return ap[n];
66 
}
 66 
}
67 
 
 67 
 
68 
/**
 68 
/**
69 
 * Derive the row index array for the inverse of L from the parent array
 69 
 * Derive the row index array for the inverse of L from the parent array
Line 112... Line 112...
112 
SEXP dtCMatrix_solve(SEXP a)
 112 
SEXP dtCMatrix_solve(SEXP a)
113 
{
 113 
{
114 
    SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
 114 
    SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
115 
    CSP A = AS_CSP(a);
 115 
    CSP A = AS_CSP(a);
116 
    int *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, (A->n) + 1)),
 116 
    int *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, (A->n) + 1)),
- 
 
 117 
	bnz = 10 * A->n,	/* initial estimate of nnz in b */
117 
	lo = uplo_P(a)[0] == 'L',
 118 
	lo = uplo_P(a)[0] == 'L', top;
118 
	bnz = 10 * A->n;	/* initial estimate of nnz in b */
 119 
    /* These arrays must use Calloc because of possible Realloc */
119 
    int *ti = Calloc(bnz, int), p, j, nz, pos = 0;
 120 
    int *ti = Calloc(bnz, int), p, j, nz, pos = 0;
120 
    double *tx = Calloc(bnz, double), *wrk = Calloc(A->n, double);
 121 
    double *tx = Calloc(bnz, double);
121 
    cs *u = cs_spalloc(A->n, 1,1,1,0);	/* Sparse unit vector */
 122 
    cs *u = cs_spalloc(A->n, 1,1,1,0);	/* Sparse unit vector */
122 
    int top;				/* top of stack */
 123 
    double  *wrk = Alloca(A->n, double);
123 
    int *xi = Calloc(2*A->n, int);	/* cs_reach uses this workspace */
 124 
    int *xi = Alloca(2*A->n, int);	/* for cs_reach */
124 
    R_CheckStack();
 125 
    R_CheckStack();
125 
 
 126 
 
126 
    SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(a, Matrix_DimSym)));
 127 
    SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(a, Matrix_DimSym)));
127 
    SET_DimNames(ans, a);
 128 
    SET_DimNames(ans, a);
128 
    SET_SLOT(ans, Matrix_uploSym, duplicate(GET_SLOT(a, Matrix_uploSym)));
 129 
    SET_SLOT(ans, Matrix_uploSym, duplicate(GET_SLOT(a, Matrix_uploSym)));
Line 155... Line 156...
155 
    }
 156 
    }
156 
    nz = bp[A->n];
 157 
    nz = bp[A->n];
157 
    Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP,  nz)), ti, nz);
 158 
    Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP,  nz)), ti, nz);
158 
    Memcpy(   REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);
 159 
    Memcpy(   REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);
159 
 
 160 
 
160 
    Free(ti); Free(tx);
 161 
    Free(ti); Free(tx); cs_spfree(u);
161 
    Free(wrk); cs_spfree(u); Free(xi);
 - 
 
162 
    UNPROTECT(1);
 162 
    UNPROTECT(1);
163 
    return ans;
 163 
    return ans;
164 
}
 164 
}
165 
 
 165 
 
166 
SEXP dtCMatrix_matrix_solve(SEXP a, SEXP b, SEXP classed)
 166 
SEXP dtCMatrix_matrix_solve(SEXP a, SEXP b, SEXP classed)
Line 194... Line 194...
194 
	n = INTEGER(GET_SLOT(a, Matrix_DimSym))[0],
 194 
	n = INTEGER(GET_SLOT(a, Matrix_DimSym))[0],
195 
	*ai = INTEGER(GET_SLOT(a,Matrix_iSym)),
 195 
	*ai = INTEGER(GET_SLOT(a,Matrix_iSym)),
196 
	*ap = INTEGER(GET_SLOT(a, Matrix_pSym)),
 196 
	*ap = INTEGER(GET_SLOT(a, Matrix_pSym)),
197 
	*bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, n + 1));
 197 
	*bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, n + 1));
198 
    int bnz = 10 * ap[n];	  /* initial estimate of nnz in b */
 198 
    int bnz = 10 * ap[n];	  /* initial estimate of nnz in b */
199 
    int *ti = Calloc(bnz, int), j, nz;
 199 
    int *ti = Alloca(bnz, int), j, nz;
200 
    double *ax = REAL(GET_SLOT(a, Matrix_xSym)), *tx = Calloc(bnz, double),
 200 
    double *ax = REAL(GET_SLOT(a, Matrix_xSym)), *tx = Alloca(bnz, double),
201 
	*tmp = Calloc(n, double);
 201 
	*tmp = Alloca(n, double);
- 
 
 202 
    R_CheckStack();
202 
 
 203 
 
203 
    if (lo || (!unit))
 204 
    if (lo || (!unit))
204 
	error(_("Code written for unit upper triangular unit matrices"));
 205 
	error(_("Code written for unit upper triangular unit matrices"));
205 
    bp[0] = 0;
 206 
    bp[0] = 0;
206 
    for (j = 0; j < n; j++) {
 207 
    for (j = 0; j < n; j++) {
Line 222... Line 223...
222 
	for (i = 0; i < n; i++) if (tmp[i]) {ti[i1] = i; tx[i1] = tmp[i]; i1++;}
 223 
	for (i = 0; i < n; i++) if (tmp[i]) {ti[i1] = i; tx[i1] = tmp[i]; i1++;}
223 
    }
 224 
    }
224 
    nz = bp[n];
 225 
    nz = bp[n];
225 
    Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)), ti, nz);
 226 
    Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)), ti, nz);
226 
    Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);
 227 
    Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);
227 
    Free(tmp); Free(tx); Free(ti);
 - 
 
228 
    SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(a, Matrix_DimSym)));
 228 
    SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(a, Matrix_DimSym)));
229 
    SET_DimNames(ans, a);
 229 
    SET_DimNames(ans, a);
230 
    SET_SLOT(ans, Matrix_uploSym, duplicate(GET_SLOT(a, Matrix_uploSym)));
 230 
    SET_SLOT(ans, Matrix_uploSym, duplicate(GET_SLOT(a, Matrix_uploSym)));
231 
    SET_SLOT(ans, Matrix_diagSym, duplicate(GET_SLOT(a, Matrix_diagSym)));
 231 
    SET_SLOT(ans, Matrix_diagSym, duplicate(GET_SLOT(a, Matrix_diagSym)));
232 
    UNPROTECT(1);
 232 
    UNPROTECT(1);