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 153, Fri May 7 19:35:35 2004 UTC revision 308, Wed Oct 27 21:39:44 2004 UTC
# Line 21  Line 21 
21    
22  SEXP sscMatrix_chol(SEXP x, SEXP pivot)  SEXP sscMatrix_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 *Flag = Calloc(n, int), *Lnz = Calloc(n, int),
33          csc_taucs_ptr(x, TAUCS_DOUBLE | TAUCS_LOWER | TAUCS_SYMMETRIC);          *P = (int *) NULL, *Pinv = (int *) NULL;
34      int nnz, piv = asLogical(pivot);      double *Ax;
35      int *dims;  
36      void *L;      if (lo) {
37            x = PROTECT(ssc_transpose(x));
38            Ai = INTEGER(GET_SLOT(x, Matrix_iSym));
39            Ap = INTEGER(GET_SLOT(x, Matrix_pSym));
40        }
41        SET_SLOT(val, Matrix_uploSym, ScalarString(mkChar("L")));
42        SET_SLOT(val, Matrix_diagSym, ScalarString(mkChar("N")));
43        SET_SLOT(val, Matrix_DimSym, duplicate(GET_SLOT(x, Matrix_DimSym)));
44        SET_SLOT(val, Matrix_ParentSym, allocVector(INTSXP, n));
45        Parent = INTEGER(GET_SLOT(val, Matrix_ParentSym));
46        SET_SLOT(val, Matrix_pSym, allocVector(INTSXP, n + 1));
47        Lp = INTEGER(GET_SLOT(val, Matrix_pSym));
48        Ax = REAL(GET_SLOT(x, Matrix_xSym));
49      if (piv) {      if (piv) {
50          int *iperm, *perm;          SEXP trip = PROTECT(sscMatrix_to_triplet(x));
51            SEXP Ti = GET_SLOT(trip, Matrix_iSym);
52    
53          SET_SLOT(val, Matrix_permSym, allocVector(INTSXP, tm->n));          /* determine the permutation with Metis */
54          SET_SLOT(val, Matrix_ipermSym, allocVector(INTSXP, tm->n));          Pinv = Calloc(n, int);
55          perm = INTEGER(GET_SLOT(val, Matrix_permSym));          SET_SLOT(val, Matrix_permSym, allocVector(INTSXP, n));
56          iperm = INTEGER(GET_SLOT(val, Matrix_ipermSym));          P = INTEGER(GET_SLOT(val, Matrix_permSym));
57          ssc_metis_order(tm->n, tm->colptr, tm->rowind, perm, iperm);          ssc_metis_order(n, Ap, Ai, P, Pinv);
58          tm = taucs_dccs_permute_symmetrically(tm, perm, iperm);          /* create a symmetrized form of x */
59      }          nnz = length(Ti);
60      if (!(L = taucs_ccs_factor_llt_mf(tm)))          Ai = Calloc(nnz, int);
61          error("Matrix is not positive definite");          Ax = Calloc(nnz, double);
62      if (piv) taucs_dccs_free(tm);          Ap = Calloc(n + 1, int);
63      tm = taucs_supernodal_factor_to_ccs(L);          triplet_to_col(n, n, nnz, INTEGER(Ti),
64      taucs_supernodal_factor_free(L);                         INTEGER(GET_SLOT(trip, Matrix_jSym)),
65      nnz = tm->colptr[tm->n];                         REAL(GET_SLOT(trip, Matrix_xSym)),
66      SET_SLOT(val, Matrix_pSym, allocVector(INTSXP, tm->n + 1));                         Ap, Ai, Ax);
67      Memcpy(INTEGER(GET_SLOT(val, Matrix_pSym)), tm->colptr, tm->n + 1);      }
68        ldl_symbolic(n, Ap, Ai, Lp, Parent, Lnz, Flag, P, Pinv);
69        nnz = Lp[n];
70      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);  
71      SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, nnz));      SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, nnz));
72      Memcpy(REAL(GET_SLOT(val, Matrix_xSym)), tm->values.d, nnz);      SET_SLOT(val, Matrix_DSym, allocVector(REALSXP, n));
73      dims = INTEGER(GET_SLOT(val, Matrix_DimSym));      info = ldl_numeric(n, Ap, Ai, Ax,
74      dims[0] = dims[1] = tm->n;                         Lp, Parent, Lnz,
75      taucs_dccs_free(tm);                         INTEGER(GET_SLOT(val, Matrix_iSym)),
76                           REAL(GET_SLOT(val, Matrix_xSym)),
77                           REAL(GET_SLOT(val, Matrix_DSym)),
78                           (double *) R_alloc(n, sizeof(double)), /* Y */
79                           (int *) R_alloc(n, sizeof(int)), /* Pattern */
80                           Flag, P, Pinv);
81        if (info != n)
82            error("Leading minor of size %d (possibly after permutation) is indefinite",
83                  info + 1);
84        Free(Flag); Free(Lnz);
85        if (piv) {
86      UNPROTECT(1);      UNPROTECT(1);
87      return set_factorization(x, val, "Cholesky");          Free(Pinv); Free(Ax); Free(Ai); Free(Ap);
88        }
89        UNPROTECT(lo ? 2 : 1);
90        return set_factorization(xorig, val, "Cholesky");
91  }  }
92    
93  SEXP sscMatrix_matrix_solve(SEXP a, SEXP b)  SEXP sscMatrix_matrix_solve(SEXP a, SEXP b)
94  {  {
95      SEXP Chol = get_factorization(a, "Cholesky"),      SEXP Chol = get_factorization(a, "Cholesky"), perm,
96          val = PROTECT(duplicate(b));          val = PROTECT(duplicate(b));
97      int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),      int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
98          *bdims = INTEGER(getAttrib(b, R_DimSymbol)),          *bdims = INTEGER(getAttrib(b, R_DimSymbol)),
99          j, n = adims[1];          *Li, *Lp, j, n = adims[1], ncol = bdims[1], piv;
100      taucs_ccs_matrix* tm;      double *Lx, *D, *in = REAL(b), *out = REAL(val), *tmp = (double *) NULL;
     double *in = REAL(b), *out = REAL(val) ;  
101    
102      if (!(isReal(b) && isMatrix(b)))      if (!(isReal(b) && isMatrix(b)))
103          error("Argument b must be a numeric matrix");          error("Argument b must be a numeric matrix");
104      if (*adims != *bdims || bdims[1] < 1 || *adims < 1)      if (*adims != *bdims || bdims[1] < 1 || *adims < 1)
105          error("Dimensions of system to be solved are inconsistent");          error("Dimensions of system to be solved are inconsistent");
106      if (Chol == R_NilValue) Chol = sscMatrix_chol(a, ScalarLogical(1.));      if (Chol == R_NilValue) Chol = sscMatrix_chol(a, ScalarLogical(1));
107      tm = csc_taucs_ptr(Chol, TAUCS_DOUBLE | TAUCS_LOWER | TAUCS_TRIANGULAR);      perm = GET_SLOT(Chol, Matrix_permSym);
108      if (!length(GET_SLOT(Chol, Matrix_permSym))) {      piv = length(perm);
109          for (j = 0; j < bdims[1]; j++, in += n, out += n) {      if (piv) tmp = Calloc(n, double);
110              int errcode = taucs_dccs_solve_llt(tm, out, in);      Li = INTEGER(GET_SLOT(Chol, Matrix_iSym));
111              if (errcode)      Lp = INTEGER(GET_SLOT(Chol, Matrix_pSym));
112                  error("taucs_solve returned error code %d for column %d",      Lx = REAL(GET_SLOT(Chol, Matrix_xSym));
113                        errcode, j + 1);      D = REAL(GET_SLOT(Chol, Matrix_DSym));
114          }      for (j = 0; j < ncol; j++, in += n, out += n) {
115      } else {          if (piv) ldl_perm(n, out, in, INTEGER(perm));
116          int *iperm = INTEGER(GET_SLOT(Chol, Matrix_ipermSym));          else Memcpy(out, in, n);
117          double *tmpIn = (double *) R_alloc(n, sizeof(double)),          ldl_lsolve(n, out, Lp, Li, Lx);
118              *tmpOut = (double *) R_alloc(n, sizeof(double));          ldl_dsolve(n, out, D);
119            ldl_ltsolve(n, out, Lp, Li, Lx);
120          for (j = 0; j < bdims[1]; j++, in += n, out += n) {          if (piv) 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]];  
         }  
121      }      }
122        if (piv) Free(tmp);
123      UNPROTECT(1);      UNPROTECT(1);
124      return val;      return val;
125  }  }
126    
127  SEXP sscMatrix_inverse_factor(SEXP A)  SEXP sscMatrix_inverse_factor(SEXP A)
128  {  {
129      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)));  
130  }  }
131    
132  SEXP ssc_transpose(SEXP x)  SEXP ssc_transpose(SEXP x)
# Line 178  Line 199 
199      return ans;      return ans;
200  }  }
201    
202  SEXP sscMatrix_ldl_symbolic(SEXP x)  SEXP sscMatrix_ldl_symbolic(SEXP x, SEXP doPerm)
203  {  {
204      SEXP ans = PROTECT(allocVector(VECSXP, 2));      SEXP Ax, Dims = GET_SLOT(x, Matrix_DimSym),
205      int lo = toupper(CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0]) == 'L',          ans = PROTECT(allocVector(VECSXP, 3)), tsc;
206          n = INTEGER(GET_SLOT(x, Matrix_DimSym))[0];      int i, n = INTEGER(Dims)[0], nz, nza,
207            *Ap, *Ai, *Lp, *Li, *Parent,
208            doperm = asLogical(doPerm),
209            *Lnz = (int *) R_alloc(n, sizeof(int)),
210            *Flag = (int *) R_alloc(n, sizeof(int)),
211            *P = (int *) NULL, *Pinv = (int *) NULL;
212    
213    
214        if (toupper(CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0]) == 'L') {
215            x = PROTECT(ssc_transpose(x));
216        } else {
217            x = PROTECT(duplicate(x));
218        }
219        Ax = GET_SLOT(x, Matrix_xSym);
220        nza = length(Ax);
221        Ap = INTEGER(GET_SLOT(x, Matrix_pSym));
222        Ai = INTEGER(GET_SLOT(x, Matrix_iSym));
223        if (doperm) {
224            int *perm;
225            SET_VECTOR_ELT(ans, 2, allocVector(INTSXP, n));
226            perm = INTEGER(VECTOR_ELT(ans, 2));
227            ssc_metis_order(n, Ap, Ai, perm, Flag);
228            ssc_symbolic_permute(n, 1, Flag, Ap, Ai);
229        }
230        SET_VECTOR_ELT(ans, 0, allocVector(INTSXP, n));
231        Parent = INTEGER(VECTOR_ELT(ans, 0));
232        SET_VECTOR_ELT(ans, 1, NEW_OBJECT(MAKE_CLASS("tscMatrix")));
233        tsc = VECTOR_ELT(ans, 1);
234        SET_SLOT(tsc, Matrix_uploSym, ScalarString(mkChar("L")));
235        SET_SLOT(tsc, Matrix_diagSym, ScalarString(mkChar("U")));
236        SET_SLOT(tsc, Matrix_DimSym, Dims);
237        SET_SLOT(tsc, Matrix_pSym, allocVector(INTSXP, n + 1));
238        Lp = INTEGER(GET_SLOT(tsc, Matrix_pSym));
239        ldl_symbolic(n, Ap, Ai, Lp, Parent, Lnz, Flag, P, Pinv);
240        nz = Lp[n];
241        SET_SLOT(tsc, Matrix_iSym, allocVector(INTSXP, nz));
242        Li = INTEGER(GET_SLOT(tsc, Matrix_iSym));
243        SET_SLOT(tsc, Matrix_xSym, allocVector(REALSXP, nz));
244        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.;
246        i = ldl_numeric(n, Ap, Ai, REAL(Ax), Lp, Parent, Lnz, Li,
247                        REAL(GET_SLOT(tsc, Matrix_xSym)),
248                        (double *) R_alloc(n, sizeof(double)), /* D */
249                        (double *) R_alloc(n, sizeof(double)), /* Y */
250                        (int *) R_alloc(n, sizeof(int)), /* Pattern */
251                        Flag, P, Pinv);
252        UNPROTECT(2);
253        return ans;
254    }
255    
256    SEXP sscMatrix_metis_perm(SEXP x)
257    {
258        SEXP pSlot = GET_SLOT(x, Matrix_pSym),
259            ans = PROTECT(allocVector(VECSXP, 2));
260        int n = length(pSlot) - 1;
261    
     if (lo) x = PROTECT(ssc_transpose(x));  
262      SET_VECTOR_ELT(ans, 0, allocVector(INTSXP, n));      SET_VECTOR_ELT(ans, 0, allocVector(INTSXP, n));
263      SET_VECTOR_ELT(ans, 1, allocVector(INTSXP, n + 1));      SET_VECTOR_ELT(ans, 1, allocVector(INTSXP, n));
264      ldl_symbolic(n, INTEGER(GET_SLOT(x, Matrix_pSym)),      ssc_metis_order(n,
265                        INTEGER(pSlot),
266                   INTEGER(GET_SLOT(x, Matrix_iSym)),                   INTEGER(GET_SLOT(x, Matrix_iSym)),
267                   INTEGER(VECTOR_ELT(ans, 1)), /* Lp */                      INTEGER(VECTOR_ELT(ans, 0)),
268                   INTEGER(VECTOR_ELT(ans, 0)), /* Parent */                      INTEGER(VECTOR_ELT(ans, 1)));
269                   (int *) R_alloc(n, sizeof(int)), /* Lnz */      UNPROTECT(1);
                  (int *) R_alloc(n, sizeof(int)), /* Flag */  
                  (int *) NULL, (int *) NULL);  /* P & Pinv */  
     UNPROTECT(lo ? 2 : 1);  
270      return ans;      return ans;
271  }  }
272    

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

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