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

pkg/src/sscMatrix.c revision 478, Wed Feb 2 14:33:51 2005 UTC pkg/src/dsCMatrix.c revision 725, Tue May 10 14:50:39 2005 UTC
# Line 2  Line 2 
2    
3  SEXP dsCMatrix_validate(SEXP obj)  SEXP dsCMatrix_validate(SEXP obj)
4  {  {
5      SEXP uplo = GET_SLOT(obj, Matrix_uploSym);      SEXP val = check_scalar_string(GET_SLOT(obj, Matrix_uploSym),
6                                       "LU", "uplo");
7      int *Dim = INTEGER(GET_SLOT(obj, Matrix_DimSym));      int *Dim = INTEGER(GET_SLOT(obj, Matrix_DimSym));
     char *val;  
8    
9      if (length(uplo) != 1)      if (isString(val)) return val;
         return ScalarString(mkChar("uplo slot must have length 1"));  
     val = CHAR(STRING_ELT(uplo, 0));  
     if (strlen(val) != 1)  
         return ScalarString(mkChar("uplo[1] must have string length 1"));  
     if (toupper(*val) != 'U' && toupper(*val) != 'L')  
         return ScalarString(mkChar("uplo[1] must be \"U\" or \"L\""));  
10      if (Dim[0] != Dim[1])      if (Dim[0] != Dim[1])
11          return ScalarString(mkChar("Symmetric matrix must be square"));          return mkString(_("Symmetric matrix must be square"));
12      csc_check_column_sorting(obj);      csc_check_column_sorting(obj);
13      return ScalarLogical(1);      return ScalarLogical(1);
14  }  }
# Line 25  Line 19 
19      int *Ai = INTEGER(GET_SLOT(x, Matrix_iSym)),      int *Ai = INTEGER(GET_SLOT(x, Matrix_iSym)),
20          *Ap = INTEGER(pSlot),          *Ap = INTEGER(pSlot),
21          *Lp, *Parent, info,          *Lp, *Parent, info,
22          lo = toupper(CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0]) == 'L',          lo = CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0] == 'L',
23          n = length(pSlot)-1,          n = length(pSlot)-1,
24          nnz, piv = asLogical(pivot);          nnz, piv = asLogical(pivot);
25      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("sscChol")));      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dCholCMatrix")));
26      int *P = (int *) NULL, *Pinv = (int *) NULL;      int *P, *Pinv;
27      double *Ax;      double *Ax;
28    
29        /* FIXME: Check if there is a Cholesky factorization.  If yes,
30           check if the permutation status matches that of the call.  If
31           so, return it. */
32    
33      if (lo) {      if (lo) {
34          x = PROTECT(ssc_transpose(x));          x = PROTECT(ssc_transpose(x));
35          Ai = INTEGER(GET_SLOT(x, Matrix_iSym));          Ai = INTEGER(GET_SLOT(x, Matrix_iSym));
36          Ap = INTEGER(GET_SLOT(x, Matrix_pSym));          Ap = INTEGER(GET_SLOT(x, Matrix_pSym));
37      }      }
38      SET_SLOT(val, Matrix_uploSym, ScalarString(mkChar("L")));      SET_SLOT(val, Matrix_uploSym, mkString("L"));
39      SET_SLOT(val, Matrix_diagSym, ScalarString(mkChar("N")));      SET_SLOT(val, Matrix_diagSym, mkString("U"));
40      SET_SLOT(val, Matrix_DimSym, duplicate(GET_SLOT(x, Matrix_DimSym)));      SET_SLOT(val, Matrix_DimSym, duplicate(GET_SLOT(x, Matrix_DimSym)));
41      SET_SLOT(val, Matrix_ParentSym, allocVector(INTSXP, n));      SET_SLOT(val, Matrix_ParentSym, allocVector(INTSXP, n));
42      Parent = INTEGER(GET_SLOT(val, Matrix_ParentSym));      Parent = INTEGER(GET_SLOT(val, Matrix_ParentSym));
43      SET_SLOT(val, Matrix_pSym, allocVector(INTSXP, n + 1));      SET_SLOT(val, Matrix_pSym, allocVector(INTSXP, n + 1));
44      Lp = INTEGER(GET_SLOT(val, Matrix_pSym));      Lp = INTEGER(GET_SLOT(val, Matrix_pSym));
45      Ax = REAL(GET_SLOT(x, Matrix_xSym));      SET_SLOT(val, Matrix_permSym, allocVector(INTSXP, n));
46        P = INTEGER(GET_SLOT(val, Matrix_permSym));
47      if (piv) {      if (piv) {
48          SEXP trip = PROTECT(dsCMatrix_to_dgTMatrix(x));          SEXP trip = PROTECT(dsCMatrix_to_dgTMatrix(x));
49          SEXP Ti = GET_SLOT(trip, Matrix_iSym);          SEXP Ti = GET_SLOT(trip, Matrix_iSym);
50    
51          /* determine the permutation with Metis */          /* determine the permutation with Metis */
52          Pinv = Calloc(n, int);          Pinv = Calloc(n, int);
         SET_SLOT(val, Matrix_permSym, allocVector(INTSXP, n));  
         P = INTEGER(GET_SLOT(val, Matrix_permSym));  
53          ssc_metis_order(n, Ap, Ai, P, Pinv);          ssc_metis_order(n, Ap, Ai, P, Pinv);
54          /* create a symmetrized form of x */          /* create a symmetrized form of x */
55          nnz = length(Ti);          nnz = length(Ti);
56          Ai = Calloc(nnz, int);          Ai = Calloc(nnz, int);
57          Ax = Calloc(nnz, double);          Ax = Calloc(nnz, double);
58          Ap = Calloc(n + 1, int);          Ap = Calloc(n + 1, int);
59          dgTMatrix_to_dgCMatrix(n, n, nnz, INTEGER(Ti),          triplet_to_col(n, n, nnz, INTEGER(Ti),
60                         INTEGER(GET_SLOT(trip, Matrix_jSym)),                         INTEGER(GET_SLOT(trip, Matrix_jSym)),
61                         REAL(GET_SLOT(trip, Matrix_xSym)),                         REAL(GET_SLOT(trip, Matrix_xSym)),
62                         Ap, Ai, Ax);                         Ap, Ai, Ax);
63            UNPROTECT(1);
64        } else {
65            int i;
66            for (i = 0; i < n; i++) P[i] = i;
67            Ax = REAL(GET_SLOT(x, Matrix_xSym));
68    
69      }      }
70      R_ldl_symbolic(n, Ap, Ai, Lp, Parent, P, Pinv);      R_ldl_symbolic(n, Ap, Ai, Lp, Parent, (piv) ? P : (int *)NULL,
71                       (piv) ? Pinv : (int *)NULL);
72      nnz = Lp[n];      nnz = Lp[n];
73      SET_SLOT(val, Matrix_iSym, allocVector(INTSXP, nnz));      SET_SLOT(val, Matrix_iSym, allocVector(INTSXP, nnz));
74      SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, nnz));      SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, nnz));
75      SET_SLOT(val, Matrix_DSym, allocVector(REALSXP, n));      SET_SLOT(val, Matrix_DSym, allocVector(REALSXP, n));
76      info = R_ldl_numeric(n, Ap, Ai, Ax,      info = R_ldl_numeric(n, Ap, Ai, Ax, Lp, Parent,
                        Lp, Parent,  
77                         INTEGER(GET_SLOT(val, Matrix_iSym)),                         INTEGER(GET_SLOT(val, Matrix_iSym)),
78                         REAL(GET_SLOT(val, Matrix_xSym)),                         REAL(GET_SLOT(val, Matrix_xSym)),
79                         REAL(GET_SLOT(val, Matrix_DSym)),                         REAL(GET_SLOT(val, Matrix_DSym)),
80                         P, Pinv);                           (piv) ? P : (int *)NULL,
81                             (piv) ? Pinv : (int *)NULL);
82      if (info != n)      if (info != n)
83          error("Leading minor of size %d (possibly after permutation) is indefinite",          error(_("Leading minor of size %d (possibly after permutation) is indefinite"),
84                info + 1);                info + 1);
85      if (piv) {      if (piv) {
         UNPROTECT(1);  
86          Free(Pinv); Free(Ax); Free(Ai); Free(Ap);          Free(Pinv); Free(Ax); Free(Ai); Free(Ap);
87      }      }
88      UNPROTECT(lo ? 2 : 1);      UNPROTECT(lo ? 2 : 1);
89      return set_factors(xorig, val, "Cholesky");      return set_factors(xorig, val, "Cholesky");
90  }  }
91    
92  SEXP dsCMatrix_matrix_solve(SEXP a, SEXP b)  SEXP dsCMatrix_matrix_solve(SEXP a, SEXP b, SEXP classed)
93  {  {
94        int cl = asLogical(classed);
95      SEXP Chol = get_factors(a, "Cholesky"), perm,      SEXP Chol = get_factors(a, "Cholesky"), perm,
96          val = PROTECT(duplicate(b));          bdP = cl ? GET_SLOT(b, Matrix_DimSym) : getAttrib(b, R_DimSymbol),
97            val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
98      int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),      int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
99          *bdims = INTEGER(getAttrib(b, R_DimSymbol)),          *bdims = INTEGER(bdP),
100          *Li, *Lp, j, n = adims[1], ncol = bdims[1], piv;          *Li, *Lp, j, piv;
101      double *Lx, *D, *in = REAL(b), *out = REAL(val), *tmp = (double *) NULL;      int n = adims[1], ncol = bdims[1];
102        double *Lx, *D, *in = REAL(cl ? GET_SLOT(b, Matrix_xSym) : b),
103      if (!(isReal(b) && isMatrix(b)))          *out = REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, n * ncol)),
104          error("Argument b must be a numeric matrix");          *tmp = (double *) NULL;
105      if (*adims != *bdims || bdims[1] < 1 || *adims < 1)  
106          error("Dimensions of system to be solved are inconsistent");      if (!cl && !(isReal(b) && isMatrix(b)))
107            error(_("Argument b must be a numeric matrix"));
108        if (*adims != *bdims || ncol < 1 || *adims < 1)
109            error(_("Dimensions of system to be solved are inconsistent"));
110      if (Chol == R_NilValue) Chol = dsCMatrix_chol(a, ScalarLogical(1));      if (Chol == R_NilValue) Chol = dsCMatrix_chol(a, ScalarLogical(1));
111        SET_SLOT(val, Matrix_DimSym, duplicate(bdP));
112      perm = GET_SLOT(Chol, Matrix_permSym);      perm = GET_SLOT(Chol, Matrix_permSym);
113      piv = length(perm);      piv = length(perm);
114      if (piv) tmp = Calloc(n, double);      if (piv) tmp = Calloc(n, double);
# Line 127  Line 136 
136    
137  SEXP ssc_transpose(SEXP x)  SEXP ssc_transpose(SEXP x)
138  {  {
139      SEXP      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dsCMatrix"))),
         ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dsCMatrix"))),  
140          islot = GET_SLOT(x, Matrix_iSym);          islot = GET_SLOT(x, Matrix_iSym);
141      int nnz = length(islot),      int nnz = length(islot), *adims,
         *adims = INTEGER(GET_SLOT(ans, Matrix_DimSym)),  
142          *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym));          *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym));
143    
144        adims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));
145      adims[0] = xdims[1]; adims[1] = xdims[0];      adims[0] = xdims[1]; adims[1] = xdims[0];
146      if (toupper(CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0]) == 'U')      if (CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0] == 'U')
147          SET_SLOT(ans, Matrix_uploSym, ScalarString(mkChar("L")));          SET_SLOT(ans, Matrix_uploSym, mkString("L"));
148      SET_SLOT(ans, Matrix_pSym, allocVector(INTSXP, xdims[0] + 1));      else
149      SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nnz));          SET_SLOT(ans, Matrix_uploSym, mkString("U"));
150      SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nnz));      csc_compTr(xdims[0], xdims[1], nnz,
151      csc_components_transpose(xdims[0], xdims[1], nnz,                 INTEGER(GET_SLOT(x, Matrix_pSym)), INTEGER(islot),
                              INTEGER(GET_SLOT(x, Matrix_pSym)),  
                              INTEGER(islot),  
152                               REAL(GET_SLOT(x, Matrix_xSym)),                               REAL(GET_SLOT(x, Matrix_xSym)),
153                               INTEGER(GET_SLOT(ans, Matrix_pSym)),                 INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, xdims[0] + 1)),
154                               INTEGER(GET_SLOT(ans, Matrix_iSym)),                 INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nnz)),
155                               REAL(GET_SLOT(ans, Matrix_xSym)));                 REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nnz)));
156      UNPROTECT(1);      UNPROTECT(1);
157      return ans;      return ans;
158  }  }
# Line 205  Line 211 
211          *P = (int *) NULL, *Pinv = (int *) NULL;          *P = (int *) NULL, *Pinv = (int *) NULL;
212    
213    
214      if (toupper(CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0]) == 'L') {      if (CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0] == 'L') {
215          x = PROTECT(ssc_transpose(x));          x = PROTECT(ssc_transpose(x));
216      } else {      } else {
217          x = PROTECT(duplicate(x));          x = PROTECT(duplicate(x));
# Line 227  Line 233 
233      Parent = INTEGER(VECTOR_ELT(ans, 0));      Parent = INTEGER(VECTOR_ELT(ans, 0));
234      SET_VECTOR_ELT(ans, 1, NEW_OBJECT(MAKE_CLASS("dtCMatrix")));      SET_VECTOR_ELT(ans, 1, NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
235      tsc = VECTOR_ELT(ans, 1);      tsc = VECTOR_ELT(ans, 1);
236      SET_SLOT(tsc, Matrix_uploSym, ScalarString(mkChar("L")));      SET_SLOT(tsc, Matrix_uploSym, mkString("L"));
237      SET_SLOT(tsc, Matrix_diagSym, ScalarString(mkChar("U")));      SET_SLOT(tsc, Matrix_diagSym, mkString("U"));
238      SET_SLOT(tsc, Matrix_DimSym, Dims);      SET_SLOT(tsc, Matrix_DimSym, Dims);
239      SET_SLOT(tsc, Matrix_pSym, allocVector(INTSXP, n + 1));      SET_SLOT(tsc, Matrix_pSym, allocVector(INTSXP, n + 1));
240      Lp = INTEGER(GET_SLOT(tsc, Matrix_pSym));      Lp = INTEGER(GET_SLOT(tsc, Matrix_pSym));

Legend:
Removed from v.478  
changed lines
  Added in v.725

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