SCM

SCM Repository

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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1959, Fri Jul 6 16:40:47 2007 UTC revision 1960, Fri Jul 6 16:54:43 2007 UTC
# Line 52  Line 52 
52   */   */
53  int parent_inv_ap(int n, int countDiag, const int pr[], int ap[])  int parent_inv_ap(int n, int countDiag, const int pr[], int ap[])
54  {  {
55      int *sz = Calloc(n, int), j;      int *sz = Alloca(n, int), j;
56        R_CheckStack();
57    
58      for (j = n - 1; j >= 0; j--) {      for (j = n - 1; j >= 0; j--) {
59          int parent = pr[j];          int parent = pr[j];
# Line 61  Line 62 
62      ap[0] = 0;      ap[0] = 0;
63      for (j = 0; j < n; j++)      for (j = 0; j < n; j++)
64          ap[j+1] = ap[j] + sz[j];          ap[j+1] = ap[j] + sz[j];
     Free(sz);  
65      return ap[n];      return ap[n];
66  }  }
67    
# Line 112  Line 112 
112  SEXP dtCMatrix_solve(SEXP a)  SEXP dtCMatrix_solve(SEXP a)
113  {  {
114      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
115      cs *A = Matrix_as_cs(a);      CSP A = AS_CSP(a);
116      int *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, (A->n) + 1)),      int *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, (A->n) + 1)),
117          lo = uplo_P(a)[0] == 'L',          bnz = 10 * A->n,        /* initial estimate of nnz in b */
118          bnz = 10 * A->n;        /* initial estimate of nnz in b */          lo = uplo_P(a)[0] == 'L', top;
119        /* These arrays must use Calloc because of possible Realloc */
120      int *ti = Calloc(bnz, int), p, j, nz, pos = 0;      int *ti = Calloc(bnz, int), p, j, nz, pos = 0;
121      double *tx = Calloc(bnz, double), *wrk = Calloc(A->n, double);      double *tx = Calloc(bnz, double);
122      cs *u = cs_spalloc(A->n, 1,1,1,0);  /* Sparse unit vector */      cs *u = cs_spalloc(A->n, 1,1,1,0);  /* Sparse unit vector */
123      int top;                            /* top of stack */      double  *wrk = Alloca(A->n, double);
124      int *xi = Calloc(2*A->n, int);      /* cs_reach uses this workspace */      int *xi = Alloca(2*A->n, int);      /* for cs_reach */
125        R_CheckStack();
126    
127      SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(a, Matrix_DimSym)));      SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(a, Matrix_DimSym)));
128      SET_DimNames(ans, a);      SET_DimNames(ans, a);
# Line 156  Line 158 
158      Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP,  nz)), ti, nz);      Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP,  nz)), ti, nz);
159      Memcpy(   REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);      Memcpy(   REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);
160    
161      Free(A); Free(ti); Free(tx);      Free(ti); Free(tx); cs_spfree(u);
     Free(wrk); cs_spfree(u); Free(xi);  
162      UNPROTECT(1);      UNPROTECT(1);
163      return ans;      return ans;
164  }  }
# Line 166  Line 167 
167  {  {
168      int cl = asLogical(classed);      int cl = asLogical(classed);
169      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
170      cs *A = Matrix_as_cs(a);      CSP A = AS_CSP(a);
171      int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),      int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
172          *bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :          *bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :
173                           getAttrib(b, R_DimSymbol));                           getAttrib(b, R_DimSymbol));
174      int j, n = bdims[0], nrhs = bdims[1], lo = (*uplo_P(a) == 'L');      int j, n = bdims[0], nrhs = bdims[1], lo = (*uplo_P(a) == 'L');
175      double *bx;      double *bx;
176        R_CheckStack();
177    
178      if (*adims != n || nrhs < 1 || *adims < 1 || *adims != adims[1])      if (*adims != n || nrhs < 1 || *adims < 1 || *adims != adims[1])
179          error(_("Dimensions of system to be solved are inconsistent"));          error(_("Dimensions of system to be solved are inconsistent"));
# Line 181  Line 183 
183                  REAL(cl ? GET_SLOT(b, Matrix_xSym):b), n * nrhs);                  REAL(cl ? GET_SLOT(b, Matrix_xSym):b), n * nrhs);
184      for (j = 0; j < nrhs; j++)      for (j = 0; j < nrhs; j++)
185          lo ? cs_lsolve(A, bx + n * j) : cs_usolve(A, bx + n * j);          lo ? cs_lsolve(A, bx + n * j) : cs_usolve(A, bx + n * j);
     Free(A);  
186      UNPROTECT(1);      UNPROTECT(1);
187      return ans;      return ans;
188  }  }
# Line 195  Line 196 
196          *ap = INTEGER(GET_SLOT(a, Matrix_pSym)),          *ap = INTEGER(GET_SLOT(a, Matrix_pSym)),
197          *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, n + 1));          *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, n + 1));
198      int bnz = 10 * ap[n];         /* initial estimate of nnz in b */      int bnz = 10 * ap[n];         /* initial estimate of nnz in b */
199      int *ti = Calloc(bnz, int), j, nz;      int *ti = Alloca(bnz, int), j, nz;
200      double *ax = REAL(GET_SLOT(a, Matrix_xSym)), *tx = Calloc(bnz, double),      double *ax = REAL(GET_SLOT(a, Matrix_xSym)), *tx = Alloca(bnz, double),
201          *tmp = Calloc(n, double);          *tmp = Alloca(n, double);
202        R_CheckStack();
203    
204      if (lo || (!unit))      if (lo || (!unit))
205          error(_("Code written for unit upper triangular unit matrices"));          error(_("Code written for unit upper triangular unit matrices"));
# Line 223  Line 225 
225      nz = bp[n];      nz = bp[n];
226      Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)), ti, nz);      Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)), ti, nz);
227      Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);      Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);
     Free(tmp); Free(tx); Free(ti);  
228      SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(a, Matrix_DimSym)));      SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(a, Matrix_DimSym)));
229      SET_DimNames(ans, a);      SET_DimNames(ans, a);
230      SET_SLOT(ans, Matrix_uploSym, duplicate(GET_SLOT(a, Matrix_uploSym)));      SET_SLOT(ans, Matrix_uploSym, duplicate(GET_SLOT(a, Matrix_uploSym)));

Legend:
Removed from v.1959  
changed lines
  Added in v.1960

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge