SCM

SCM Repository

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

Diff of /pkg/src/dsCMatrix.c

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

revision 308, Wed Oct 27 21:39:44 2004 UTC revision 366, Sun Dec 5 13:49:51 2004 UTC
# Line 29  Line 29 
29          n = length(pSlot)-1,          n = length(pSlot)-1,
30          nnz, piv = asLogical(pivot);          nnz, piv = asLogical(pivot);
31      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("sscChol")));      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("sscChol")));
32      int *Flag = Calloc(n, int), *Lnz = Calloc(n, int),      int *P = (int *) NULL, *Pinv = (int *) NULL;
         *P = (int *) NULL, *Pinv = (int *) NULL;  
33      double *Ax;      double *Ax;
34    
35      if (lo) {      if (lo) {
# Line 65  Line 64 
64                         REAL(GET_SLOT(trip, Matrix_xSym)),                         REAL(GET_SLOT(trip, Matrix_xSym)),
65                         Ap, Ai, Ax);                         Ap, Ai, Ax);
66      }      }
67      ldl_symbolic(n, Ap, Ai, Lp, Parent, Lnz, Flag, P, Pinv);      R_ldl_symbolic(n, Ap, Ai, Lp, Parent, P, Pinv);
68      nnz = Lp[n];      nnz = Lp[n];
69      SET_SLOT(val, Matrix_iSym, allocVector(INTSXP, nnz));      SET_SLOT(val, Matrix_iSym, allocVector(INTSXP, nnz));
70      SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, nnz));      SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, nnz));
71      SET_SLOT(val, Matrix_DSym, allocVector(REALSXP, n));      SET_SLOT(val, Matrix_DSym, allocVector(REALSXP, n));
72      info = ldl_numeric(n, Ap, Ai, Ax,      info = R_ldl_numeric(n, Ap, Ai, Ax,
73                         Lp, Parent, Lnz,                         Lp, Parent,
74                         INTEGER(GET_SLOT(val, Matrix_iSym)),                         INTEGER(GET_SLOT(val, Matrix_iSym)),
75                         REAL(GET_SLOT(val, Matrix_xSym)),                         REAL(GET_SLOT(val, Matrix_xSym)),
76                         REAL(GET_SLOT(val, Matrix_DSym)),                         REAL(GET_SLOT(val, Matrix_DSym)),
77                         (double *) R_alloc(n, sizeof(double)), /* Y */                         P, Pinv);
                        (int *) R_alloc(n, sizeof(int)), /* Pattern */  
                        Flag, P, Pinv);  
78      if (info != n)      if (info != n)
79          error("Leading minor of size %d (possibly after permutation) is indefinite",          error("Leading minor of size %d (possibly after permutation) is indefinite",
80                info + 1);                info + 1);
     Free(Flag); Free(Lnz);  
81      if (piv) {      if (piv) {
82          UNPROTECT(1);          UNPROTECT(1);
83          Free(Pinv); Free(Ax); Free(Ai); Free(Ap);          Free(Pinv); Free(Ax); Free(Ai); Free(Ap);
# Line 112  Line 108 
108      Lx = REAL(GET_SLOT(Chol, Matrix_xSym));      Lx = REAL(GET_SLOT(Chol, Matrix_xSym));
109      D = REAL(GET_SLOT(Chol, Matrix_DSym));      D = REAL(GET_SLOT(Chol, Matrix_DSym));
110      for (j = 0; j < ncol; j++, in += n, out += n) {      for (j = 0; j < ncol; j++, in += n, out += n) {
111          if (piv) ldl_perm(n, out, in, INTEGER(perm));          if (piv) R_ldl_perm(n, out, in, INTEGER(perm));
112          else Memcpy(out, in, n);          else Memcpy(out, in, n);
113          ldl_lsolve(n, out, Lp, Li, Lx);          R_ldl_lsolve(n, out, Lp, Li, Lx);
114          ldl_dsolve(n, out, D);          R_ldl_dsolve(n, out, D);
115          ldl_ltsolve(n, out, Lp, Li, Lx);          R_ldl_ltsolve(n, out, Lp, Li, Lx);
116          if (piv) ldl_permt(n, out, Memcpy(tmp, out, n), INTEGER(perm));          if (piv) R_ldl_permt(n, out, Memcpy(tmp, out, n), INTEGER(perm));
117      }      }
118      if (piv) Free(tmp);      if (piv) Free(tmp);
119      UNPROTECT(1);      UNPROTECT(1);
# Line 206  Line 202 
202      int i, n = INTEGER(Dims)[0], nz, nza,      int i, n = INTEGER(Dims)[0], nz, nza,
203          *Ap, *Ai, *Lp, *Li, *Parent,          *Ap, *Ai, *Lp, *Li, *Parent,
204          doperm = asLogical(doPerm),          doperm = asLogical(doPerm),
         *Lnz = (int *) R_alloc(n, sizeof(int)),  
         *Flag = (int *) R_alloc(n, sizeof(int)),  
205          *P = (int *) NULL, *Pinv = (int *) NULL;          *P = (int *) NULL, *Pinv = (int *) NULL;
206    
207    
# Line 221  Line 215 
215      Ap = INTEGER(GET_SLOT(x, Matrix_pSym));      Ap = INTEGER(GET_SLOT(x, Matrix_pSym));
216      Ai = INTEGER(GET_SLOT(x, Matrix_iSym));      Ai = INTEGER(GET_SLOT(x, Matrix_iSym));
217      if (doperm) {      if (doperm) {
218          int *perm;          int *perm, *iperm = Calloc(n, int);
219    
220          SET_VECTOR_ELT(ans, 2, allocVector(INTSXP, n));          SET_VECTOR_ELT(ans, 2, allocVector(INTSXP, n));
221          perm = INTEGER(VECTOR_ELT(ans, 2));          perm = INTEGER(VECTOR_ELT(ans, 2));
222          ssc_metis_order(n, Ap, Ai, perm, Flag);          ssc_metis_order(n, Ap, Ai, perm, iperm);
223          ssc_symbolic_permute(n, 1, Flag, Ap, Ai);          ssc_symbolic_permute(n, 1, iperm, Ap, Ai);
224            Free(iperm);
225      }      }
226      SET_VECTOR_ELT(ans, 0, allocVector(INTSXP, n));      SET_VECTOR_ELT(ans, 0, allocVector(INTSXP, n));
227      Parent = INTEGER(VECTOR_ELT(ans, 0));      Parent = INTEGER(VECTOR_ELT(ans, 0));
# Line 236  Line 232 
232      SET_SLOT(tsc, Matrix_DimSym, Dims);      SET_SLOT(tsc, Matrix_DimSym, Dims);
233      SET_SLOT(tsc, Matrix_pSym, allocVector(INTSXP, n + 1));      SET_SLOT(tsc, Matrix_pSym, allocVector(INTSXP, n + 1));
234      Lp = INTEGER(GET_SLOT(tsc, Matrix_pSym));      Lp = INTEGER(GET_SLOT(tsc, Matrix_pSym));
235      ldl_symbolic(n, Ap, Ai, Lp, Parent, Lnz, Flag, P, Pinv);      R_ldl_symbolic(n, Ap, Ai, Lp, Parent, P, Pinv);
236      nz = Lp[n];      nz = Lp[n];
237      SET_SLOT(tsc, Matrix_iSym, allocVector(INTSXP, nz));      SET_SLOT(tsc, Matrix_iSym, allocVector(INTSXP, nz));
238      Li = INTEGER(GET_SLOT(tsc, Matrix_iSym));      Li = INTEGER(GET_SLOT(tsc, Matrix_iSym));
239      SET_SLOT(tsc, Matrix_xSym, allocVector(REALSXP, nz));      SET_SLOT(tsc, Matrix_xSym, allocVector(REALSXP, nz));
240      for (i = 0; i < nza; i++) REAL(Ax)[i] = 0.00001;      for (i = 0; i < nza; i++) REAL(Ax)[i] = 0.00001;
241      for (i = 0; i < n; i++) REAL(Ax)[Ap[i+1]-1] = 10000.;      for (i = 0; i < n; i++) REAL(Ax)[Ap[i+1]-1] = 10000.;
242      i = ldl_numeric(n, Ap, Ai, REAL(Ax), Lp, Parent, Lnz, Li,      i = R_ldl_numeric(n, Ap, Ai, REAL(Ax), Lp, Parent, Li,
243                      REAL(GET_SLOT(tsc, Matrix_xSym)),                      REAL(GET_SLOT(tsc, Matrix_xSym)),
244                      (double *) R_alloc(n, sizeof(double)), /* D */                      (double *) R_alloc(n, sizeof(double)), /* D */
245                      (double *) R_alloc(n, sizeof(double)), /* Y */                      P, Pinv);
                     (int *) R_alloc(n, sizeof(int)), /* Pattern */  
                     Flag, P, Pinv);  
246      UNPROTECT(2);      UNPROTECT(2);
247      return ans;      return ans;
248  }  }

Legend:
Removed from v.308  
changed lines
  Added in v.366

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