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 296, Mon Oct 4 17:13:29 2004 UTC pkg/src/dsCMatrix.c revision 587, Wed Mar 2 18:19:15 2005 UTC
# Line 1  Line 1 
1  #include "sscMatrix.h"  #include "dsCMatrix.h"
2    
3  SEXP sscMatrix_validate(SEXP obj)  SEXP dsCMatrix_validate(SEXP obj)
4  {  {
5      SEXP uplo = GET_SLOT(obj, Matrix_uploSym);      SEXP uplo = GET_SLOT(obj, Matrix_uploSym);
6      int *Dim = INTEGER(GET_SLOT(obj, Matrix_DimSym));      int *Dim = INTEGER(GET_SLOT(obj, Matrix_DimSym));
7      char *val;      char *val;
8    
9      if (length(uplo) != 1)      if (length(uplo) != 1)
10          return ScalarString(mkChar("uplo slot must have length 1"));          return mkString(_("uplo slot must have length 1"));
11      val = CHAR(STRING_ELT(uplo, 0));      val = CHAR(STRING_ELT(uplo, 0));
12      if (strlen(val) != 1)      if (strlen(val) != 1)
13          return ScalarString(mkChar("uplo[1] must have string length 1"));          return mkString(_("uplo[1] must have string length 1"));
14      if (toupper(*val) != 'U' && toupper(*val) != 'L')      if (*val != 'U' && *val != 'L')
15          return ScalarString(mkChar("uplo[1] must be \"U\" or \"L\""));          return mkString(_("uplo[1] must be \"U\" or \"L\""));
16      if (Dim[0] != Dim[1])      if (Dim[0] != Dim[1])
17          return ScalarString(mkChar("Symmetric matrix must be square"));          return mkString(_("Symmetric matrix must be square"));
18      csc_check_column_sorting(obj);      csc_check_column_sorting(obj);
19      return ScalarLogical(1);      return ScalarLogical(1);
20  }  }
21    
22  SEXP sscMatrix_chol(SEXP x, SEXP pivot)  SEXP dsCMatrix_chol(SEXP x, SEXP pivot)
23  {  {
24      SEXP pSlot = GET_SLOT(x, Matrix_pSym), xorig = x;      SEXP pSlot = GET_SLOT(x, Matrix_pSym), xorig = x;
25      int *Ai = INTEGER(GET_SLOT(x, Matrix_iSym)),      int *Ai = INTEGER(GET_SLOT(x, Matrix_iSym)),
26          *Ap = INTEGER(pSlot),          *Ap = INTEGER(pSlot),
27          *Lp, *Parent, info,          *Lp, *Parent, info,
28          lo = toupper(CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0]) == 'L',          lo = CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0] == 'L',
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("dCholCMatrix")));
32      int *Flag = Calloc(n, int), *Lnz = Calloc(n, int),      int *P, *Pinv;
         *P = (int *) NULL, *Pinv = (int *) NULL;  
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(sscMatrix_to_triplet(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 64  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      ldl_symbolic(n, Ap, Ai, Lp, Parent, Lnz, Flag, 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 = ldl_numeric(n, Ap, Ai, Ax,      info = R_ldl_numeric(n, Ap, Ai, Ax, Lp, Parent,
                        Lp, Parent, Lnz,  
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                         (double *) R_alloc(n, sizeof(double)), /* Y */                           (piv) ? P : (int *)NULL,
87                         (int *) R_alloc(n, sizeof(int)), /* Pattern */                           (piv) ? Pinv : (int *)NULL);
                        Flag, P, Pinv);  
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);
     Free(Flag); Free(Lnz);  
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);
95      return set_factorization(xorig, val, "Cholesky");      return set_factors(xorig, val, "Cholesky");
96  }  }
97    
98  SEXP sscMatrix_matrix_solve(SEXP a, SEXP b)  SEXP dsCMatrix_matrix_solve(SEXP a, SEXP b)
99  {  {
100      SEXP Chol = get_factorization(a, "Cholesky"), perm,      SEXP Chol = get_factors(a, "Cholesky"), perm,
101          val = PROTECT(duplicate(b));          val = PROTECT(duplicate(b));
102      int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),      int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
103          *bdims = INTEGER(getAttrib(b, R_DimSymbol)),          *bdims = INTEGER(getAttrib(b, R_DimSymbol)),
# Line 100  Line 105 
105      double *Lx, *D, *in = REAL(b), *out = REAL(val), *tmp = (double *) NULL;      double *Lx, *D, *in = REAL(b), *out = REAL(val), *tmp = (double *) NULL;
106    
107      if (!(isReal(b) && isMatrix(b)))      if (!(isReal(b) && isMatrix(b)))
108          error("Argument b must be a numeric matrix");          error(_("Argument b must be a numeric matrix"));
109      if (*adims != *bdims || bdims[1] < 1 || *adims < 1)      if (*adims != *bdims || bdims[1] < 1 || *adims < 1)
110          error("Dimensions of system to be solved are inconsistent");          error(_("Dimensions of system to be solved are inconsistent"));
111      if (Chol == R_NilValue) Chol = sscMatrix_chol(a, ScalarLogical(1));      if (Chol == R_NilValue) Chol = dsCMatrix_chol(a, ScalarLogical(1));
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 111  Line 116 
116      Lp = INTEGER(GET_SLOT(Chol, Matrix_pSym));      Lp = INTEGER(GET_SLOT(Chol, Matrix_pSym));
117      Lx = REAL(GET_SLOT(Chol, Matrix_xSym));      Lx = REAL(GET_SLOT(Chol, Matrix_xSym));
118      D = REAL(GET_SLOT(Chol, Matrix_DSym));      D = REAL(GET_SLOT(Chol, Matrix_DSym));
119      for (j = 0; j < bdims[1]; j++, in += n, out += n) {      for (j = 0; j < ncol; j++, in += n, out += n) {
120          if (piv) ldl_perm(n, out, in, INTEGER(perm));          if (piv) R_ldl_perm(n, out, in, INTEGER(perm));
121          else Memcpy(out, in, n);          else Memcpy(out, in, n);
122          ldl_lsolve(n, out, Lp, Li, Lx);          R_ldl_lsolve(n, out, Lp, Li, Lx);
123          ldl_dsolve(n, out, D);          R_ldl_dsolve(n, out, D);
124          ldl_ltsolve(n, out, Lp, Li, Lx);          R_ldl_ltsolve(n, out, Lp, Li, Lx);
125          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));
126      }      }
127      if (piv) Free(tmp);      if (piv) Free(tmp);
128      UNPROTECT(1);      UNPROTECT(1);
129      return val;      return val;
130  }  }
131    
132  SEXP sscMatrix_inverse_factor(SEXP A)  SEXP dsCMatrix_inverse_factor(SEXP A)
133  {  {
134      return R_NilValue;          /* FIXME: Write this function. */      return R_NilValue;          /* FIXME: Write this function. */
135  }  }
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("sscMatrix"))),  
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));      csc_compTr(xdims[0], xdims[1], nnz,
149      SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nnz));                 INTEGER(GET_SLOT(x, Matrix_pSym)), INTEGER(islot),
     SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nnz));  
     csc_components_transpose(xdims[0], xdims[1], nnz,  
                              INTEGER(GET_SLOT(x, Matrix_pSym)),  
                              INTEGER(islot),  
150                               REAL(GET_SLOT(x, Matrix_xSym)),                               REAL(GET_SLOT(x, Matrix_xSym)),
151                               INTEGER(GET_SLOT(ans, Matrix_pSym)),                 INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, xdims[0] + 1)),
152                               INTEGER(GET_SLOT(ans, Matrix_iSym)),                 INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nnz)),
153                               REAL(GET_SLOT(ans, Matrix_xSym)));                 REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nnz)));
154      UNPROTECT(1);      UNPROTECT(1);
155      return ans;      return ans;
156  }  }
157    
158  SEXP sscMatrix_to_triplet(SEXP x)  SEXP dsCMatrix_to_dgTMatrix(SEXP x)
159  {  {
160      SEXP      SEXP
161          ans = PROTECT(NEW_OBJECT(MAKE_CLASS("tripletMatrix"))),          ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgTMatrix"))),
162          islot = GET_SLOT(x, Matrix_iSym),          islot = GET_SLOT(x, Matrix_iSym),
163          pslot = GET_SLOT(x, Matrix_pSym);          pslot = GET_SLOT(x, Matrix_pSym);
164      int *ai, *aj, *iv = INTEGER(islot),      int *ai, *aj, *iv = INTEGER(islot),
# Line 199  Line 199 
199      return ans;      return ans;
200  }  }
201    
202  SEXP sscMatrix_ldl_symbolic(SEXP x, SEXP doPerm)  SEXP dsCMatrix_ldl_symbolic(SEXP x, SEXP doPerm)
203  {  {
204      SEXP Ax, Dims = GET_SLOT(x, Matrix_DimSym),      SEXP Ax, Dims = GET_SLOT(x, Matrix_DimSym),
205          ans = PROTECT(allocVector(VECSXP, 3)), tsc;          ans = PROTECT(allocVector(VECSXP, 3)), tsc;
206      int i, n = INTEGER(Dims)[0], nz, nza,      int i, n = INTEGER(Dims)[0], nz, nza,
207          *Ap, *Ai, *Lp, *Li, *Parent,          *Ap, *Ai, *Lp, *Li, *Parent,
208          doperm = asLogical(doPerm),          doperm = asLogical(doPerm),
         *Lnz = (int *) R_alloc(n, sizeof(int)),  
         *Flag = (int *) R_alloc(n, sizeof(int)),  
209          *P = (int *) NULL, *Pinv = (int *) NULL;          *P = (int *) NULL, *Pinv = (int *) NULL;
210    
211    
212      if (toupper(CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0]) == 'L') {      if (CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0] == 'L') {
213          x = PROTECT(ssc_transpose(x));          x = PROTECT(ssc_transpose(x));
214      } else {      } else {
215          x = PROTECT(duplicate(x));          x = PROTECT(duplicate(x));
# Line 221  Line 219 
219      Ap = INTEGER(GET_SLOT(x, Matrix_pSym));      Ap = INTEGER(GET_SLOT(x, Matrix_pSym));
220      Ai = INTEGER(GET_SLOT(x, Matrix_iSym));      Ai = INTEGER(GET_SLOT(x, Matrix_iSym));
221      if (doperm) {      if (doperm) {
222          int *perm;          int *perm, *iperm = Calloc(n, int);
223    
224          SET_VECTOR_ELT(ans, 2, allocVector(INTSXP, n));          SET_VECTOR_ELT(ans, 2, allocVector(INTSXP, n));
225          perm = INTEGER(VECTOR_ELT(ans, 2));          perm = INTEGER(VECTOR_ELT(ans, 2));
226          ssc_metis_order(n, Ap, Ai, perm, Flag);          ssc_metis_order(n, Ap, Ai, perm, iperm);
227          ssc_symbolic_permute(n, 1, Flag, Ap, Ai);          ssc_symbolic_permute(n, 1, iperm, Ap, Ai);
228            Free(iperm);
229      }      }
230      SET_VECTOR_ELT(ans, 0, allocVector(INTSXP, n));      SET_VECTOR_ELT(ans, 0, allocVector(INTSXP, n));
231      Parent = INTEGER(VECTOR_ELT(ans, 0));      Parent = INTEGER(VECTOR_ELT(ans, 0));
232      SET_VECTOR_ELT(ans, 1, NEW_OBJECT(MAKE_CLASS("tscMatrix")));      SET_VECTOR_ELT(ans, 1, NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
233      tsc = VECTOR_ELT(ans, 1);      tsc = VECTOR_ELT(ans, 1);
234      SET_SLOT(tsc, Matrix_uploSym, ScalarString(mkChar("L")));      SET_SLOT(tsc, Matrix_uploSym, mkString("L"));
235      SET_SLOT(tsc, Matrix_diagSym, ScalarString(mkChar("U")));      SET_SLOT(tsc, Matrix_diagSym, mkString("U"));
236      SET_SLOT(tsc, Matrix_DimSym, Dims);      SET_SLOT(tsc, Matrix_DimSym, Dims);
237      SET_SLOT(tsc, Matrix_pSym, allocVector(INTSXP, n + 1));      SET_SLOT(tsc, Matrix_pSym, allocVector(INTSXP, n + 1));
238      Lp = INTEGER(GET_SLOT(tsc, Matrix_pSym));      Lp = INTEGER(GET_SLOT(tsc, Matrix_pSym));
239      ldl_symbolic(n, Ap, Ai, Lp, Parent, Lnz, Flag, P, Pinv);      R_ldl_symbolic(n, Ap, Ai, Lp, Parent, P, Pinv);
240      nz = Lp[n];      nz = Lp[n];
241      SET_SLOT(tsc, Matrix_iSym, allocVector(INTSXP, nz));      SET_SLOT(tsc, Matrix_iSym, allocVector(INTSXP, nz));
242      Li = INTEGER(GET_SLOT(tsc, Matrix_iSym));      Li = INTEGER(GET_SLOT(tsc, Matrix_iSym));
243      SET_SLOT(tsc, Matrix_xSym, allocVector(REALSXP, nz));      SET_SLOT(tsc, Matrix_xSym, allocVector(REALSXP, nz));
244      for (i = 0; i < nza; i++) REAL(Ax)[i] = 0.00001;      for (i = 0; i < nza; i++) REAL(Ax)[i] = 0.00001;
245      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.;
246      i = ldl_numeric(n, Ap, Ai, REAL(Ax), Lp, Parent, Lnz, Li,      i = R_ldl_numeric(n, Ap, Ai, REAL(Ax), Lp, Parent, Li,
247                      REAL(GET_SLOT(tsc, Matrix_xSym)),                      REAL(GET_SLOT(tsc, Matrix_xSym)),
248                      (double *) R_alloc(n, sizeof(double)), /* D */                      (double *) R_alloc(n, sizeof(double)), /* D */
249                      (double *) R_alloc(n, sizeof(double)), /* Y */                      P, Pinv);
                     (int *) R_alloc(n, sizeof(int)), /* Pattern */  
                     Flag, P, Pinv);  
250      UNPROTECT(2);      UNPROTECT(2);
251      return ans;      return ans;
252  }  }
253    
254  SEXP sscMatrix_metis_perm(SEXP x)  SEXP dsCMatrix_metis_perm(SEXP x)
255  {  {
256      SEXP pSlot = GET_SLOT(x, Matrix_pSym),      SEXP pSlot = GET_SLOT(x, Matrix_pSym),
257          ans = PROTECT(allocVector(VECSXP, 2));          ans = PROTECT(allocVector(VECSXP, 2));

Legend:
Removed from v.296  
changed lines
  Added in v.587

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