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 488, Wed Feb 2 19:11:51 2005 UTC revision 492, Thu Feb 3 14:24:03 2005 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("dCholCMatrix")));      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dCholCMatrix")));
32      int *P = (int *) NULL, *Pinv = (int *) NULL;      int *P, *Pinv;
33      double *Ax;      double *Ax;
34    
35        /* FIXME: Check if there is a Cholesky factorization.  If yes,
36           check if the permutation status matches that of the call.  If
37           so, return it. */
38    
39      if (lo) {      if (lo) {
40          x = PROTECT(ssc_transpose(x));          x = PROTECT(ssc_transpose(x));
41          Ai = INTEGER(GET_SLOT(x, Matrix_iSym));          Ai = INTEGER(GET_SLOT(x, Matrix_iSym));
42          Ap = INTEGER(GET_SLOT(x, Matrix_pSym));          Ap = INTEGER(GET_SLOT(x, Matrix_pSym));
43      }      }
44      SET_SLOT(val, Matrix_uploSym, ScalarString(mkChar("L")));      SET_SLOT(val, Matrix_uploSym, mkString("L"));
45      SET_SLOT(val, Matrix_diagSym, ScalarString(mkChar("N")));      SET_SLOT(val, Matrix_diagSym, mkString("U"));
46      SET_SLOT(val, Matrix_DimSym, duplicate(GET_SLOT(x, Matrix_DimSym)));      SET_SLOT(val, Matrix_DimSym, duplicate(GET_SLOT(x, Matrix_DimSym)));
47      SET_SLOT(val, Matrix_ParentSym, allocVector(INTSXP, n));      SET_SLOT(val, Matrix_ParentSym, allocVector(INTSXP, n));
48      Parent = INTEGER(GET_SLOT(val, Matrix_ParentSym));      Parent = INTEGER(GET_SLOT(val, Matrix_ParentSym));
49      SET_SLOT(val, Matrix_pSym, allocVector(INTSXP, n + 1));      SET_SLOT(val, Matrix_pSym, allocVector(INTSXP, n + 1));
50      Lp = INTEGER(GET_SLOT(val, Matrix_pSym));      Lp = INTEGER(GET_SLOT(val, Matrix_pSym));
51      Ax = REAL(GET_SLOT(x, Matrix_xSym));      SET_SLOT(val, Matrix_permSym, allocVector(INTSXP, n));
52        P = INTEGER(GET_SLOT(val, Matrix_permSym));
53      if (piv) {      if (piv) {
54          SEXP trip = PROTECT(dsCMatrix_to_dgTMatrix(x));          SEXP trip = PROTECT(dsCMatrix_to_dgTMatrix(x));
55          SEXP Ti = GET_SLOT(trip, Matrix_iSym);          SEXP Ti = GET_SLOT(trip, Matrix_iSym);
56    
57          /* determine the permutation with Metis */          /* determine the permutation with Metis */
58          Pinv = Calloc(n, int);          Pinv = Calloc(n, int);
         SET_SLOT(val, Matrix_permSym, allocVector(INTSXP, n));  
         P = INTEGER(GET_SLOT(val, Matrix_permSym));  
59          ssc_metis_order(n, Ap, Ai, P, Pinv);          ssc_metis_order(n, Ap, Ai, P, Pinv);
60          /* create a symmetrized form of x */          /* create a symmetrized form of x */
61          nnz = length(Ti);          nnz = length(Ti);
# Line 63  Line 66 
66                         INTEGER(GET_SLOT(trip, Matrix_jSym)),                         INTEGER(GET_SLOT(trip, Matrix_jSym)),
67                         REAL(GET_SLOT(trip, Matrix_xSym)),                         REAL(GET_SLOT(trip, Matrix_xSym)),
68                         Ap, Ai, Ax);                         Ap, Ai, Ax);
69            UNPROTECT(1);
70        } else {
71            int i;
72            for (i = 0; i < n; i++) P[i] = i;
73            Ax = REAL(GET_SLOT(x, Matrix_xSym));
74    
75      }      }
76      R_ldl_symbolic(n, Ap, Ai, Lp, Parent, P, Pinv);      R_ldl_symbolic(n, Ap, Ai, Lp, Parent, (piv) ? P : (int *)NULL,
77                       (piv) ? Pinv : (int *)NULL);
78      nnz = Lp[n];      nnz = Lp[n];
79      SET_SLOT(val, Matrix_iSym, allocVector(INTSXP, nnz));      SET_SLOT(val, Matrix_iSym, allocVector(INTSXP, nnz));
80      SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, nnz));      SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, nnz));
81      SET_SLOT(val, Matrix_DSym, allocVector(REALSXP, n));      SET_SLOT(val, Matrix_DSym, allocVector(REALSXP, n));
82      info = R_ldl_numeric(n, Ap, Ai, Ax,      info = R_ldl_numeric(n, Ap, Ai, Ax, Lp, Parent,
                        Lp, Parent,  
83                         INTEGER(GET_SLOT(val, Matrix_iSym)),                         INTEGER(GET_SLOT(val, Matrix_iSym)),
84                         REAL(GET_SLOT(val, Matrix_xSym)),                         REAL(GET_SLOT(val, Matrix_xSym)),
85                         REAL(GET_SLOT(val, Matrix_DSym)),                         REAL(GET_SLOT(val, Matrix_DSym)),
86                         P, Pinv);                           (piv) ? P : (int *)NULL,
87                             (piv) ? Pinv : (int *)NULL);
88      if (info != n)      if (info != n)
89          error("Leading minor of size %d (possibly after permutation) is indefinite",          error("Leading minor of size %d (possibly after permutation) is indefinite",
90                info + 1);                info + 1);
91      if (piv) {      if (piv) {
         UNPROTECT(1);  
92          Free(Pinv); Free(Ax); Free(Ai); Free(Ap);          Free(Pinv); Free(Ax); Free(Ai); Free(Ap);
93      }      }
94      UNPROTECT(lo ? 2 : 1);      UNPROTECT(lo ? 2 : 1);

Legend:
Removed from v.488  
changed lines
  Added in v.492

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