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 70, Mon Apr 12 12:10:01 2004 UTC revision 478, Wed Feb 2 14:33:51 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));
# Line 19  Line 19 
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;
25        int *Ai = INTEGER(GET_SLOT(x, Matrix_iSym)),
26            *Ap = INTEGER(pSlot),
27            *Lp, *Parent, info,
28            lo = toupper(CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0]) == 'L',
29            n = length(pSlot)-1,
30            nnz, piv = asLogical(pivot);
31      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("sscChol")));      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("sscChol")));
32      taucs_ccs_matrix *tm =      int *P = (int *) NULL, *Pinv = (int *) NULL;
33          csc_taucs_ptr(x, TAUCS_DOUBLE | TAUCS_LOWER | TAUCS_SYMMETRIC);      double *Ax;
     int nnz, piv = asLogical(pivot);  
     int *dims;  
     void *L;  
34    
35        if (lo) {
36            x = PROTECT(ssc_transpose(x));
37            Ai = INTEGER(GET_SLOT(x, Matrix_iSym));
38            Ap = INTEGER(GET_SLOT(x, Matrix_pSym));
39        }
40        SET_SLOT(val, Matrix_uploSym, ScalarString(mkChar("L")));
41        SET_SLOT(val, Matrix_diagSym, ScalarString(mkChar("N")));
42        SET_SLOT(val, Matrix_DimSym, duplicate(GET_SLOT(x, Matrix_DimSym)));
43        SET_SLOT(val, Matrix_ParentSym, allocVector(INTSXP, n));
44        Parent = INTEGER(GET_SLOT(val, Matrix_ParentSym));
45        SET_SLOT(val, Matrix_pSym, allocVector(INTSXP, n + 1));
46        Lp = INTEGER(GET_SLOT(val, Matrix_pSym));
47        Ax = REAL(GET_SLOT(x, Matrix_xSym));
48      if (piv) {      if (piv) {
49          int *iperm, *perm;          SEXP trip = PROTECT(dsCMatrix_to_dgTMatrix(x));
50            SEXP Ti = GET_SLOT(trip, Matrix_iSym);
51    
52          SET_SLOT(val, Matrix_permSym, allocVector(INTSXP, tm->n));          /* determine the permutation with Metis */
53          SET_SLOT(val, Matrix_ipermSym, allocVector(INTSXP, tm->n));          Pinv = Calloc(n, int);
54          perm = INTEGER(GET_SLOT(val, Matrix_permSym));          SET_SLOT(val, Matrix_permSym, allocVector(INTSXP, n));
55          iperm = INTEGER(GET_SLOT(val, Matrix_ipermSym));          P = INTEGER(GET_SLOT(val, Matrix_permSym));
56          ssc_metis_order(tm->n, tm->colptr, tm->rowind, perm, iperm);          ssc_metis_order(n, Ap, Ai, P, Pinv);
57          tm = taucs_dccs_permute_symmetrically(tm, perm, iperm);          /* create a symmetrized form of x */
58      }          nnz = length(Ti);
59      if (!(L = taucs_ccs_factor_llt_mf(tm)))          Ai = Calloc(nnz, int);
60          error("Matrix is not positive definite");          Ax = Calloc(nnz, double);
61      if (piv) taucs_dccs_free(tm);          Ap = Calloc(n + 1, int);
62      tm = taucs_supernodal_factor_to_ccs(L);          dgTMatrix_to_dgCMatrix(n, n, nnz, INTEGER(Ti),
63      taucs_supernodal_factor_free(L);                         INTEGER(GET_SLOT(trip, Matrix_jSym)),
64      nnz = tm->colptr[tm->n];                         REAL(GET_SLOT(trip, Matrix_xSym)),
65      SET_SLOT(val, Matrix_pSym, allocVector(INTSXP, tm->n + 1));                         Ap, Ai, Ax);
66      Memcpy(INTEGER(GET_SLOT(val, Matrix_pSym)), tm->colptr, tm->n + 1);      }
67        R_ldl_symbolic(n, Ap, Ai, Lp, Parent, P, Pinv);
68        nnz = Lp[n];
69      SET_SLOT(val, Matrix_iSym, allocVector(INTSXP, nnz));      SET_SLOT(val, Matrix_iSym, allocVector(INTSXP, nnz));
     Memcpy(INTEGER(GET_SLOT(val, Matrix_iSym)), tm->rowind, nnz);  
70      SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, nnz));      SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, nnz));
71      Memcpy(REAL(GET_SLOT(val, Matrix_xSym)), tm->values.d, nnz);      SET_SLOT(val, Matrix_DSym, allocVector(REALSXP, n));
72      dims = INTEGER(GET_SLOT(val, Matrix_DimSym));      info = R_ldl_numeric(n, Ap, Ai, Ax,
73      dims[0] = dims[1] = tm->n;                         Lp, Parent,
74      taucs_dccs_free(tm);                         INTEGER(GET_SLOT(val, Matrix_iSym)),
75                           REAL(GET_SLOT(val, Matrix_xSym)),
76                           REAL(GET_SLOT(val, Matrix_DSym)),
77                           P, Pinv);
78        if (info != n)
79            error("Leading minor of size %d (possibly after permutation) is indefinite",
80                  info + 1);
81        if (piv) {
82      UNPROTECT(1);      UNPROTECT(1);
83      return set_factorization(x, val, "Cholesky");          Free(Pinv); Free(Ax); Free(Ai); Free(Ap);
84        }
85        UNPROTECT(lo ? 2 : 1);
86        return set_factors(xorig, val, "Cholesky");
87  }  }
88    
89  SEXP sscMatrix_matrix_solve(SEXP a, SEXP b)  SEXP dsCMatrix_matrix_solve(SEXP a, SEXP b)
90  {  {
91      SEXP Chol = get_factorization(a, "Cholesky"),      SEXP Chol = get_factors(a, "Cholesky"), perm,
92          val = PROTECT(duplicate(b));          val = PROTECT(duplicate(b));
93      int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),      int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
94          *bdims = INTEGER(getAttrib(b, R_DimSymbol)),          *bdims = INTEGER(getAttrib(b, R_DimSymbol)),
95          j, n = adims[1];          *Li, *Lp, j, n = adims[1], ncol = bdims[1], piv;
96      taucs_ccs_matrix* tm;      double *Lx, *D, *in = REAL(b), *out = REAL(val), *tmp = (double *) NULL;
     double *in = REAL(b), *out = REAL(val) ;  
97    
98      if (!(isReal(b) && isMatrix(b)))      if (!(isReal(b) && isMatrix(b)))
99          error("Argument b must be a numeric matrix");          error("Argument b must be a numeric matrix");
100      if (*adims != *bdims || bdims[1] < 1 || *adims < 1)      if (*adims != *bdims || bdims[1] < 1 || *adims < 1)
101          error("Dimensions of system to be solved are inconsistent");          error("Dimensions of system to be solved are inconsistent");
102      if (Chol == R_NilValue) Chol = sscMatrix_chol(a, ScalarLogical(1.));      if (Chol == R_NilValue) Chol = dsCMatrix_chol(a, ScalarLogical(1));
103      tm = csc_taucs_ptr(Chol, TAUCS_DOUBLE | TAUCS_LOWER | TAUCS_TRIANGULAR);      perm = GET_SLOT(Chol, Matrix_permSym);
104      if (!length(GET_SLOT(Chol, Matrix_permSym))) {      piv = length(perm);
105          for (j = 0; j < bdims[1]; j++, in += n, out += n) {      if (piv) tmp = Calloc(n, double);
106              int errcode = taucs_dccs_solve_llt(tm, out, in);      Li = INTEGER(GET_SLOT(Chol, Matrix_iSym));
107              if (errcode)      Lp = INTEGER(GET_SLOT(Chol, Matrix_pSym));
108                  error("taucs_solve returned error code %d for column %d",      Lx = REAL(GET_SLOT(Chol, Matrix_xSym));
109                        errcode, j + 1);      D = REAL(GET_SLOT(Chol, Matrix_DSym));
110          }      for (j = 0; j < ncol; j++, in += n, out += n) {
111      } else {          if (piv) R_ldl_perm(n, out, in, INTEGER(perm));
112          int *iperm = INTEGER(GET_SLOT(Chol, Matrix_ipermSym));          else Memcpy(out, in, n);
113          double *tmpIn = (double *) R_alloc(n, sizeof(double)),          R_ldl_lsolve(n, out, Lp, Li, Lx);
114              *tmpOut = (double *) R_alloc(n, sizeof(double));          R_ldl_dsolve(n, out, D);
115            R_ldl_ltsolve(n, out, Lp, Li, Lx);
116          for (j = 0; j < bdims[1]; j++, in += n, out += n) {          if (piv) R_ldl_permt(n, out, Memcpy(tmp, out, n), INTEGER(perm));
             int errcode, i;  
                                 /* permute y */  
             for (i = 0; i < n; i++) tmpIn[iperm[i]] = in[i];  
                                 /* solve */  
             errcode = taucs_dccs_solve_llt(tm, tmpOut, tmpIn);  
             if (errcode)  
                 error("taucs_solve returned error code %d for column %d",  
                       errcode, j + 1);  
                                 /* inverse permute b */  
             for (i = 0; i < n; i++) out[i] = tmpOut[iperm[i]];  
         }  
117      }      }
118        if (piv) Free(tmp);
119      UNPROTECT(1);      UNPROTECT(1);
120      return val;      return val;
121  }  }
122    
123  SEXP sscMatrix_inverse_factor(SEXP A)  SEXP dsCMatrix_inverse_factor(SEXP A)
124  {  {
125      return mat_from_taucs(taucs_ccs_factor_xxt(      return R_NilValue;          /* FIXME: Write this function. */
         csc_taucs_ptr(A, TAUCS_DOUBLE | TAUCS_LOWER | TAUCS_SYMMETRIC)));  
126  }  }
127    
128  SEXP ssc_transpose(SEXP x)  SEXP ssc_transpose(SEXP x)
129  {  {
130      SEXP      SEXP
131          ans = PROTECT(NEW_OBJECT(MAKE_CLASS("sscMatrix"))),          ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dsCMatrix"))),
132          islot = GET_SLOT(x, Matrix_iSym);          islot = GET_SLOT(x, Matrix_iSym);
133      int nnz = length(islot),      int nnz = length(islot),
134          *adims = INTEGER(GET_SLOT(ans, Matrix_DimSym)),          *adims = INTEGER(GET_SLOT(ans, Matrix_DimSym)),
# Line 134  Line 151 
151      return ans;      return ans;
152  }  }
153    
154  SEXP sscMatrix_to_triplet(SEXP x)  SEXP dsCMatrix_to_dgTMatrix(SEXP x)
155  {  {
156      SEXP      SEXP
157          ans = PROTECT(NEW_OBJECT(MAKE_CLASS("tripletMatrix"))),          ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgTMatrix"))),
158          islot = GET_SLOT(x, Matrix_iSym),          islot = GET_SLOT(x, Matrix_iSym),
159          pslot = GET_SLOT(x, Matrix_pSym);          pslot = GET_SLOT(x, Matrix_pSym);
160      int *ai, *aj, *iv = INTEGER(islot),      int *ai, *aj, *iv = INTEGER(islot),
# Line 177  Line 194 
194      UNPROTECT(1);      UNPROTECT(1);
195      return ans;      return ans;
196  }  }
197    
198    SEXP dsCMatrix_ldl_symbolic(SEXP x, SEXP doPerm)
199    {
200        SEXP Ax, Dims = GET_SLOT(x, Matrix_DimSym),
201            ans = PROTECT(allocVector(VECSXP, 3)), tsc;
202        int i, n = INTEGER(Dims)[0], nz, nza,
203            *Ap, *Ai, *Lp, *Li, *Parent,
204            doperm = asLogical(doPerm),
205            *P = (int *) NULL, *Pinv = (int *) NULL;
206    
207    
208        if (toupper(CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0]) == 'L') {
209            x = PROTECT(ssc_transpose(x));
210        } else {
211            x = PROTECT(duplicate(x));
212        }
213        Ax = GET_SLOT(x, Matrix_xSym);
214        nza = length(Ax);
215        Ap = INTEGER(GET_SLOT(x, Matrix_pSym));
216        Ai = INTEGER(GET_SLOT(x, Matrix_iSym));
217        if (doperm) {
218            int *perm, *iperm = Calloc(n, int);
219    
220            SET_VECTOR_ELT(ans, 2, allocVector(INTSXP, n));
221            perm = INTEGER(VECTOR_ELT(ans, 2));
222            ssc_metis_order(n, Ap, Ai, perm, iperm);
223            ssc_symbolic_permute(n, 1, iperm, Ap, Ai);
224            Free(iperm);
225        }
226        SET_VECTOR_ELT(ans, 0, allocVector(INTSXP, n));
227        Parent = INTEGER(VECTOR_ELT(ans, 0));
228        SET_VECTOR_ELT(ans, 1, NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
229        tsc = VECTOR_ELT(ans, 1);
230        SET_SLOT(tsc, Matrix_uploSym, ScalarString(mkChar("L")));
231        SET_SLOT(tsc, Matrix_diagSym, ScalarString(mkChar("U")));
232        SET_SLOT(tsc, Matrix_DimSym, Dims);
233        SET_SLOT(tsc, Matrix_pSym, allocVector(INTSXP, n + 1));
234        Lp = INTEGER(GET_SLOT(tsc, Matrix_pSym));
235        R_ldl_symbolic(n, Ap, Ai, Lp, Parent, P, Pinv);
236        nz = Lp[n];
237        SET_SLOT(tsc, Matrix_iSym, allocVector(INTSXP, nz));
238        Li = INTEGER(GET_SLOT(tsc, Matrix_iSym));
239        SET_SLOT(tsc, Matrix_xSym, allocVector(REALSXP, nz));
240        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.;
242        i = R_ldl_numeric(n, Ap, Ai, REAL(Ax), Lp, Parent, Li,
243                        REAL(GET_SLOT(tsc, Matrix_xSym)),
244                        (double *) R_alloc(n, sizeof(double)), /* D */
245                        P, Pinv);
246        UNPROTECT(2);
247        return ans;
248    }
249    
250    SEXP dsCMatrix_metis_perm(SEXP x)
251    {
252        SEXP pSlot = GET_SLOT(x, Matrix_pSym),
253            ans = PROTECT(allocVector(VECSXP, 2));
254        int n = length(pSlot) - 1;
255    
256        SET_VECTOR_ELT(ans, 0, allocVector(INTSXP, n));
257        SET_VECTOR_ELT(ans, 1, allocVector(INTSXP, n));
258        ssc_metis_order(n,
259                        INTEGER(pSlot),
260                        INTEGER(GET_SLOT(x, Matrix_iSym)),
261                        INTEGER(VECTOR_ELT(ans, 0)),
262                        INTEGER(VECTOR_ELT(ans, 1)));
263        UNPROTECT(1);
264        return ans;
265    }
266    

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

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