| 2 |
#include "dtCMatrix.h" |
#include "dtCMatrix.h" |
| 3 |
#include "cs_utils.h" |
#include "cs_utils.h" |
| 4 |
|
|
| 5 |
SEXP tsc_validate(SEXP x) |
SEXP dtCMatrix_validate(SEXP x) |
| 6 |
{ |
{ |
| 7 |
return triangularMatrix_validate(x); |
return triangularMatrix_validate(x); |
| 8 |
/* see ./dsCMatrix.c or ./dtpMatrix.c on how to do more testing here */ |
/* see ./dsCMatrix.c or ./dtpMatrix.c on how to do more testing here */ |
| 9 |
} |
} |
| 10 |
|
|
| 11 |
|
#if 0 /* no longer used */ |
| 12 |
SEXP tsc_to_dgTMatrix(SEXP x) |
SEXP tsc_to_dgTMatrix(SEXP x) |
| 13 |
{ |
{ |
| 14 |
SEXP ans; |
SEXP ans; |
| 44 |
} |
} |
| 45 |
return ans; |
return ans; |
| 46 |
} |
} |
| 47 |
|
#endif |
| 48 |
|
|
| 49 |
/** |
/** |
| 50 |
* Derive the column pointer vector for the inverse of L from the parent array |
* Derive the column pointer vector for the inverse of L from the parent array |
| 115 |
return ans; |
return ans; |
| 116 |
} |
} |
| 117 |
|
|
|
#if 0 |
|
|
SEXP dtCMatrix_solve(SEXP a) |
|
|
{ |
|
|
SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix"))); |
|
|
int lo = uplo_P(a)[0] == 'L', unit = diag_P(a)[0] == 'U', |
|
|
n = INTEGER(GET_SLOT(a, Matrix_DimSym))[0], nnz, |
|
|
*ai = INTEGER(GET_SLOT(a,Matrix_iSym)), |
|
|
*ap = INTEGER(GET_SLOT(a, Matrix_pSym)), *bi, |
|
|
*bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, n + 1)); |
|
|
int bnz = 10 * ap[n]; /* initial estimate of nnz in b */ |
|
|
int *ri = Calloc(bnz, int), *ind = Calloc(n, int), j; |
|
|
double *ax = REAL(GET_SLOT(a, Matrix_xSym)), *bx; |
|
|
|
|
|
SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(a, Matrix_DimSym))); |
|
|
SET_SLOT(ans, Matrix_DimNamesSym, |
|
|
duplicate(GET_SLOT(a, Matrix_DimNamesSym))); |
|
|
SET_SLOT(ans, Matrix_uploSym, duplicate(GET_SLOT(a, Matrix_uploSym))); |
|
|
SET_SLOT(ans, Matrix_diagSym, duplicate(GET_SLOT(a, Matrix_diagSym))); |
|
|
|
|
|
if (!(lo && unit)) |
|
|
error("code for non-unit or upper triangular not yet written"); |
|
|
/* Initially bp will contain increasing negative values ending at zero. */ |
|
|
/* Later we add the negative of bp[0] to all values. */ |
|
|
bp[n] = 0; |
|
|
for (j = n - 1; j >= 0; j--) { /* columns in reverse order */ |
|
|
int i, i1 = ap[j], i2 = ap[j + 1], k, nr; |
|
|
if (i1 < i2) AZERO(ind, n); |
|
|
for (i = i1; i < i2; i++) { |
|
|
ind[ai[i]] = 1; |
|
|
for (k = -bp[ai[i] + 1]; k < -bp[ai[i]]; k++) ind[ri[k]] = 1; |
|
|
} |
|
|
for (k = 0, nr = 0; k < n; k++) if (ind[k]) nr++; |
|
|
if ((nr - bp[j + 1]) > bnz) { |
|
|
while (nr > (bnz + bp[j + 1])) bnz *= 2; |
|
|
ri = Realloc(ri, bnz, int); |
|
|
} |
|
|
bp[j] = bp[j + 1] - nr; |
|
|
for (k = 0, i = -bp[j + 1]; k < n; k++) if (ind[k]) ri[i++] = k; |
|
|
} |
|
|
bnz = -bp[0]; |
|
|
bi = INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, bnz)); |
|
|
bx = REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, bnz)); |
|
|
for (j = 0; j < n; j++) { |
|
|
int bpnew = bp[j] + bnz; |
|
|
Memcpy(bi + bpnew, ri - bp[j], bp[j + 1] - bp[j]); |
|
|
bp[j] = bpnew; |
|
|
} |
|
|
/* insert code for calculating the actual values here */ |
|
|
for (j = 0; j < bnz; j++) bx[j] = 1; |
|
|
|
|
|
Free(ind); Free(ri); |
|
|
UNPROTECT(1); |
|
|
return ans; |
|
|
} |
|
|
#else |
|
| 118 |
SEXP dtCMatrix_solve(SEXP a) |
SEXP dtCMatrix_solve(SEXP a) |
| 119 |
{ |
{ |
| 120 |
SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix"))); |
SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix"))); |
| 153 |
UNPROTECT(1); |
UNPROTECT(1); |
| 154 |
return ans; |
return ans; |
| 155 |
} |
} |
|
#endif |
|
| 156 |
|
|
| 157 |
SEXP dtCMatrix_matrix_solve(SEXP a, SEXP b, SEXP classed) |
SEXP dtCMatrix_matrix_solve(SEXP a, SEXP b, SEXP classed) |
| 158 |
{ |
{ |
Parent Directory
| 
Revision Log
| 
Patch
