SCM

SCM Repository

[matrix] Annotation of /pkg/src/dsCMatrix.c
ViewVC logotype

Annotation of /pkg/src/dsCMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 10 - (view) (download) (as text)
Original Path: pkg/src/sscMatrix.c

1 : bates 10 #include "sscMatrix.h"
2 :    
3 :     SEXP sscMatrix_validate(SEXP obj)
4 :     {
5 :     SEXP uplo = GET_SLOT(obj, Matrix_uploSym);
6 :     int *Dim = INTEGER(GET_SLOT(obj, Matrix_DimSym));
7 :     char *val;
8 :    
9 :     if (length(uplo) != 1)
10 :     return ScalarString(mkChar("uplo slot must have length 1"));
11 :     val = CHAR(STRING_ELT(uplo, 0));
12 :     if (strlen(val) != 1)
13 :     return ScalarString(mkChar("uplo[1] must have string length 1"));
14 :     if (toupper(*val) != 'U' && toupper(*val) != 'L')
15 :     return ScalarString(mkChar("uplo[1] must be \"U\" or \"L\""));
16 :     if (Dim[0] != Dim[1])
17 :     return ScalarString(mkChar("Symmetric matrix must be square"));
18 :     csc_check_column_sorting(obj);
19 :     return ScalarLogical(1);
20 :     }
21 :    
22 :     SEXP sscMatrix_chol(SEXP x, SEXP pivot)
23 :     {
24 :     SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("sscChol")));
25 :     taucs_ccs_matrix *tm =
26 :     csc_taucs_ptr(x, TAUCS_DOUBLE | TAUCS_LOWER | TAUCS_SYMMETRIC);
27 :     int nnz, piv = asLogical(pivot);
28 :     int *dims;
29 :     void *L;
30 :    
31 :     if (piv) {
32 :     int *iperm, *perm;
33 :    
34 :     SET_SLOT(val, Matrix_permSym, allocVector(INTSXP, tm->n));
35 :     SET_SLOT(val, Matrix_ipermSym, allocVector(INTSXP, tm->n));
36 :     perm = INTEGER(GET_SLOT(val, Matrix_permSym));
37 :     iperm = INTEGER(GET_SLOT(val, Matrix_ipermSym));
38 :     ssc_metis_order(tm->n, (tm->colptr)[tm->n], tm->colptr,
39 :     tm->rowind, perm, iperm);
40 :     tm = taucs_dccs_permute_symmetrically(tm, perm, iperm);
41 :     }
42 :     if (!(L = taucs_ccs_factor_llt_mf(tm)))
43 :     error("Matrix is not positive definite");
44 :     if (piv) taucs_dccs_free(tm);
45 :     tm = taucs_supernodal_factor_to_ccs(L);
46 :     taucs_supernodal_factor_free(L);
47 :     nnz = tm->colptr[tm->n];
48 :     SET_SLOT(val, Matrix_pSym, allocVector(INTSXP, tm->n + 1));
49 :     Memcpy(INTEGER(GET_SLOT(val, Matrix_pSym)), tm->colptr, tm->n + 1);
50 :     SET_SLOT(val, Matrix_iSym, allocVector(INTSXP, nnz));
51 :     Memcpy(INTEGER(GET_SLOT(val, Matrix_iSym)), tm->rowind, nnz);
52 :     SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, nnz));
53 :     Memcpy(REAL(GET_SLOT(val, Matrix_xSym)), tm->values.d, nnz);
54 :     dims = INTEGER(GET_SLOT(val, Matrix_DimSym));
55 :     dims[0] = dims[1] = tm->n;
56 :     taucs_dccs_free(tm);
57 :     UNPROTECT(1);
58 :     return set_factorization(x, val, "Cholesky");
59 :     }
60 :    
61 :     SEXP sscMatrix_matrix_solve(SEXP a, SEXP b)
62 :     {
63 :     SEXP Chol = get_factorization(a, "Cholesky"),
64 :     val = PROTECT(duplicate(b));
65 :     int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
66 :     *bdims = INTEGER(getAttrib(b, R_DimSymbol)),
67 :     j, n = adims[1];
68 :     taucs_ccs_matrix* tm;
69 :     double *in = REAL(b), *out = REAL(val) ;
70 :    
71 :     if (!(isReal(b) && isMatrix(b)))
72 :     error("Argument b must be a numeric matrix");
73 :     if (*adims != *bdims || bdims[1] < 1 || *adims < 1)
74 :     error("Dimensions of system to be solved are inconsistent");
75 :     if (Chol == R_NilValue) Chol = sscMatrix_chol(a, ScalarLogical(1.));
76 :     tm = csc_taucs_ptr(Chol, TAUCS_DOUBLE | TAUCS_LOWER | TAUCS_TRIANGULAR);
77 :     if (!length(GET_SLOT(Chol, Matrix_permSym))) {
78 :     for (j = 0; j < bdims[1]; j++, in += n, out += n) {
79 :     int errcode = taucs_dccs_solve_llt(tm, out, in);
80 :     if (errcode)
81 :     error("taucs_solve returned error code %d for column %d",
82 :     errcode, j + 1);
83 :     }
84 :     } else {
85 :     int *iperm = INTEGER(GET_SLOT(Chol, Matrix_ipermSym));
86 :     double *tmpIn = (double *) R_alloc(n, sizeof(double)),
87 :     *tmpOut = (double *) R_alloc(n, sizeof(double));
88 :    
89 :     for (j = 0; j < bdims[1]; j++, in += n, out += n) {
90 :     int errcode, i;
91 :     /* permute y */
92 :     for (i = 0; i < n; i++) tmpIn[iperm[i]] = in[i];
93 :     /* solve */
94 :     errcode = taucs_dccs_solve_llt(tm, tmpOut, tmpIn);
95 :     if (errcode)
96 :     error("taucs_solve returned error code %d for column %d",
97 :     errcode, j + 1);
98 :     /* inverse permute b */
99 :     for (i = 0; i < n; i++) out[i] = tmpOut[iperm[i]];
100 :     }
101 :     }
102 :     UNPROTECT(1);
103 :     return val;
104 :     }
105 :    
106 :     SEXP sscMatrix_inverse_factor(SEXP A)
107 :     {
108 :     return mat_from_taucs(taucs_ccs_factor_xxt(
109 :     csc_taucs_ptr(A, TAUCS_DOUBLE | TAUCS_LOWER | TAUCS_SYMMETRIC)));
110 :     }
111 :    
112 :     SEXP ssc_transpose(SEXP x)
113 :     {
114 :     SEXP
115 :     ans = PROTECT(NEW_OBJECT(MAKE_CLASS("sscMatrix"))),
116 :     islot = GET_SLOT(x, Matrix_iSym);
117 :     int nnz = length(islot),
118 :     *adims = INTEGER(GET_SLOT(ans, Matrix_DimSym)),
119 :     *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym));
120 :    
121 :     adims[0] = xdims[1]; adims[1] = xdims[0];
122 :     if (toupper(CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0]) == 'U')
123 :     SET_SLOT(ans, Matrix_uploSym, ScalarString(mkChar("L")));
124 :     SET_SLOT(ans, Matrix_pSym, allocVector(INTSXP, xdims[0] + 1));
125 :     SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nnz));
126 :     SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nnz));
127 :     csc_components_transpose(xdims[0], xdims[1], nnz,
128 :     INTEGER(GET_SLOT(x, Matrix_pSym)),
129 :     INTEGER(islot),
130 :     REAL(GET_SLOT(x, Matrix_xSym)),
131 :     INTEGER(GET_SLOT(ans, Matrix_pSym)),
132 :     INTEGER(GET_SLOT(ans, Matrix_iSym)),
133 :     REAL(GET_SLOT(ans, Matrix_xSym)));
134 :     UNPROTECT(1);
135 :     return ans;
136 :     }

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