The R Project SVN R-packages

Rev

Rev 4529 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 4529 Rev 4543
Line 110... Line 110...
110 
}
 110 
}
111 
 
 111 
 
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 
    cs *A = Matrix_as_cs(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 
	lo = uplo_P(a)[0] == 'L',
 117 
	lo = uplo_P(a)[0] == 'L',
118 
	bnz = 10 * A->n;	/* initial estimate of nnz in b */
 118 
	bnz = 10 * A->n;	/* initial estimate of nnz in b */
119 
    int *ti = Calloc(bnz, int), p, j, nz, pos = 0;
 119 
    int *ti = Calloc(bnz, int), p, j, nz, pos = 0;
120 
    double *tx = Calloc(bnz, double), *wrk = Calloc(A->n, double);
 120 
    double *tx = Calloc(bnz, double), *wrk = Calloc(A->n, double);
121 
    cs *u = cs_spalloc(A->n, 1,1,1,0);	/* Sparse unit vector */
 121 
    cs *u = cs_spalloc(A->n, 1,1,1,0);	/* Sparse unit vector */
122 
    int top;				/* top of stack */
 122 
    int top;				/* top of stack */
123 
    int *xi = Calloc(2*A->n, int);	/* cs_reach uses this workspace */
 123 
    int *xi = Calloc(2*A->n, int);	/* cs_reach uses this workspace */
- 
 
 124 
    R_CheckStack();
124 
 
 125 
 
125 
    SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(a, Matrix_DimSym)));
 126 
    SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(a, Matrix_DimSym)));
126 
    SET_DimNames(ans, a);
 127 
    SET_DimNames(ans, a);
127 
    SET_SLOT(ans, Matrix_uploSym, duplicate(GET_SLOT(a, Matrix_uploSym)));
 128 
    SET_SLOT(ans, Matrix_uploSym, duplicate(GET_SLOT(a, Matrix_uploSym)));
128 
    SET_SLOT(ans, Matrix_diagSym, duplicate(GET_SLOT(a, Matrix_diagSym)));
 129 
    SET_SLOT(ans, Matrix_diagSym, duplicate(GET_SLOT(a, Matrix_diagSym)));
Line 154... Line 155...
154 
    }
 155 
    }
155 
    nz = bp[A->n];
 156 
    nz = bp[A->n];
156 
    Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP,  nz)), ti, nz);
 157 
    Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP,  nz)), ti, nz);
157 
    Memcpy(   REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);
 158 
    Memcpy(   REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);
158 
 
 159 
 
159 
    Free(A); Free(ti); Free(tx);
 160 
    Free(ti); Free(tx);
160 
    Free(wrk); cs_spfree(u); Free(xi);
 161 
    Free(wrk); cs_spfree(u); Free(xi);
161 
    UNPROTECT(1);
 162 
    UNPROTECT(1);
162 
    return ans;
 163 
    return ans;
163 
}
 164 
}
164 
 
 165 
 
165 
SEXP dtCMatrix_matrix_solve(SEXP a, SEXP b, SEXP classed)
 166 
SEXP dtCMatrix_matrix_solve(SEXP a, SEXP b, SEXP classed)
166 
{
 167 
{
167 
    int cl = asLogical(classed);
 168 
    int cl = asLogical(classed);
168 
    SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
 169 
    SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
169 
    cs *A = Matrix_as_cs(a);
 170 
    CSP A = AS_CSP(a);
170 
    int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
 171 
    int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
171 
	*bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :
 172 
	*bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :
172 
			 getAttrib(b, R_DimSymbol));
 173 
			 getAttrib(b, R_DimSymbol));
173 
    int j, n = bdims[0], nrhs = bdims[1], lo = (*uplo_P(a) == 'L');
 174 
    int j, n = bdims[0], nrhs = bdims[1], lo = (*uplo_P(a) == 'L');
174 
    double *bx;
 175 
    double *bx;
- 
 
 176 
    R_CheckStack();
175 
 
 177 
 
176 
    if (*adims != n || nrhs < 1 || *adims < 1 || *adims != adims[1])
 178 
    if (*adims != n || nrhs < 1 || *adims < 1 || *adims != adims[1])
177 
	error(_("Dimensions of system to be solved are inconsistent"));
 179 
	error(_("Dimensions of system to be solved are inconsistent"));
178 
    Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2)), bdims, 2);
 180 
    Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2)), bdims, 2);
179 
    /* FIXME: copy dimnames or Dimnames as well */
 181 
    /* FIXME: copy dimnames or Dimnames as well */
180 
    bx = Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, n * nrhs)),
 182 
    bx = Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, n * nrhs)),
181 
		REAL(cl ? GET_SLOT(b, Matrix_xSym):b), n * nrhs);
 183 
		REAL(cl ? GET_SLOT(b, Matrix_xSym):b), n * nrhs);
182 
    for (j = 0; j < nrhs; j++)
 184 
    for (j = 0; j < nrhs; j++)
183 
	lo ? cs_lsolve(A, bx + n * j) : cs_usolve(A, bx + n * j);
 185 
	lo ? cs_lsolve(A, bx + n * j) : cs_usolve(A, bx + n * j);
184 
    Free(A);
 - 
 
185 
    UNPROTECT(1);
 186 
    UNPROTECT(1);
186 
    return ans;
 187 
    return ans;
187 
}
 188 
}
188 
 
 189 
 
189 
SEXP dtCMatrix_upper_solve(SEXP a)
 190 
SEXP dtCMatrix_upper_solve(SEXP a)