# SCM Repository

[matrix] Diff of /pkg/Matrix/src/dtCMatrix.c
 [matrix] / pkg / Matrix / src / dtCMatrix.c

# Diff of /pkg/Matrix/src/dtCMatrix.c

revision 2208, Sun Jul 13 15:51:45 2008 UTC revision 2210, Wed Jul 16 22:20:22 2008 UTC
# Line 74  Line 74
74      }      }
75  }  }
76

77  #undef RETURN  #undef RETURN
78
/**
* Derive the column pointer vector for the inverse of L from the parent array
*
* @param n length of parent array
* @param countDiag 0 for a unit triangular matrix with implicit diagonal, otherwise 1
* @param pr parent vector describing the elimination tree
* @param ap array of length n+1 to be filled with the column pointers
*
* @return the number of non-zero entries (ap[n])
*/
int parent_inv_ap(int n, int countDiag, const int pr[], int ap[])
{
int *sz = Alloca(n, int), j;
R_CheckStack();

for (j = n - 1; j >= 0; j--) {
int parent = pr[j];
sz[j] = (parent < 0) ?  countDiag : (1 + sz[parent]);
}
ap[0] = 0;
for (j = 0; j < n; j++)
ap[j+1] = ap[j] + sz[j];
return ap[n];
}

/**
* Derive the row index array for the inverse of L from the parent array
*
* @param n length of parent array
* @param countDiag 0 for a unit triangular matrix with implicit diagonal, otherwise 1
* @param pr parent vector describing the elimination tree
* @param ai row index vector of length ap[n]
*/
void parent_inv_ai(int n, int countDiag, const int pr[], int ai[])
{
int j, k, pos = 0;
for (j = 0; j < n; j++) {
if (countDiag) ai[pos++] = j;
for (k = pr[j]; k >= 0; k = pr[k]) ai[pos++] = k;
}
}

SEXP Parent_inverse(SEXP par, SEXP unitdiag)
{
SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
int *ap, *ai, *dims, *pr = INTEGER(par),
countDiag = 1 - asLogical(unitdiag),
j, n = length(par), nnz;
double *ax;

if (!isInteger(par)) error(_("par argument must be an integer vector"));
SET_SLOT(ans, Matrix_pSym, allocVector(INTSXP, n + 1));
ap = INTEGER(GET_SLOT(ans, Matrix_pSym));
nnz = parent_inv_ap(n, countDiag, pr, ap);
SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nnz));
ai = INTEGER(GET_SLOT(ans, Matrix_iSym));
SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nnz));
ax = REAL(GET_SLOT(ans, Matrix_xSym));
for (j = 0; j < nnz; j++) ax[j] = 1.;
SET_SLOT(ans, Matrix_DimSym, allocVector(INTSXP, 2));
dims = INTEGER(GET_SLOT(ans, Matrix_DimSym));
dims[0] = dims[1] = n;
SET_SLOT(ans, Matrix_uploSym, mkString("L"));
SET_SLOT(ans, Matrix_diagSym, (countDiag ? mkString("N") : mkString("U")));
parent_inv_ai(n, countDiag, pr, ai);
UNPROTECT(1);
return ans;
}

79  SEXP dtCMatrix_solve(SEXP a)  SEXP dtCMatrix_solve(SEXP a)
80  {  {
81      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
82      CSP A = AS_CSP(a);      CSP A = AS_CSP(Csparse_diagU2N(a));
83      int *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, (A->n) + 1)),      int *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, (A->n) + 1)),
84          bnz = 10 * A->n,        /* initial estimate of nnz in b */          bnz = 10 * A->n,        /* initial estimate of nnz in b */
85          lo = uplo_P(a)[0] == 'L', top;          lo = uplo_P(a)[0] == 'L', top;
# Line 204  Line 134
134  {  {
135      int cl = asLogical(classed);      int cl = asLogical(classed);
136      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
137      CSP A = AS_CSP(a);      CSP A = AS_CSP(Csparse_diagU2N(a));
138      int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),      int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
139          *bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :          *bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :
140                           getAttrib(b, R_DimSymbol));                           getAttrib(b, R_DimSymbol));
# Line 224  Line 154
154      return ans;      return ans;
155  }  }
156
157  SEXP dtCMatrix_upper_solve(SEXP a)  SEXP dtCMatrix_sparse_solve(SEXP a, SEXP b)
158  {  {
159      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
160      int lo = uplo_P(a)[0] == 'L', unit = diag_P(a)[0] == 'U',      CSP A = AS_CSP(Csparse_diagU2N(a)), B = AS_CSP(Csparse_diagU2N(b));
161          n = INTEGER(GET_SLOT(a, Matrix_DimSym))[0],      int *xp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, (B->n) + 1)),
162          *ai = INTEGER(GET_SLOT(a,Matrix_iSym)),          xnz = 10 * B->p[B->n],  /* initial estimate of nnz in x */
163          *ap = INTEGER(GET_SLOT(a, Matrix_pSym)),          lo = uplo_P(a)[0] == 'L', top;
164          *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, n + 1));      /* These arrays must use Calloc because of possible Realloc */
165      int bnz = 10 * ap[n];         /* initial estimate of nnz in b */      int *ti = Calloc(xnz, int), p, j, nz, pos = 0;
166      int *ti = Alloca(bnz, int), j, nz;      double *tx = Calloc(xnz, double);
167      double *ax = REAL(GET_SLOT(a, Matrix_xSym)), *tx = Alloca(bnz, double),      double  *wrk = Alloca(A->n, double);
168          *tmp = Alloca(n, double);      int *xi = Alloca(2*A->n, int);      /* for cs_reach */
169      R_CheckStack();      R_CheckStack();
170
171      if (lo || (!unit))      if (A->m != A->n || B->n < 1 || A->n < 1 || A->n != B->m)
172          error(_("Code written for unit upper triangular unit matrices"));          error(_("Dimensions of system to be solved are inconsistent"));
173      bp[0] = 0;      slot_dup(ans, b, Matrix_DimSym);
174      for (j = 0; j < n; j++) {      SET_DimNames(ans, b);
175          int i, i1 = ap[j + 1];      xp[0] = 0;
176          AZERO(tmp, n);      for (j = 0; j < B->n; j++) {
177          for (i = ap[j]; i < i1; i++) {          /* (wrk[top:n],xi[top:n]) :=  A^{-1} B[,j] */
178              int ii = ai[i], k;          top = cs_spsolve (A, B, j, xi, wrk, 0, lo);
179              if (unit) tmp[ii] -= ax[i];          nz = A->n - top;
180              for (k = bp[ii]; k < bp[ii + 1]; k++) tmp[ti[k]] -= ax[i] * tx[k];          xp[j + 1] = nz + xp[j];
181            if (xp[j + 1] > xnz) {
182                while (xp[j + 1] > xnz) xnz *= 2;
183                ti = Realloc(ti, xnz, int);
184                tx = Realloc(tx, xnz, double);
185          }          }
186          for (i = 0, nz = 0; i < n; i++) if (tmp[i]) nz++;          if (lo)
187          bp[j + 1] = bp[j] + nz;              for(p = top; p < A->n; p++, pos++) {
188          if (bp[j + 1] > bnz) {                  ti[pos] = xi[p];
189              while (bp[j + 1] > bnz) bnz *= 2;                  tx[pos] = wrk[xi[p]];
ti = Realloc(ti, bnz, int);
tx = Realloc(tx, bnz, double);
190          }          }
191          i1 = bp[j];          else /* upper triagonal */
192          for (i = 0; i < n; i++) if (tmp[i]) {ti[i1] = i; tx[i1] = tmp[i]; i1++;}              for(p = A->n - 1; p >= top; p--, pos++) {
193                    ti[pos] = xi[p];
194                    tx[pos] = wrk[xi[p]];
195      }      }
196      nz = bp[n];      }
197        nz = xp[A->n];
198      Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)), ti, nz);      Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)), ti, nz);
199      Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);      Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);
200      slot_dup(ans, a, Matrix_DimSym);
201      SET_DimNames(ans, a);      Free(ti); Free(tx);
slot_dup(ans, a, Matrix_uploSym);
slot_dup(ans, a, Matrix_diagSym);
202      UNPROTECT(1);      UNPROTECT(1);
203      return ans;      return ans;
204  }  }

Legend:
 Removed from v.2208 changed lines Added in v.2210

 root@r-forge.r-project.org ViewVC Help Powered by ViewVC 1.0.0
Thanks to: