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 70 - (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 : bates 70 ssc_metis_order(tm->n, tm->colptr, tm->rowind, perm, iperm);
39 : bates 10 tm = taucs_dccs_permute_symmetrically(tm, perm, iperm);
40 :     }
41 :     if (!(L = taucs_ccs_factor_llt_mf(tm)))
42 :     error("Matrix is not positive definite");
43 :     if (piv) taucs_dccs_free(tm);
44 :     tm = taucs_supernodal_factor_to_ccs(L);
45 :     taucs_supernodal_factor_free(L);
46 :     nnz = tm->colptr[tm->n];
47 :     SET_SLOT(val, Matrix_pSym, allocVector(INTSXP, tm->n + 1));
48 :     Memcpy(INTEGER(GET_SLOT(val, Matrix_pSym)), tm->colptr, tm->n + 1);
49 :     SET_SLOT(val, Matrix_iSym, allocVector(INTSXP, nnz));
50 :     Memcpy(INTEGER(GET_SLOT(val, Matrix_iSym)), tm->rowind, nnz);
51 :     SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, nnz));
52 :     Memcpy(REAL(GET_SLOT(val, Matrix_xSym)), tm->values.d, nnz);
53 :     dims = INTEGER(GET_SLOT(val, Matrix_DimSym));
54 :     dims[0] = dims[1] = tm->n;
55 :     taucs_dccs_free(tm);
56 :     UNPROTECT(1);
57 :     return set_factorization(x, val, "Cholesky");
58 :     }
59 :    
60 :     SEXP sscMatrix_matrix_solve(SEXP a, SEXP b)
61 :     {
62 :     SEXP Chol = get_factorization(a, "Cholesky"),
63 :     val = PROTECT(duplicate(b));
64 :     int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
65 :     *bdims = INTEGER(getAttrib(b, R_DimSymbol)),
66 :     j, n = adims[1];
67 :     taucs_ccs_matrix* tm;
68 :     double *in = REAL(b), *out = REAL(val) ;
69 :    
70 :     if (!(isReal(b) && isMatrix(b)))
71 :     error("Argument b must be a numeric matrix");
72 :     if (*adims != *bdims || bdims[1] < 1 || *adims < 1)
73 :     error("Dimensions of system to be solved are inconsistent");
74 :     if (Chol == R_NilValue) Chol = sscMatrix_chol(a, ScalarLogical(1.));
75 :     tm = csc_taucs_ptr(Chol, TAUCS_DOUBLE | TAUCS_LOWER | TAUCS_TRIANGULAR);
76 :     if (!length(GET_SLOT(Chol, Matrix_permSym))) {
77 :     for (j = 0; j < bdims[1]; j++, in += n, out += n) {
78 :     int errcode = taucs_dccs_solve_llt(tm, out, in);
79 :     if (errcode)
80 :     error("taucs_solve returned error code %d for column %d",
81 :     errcode, j + 1);
82 :     }
83 :     } else {
84 :     int *iperm = INTEGER(GET_SLOT(Chol, Matrix_ipermSym));
85 :     double *tmpIn = (double *) R_alloc(n, sizeof(double)),
86 :     *tmpOut = (double *) R_alloc(n, sizeof(double));
87 :    
88 :     for (j = 0; j < bdims[1]; j++, in += n, out += n) {
89 :     int errcode, i;
90 :     /* permute y */
91 :     for (i = 0; i < n; i++) tmpIn[iperm[i]] = in[i];
92 :     /* solve */
93 :     errcode = taucs_dccs_solve_llt(tm, tmpOut, tmpIn);
94 :     if (errcode)
95 :     error("taucs_solve returned error code %d for column %d",
96 :     errcode, j + 1);
97 :     /* inverse permute b */
98 :     for (i = 0; i < n; i++) out[i] = tmpOut[iperm[i]];
99 :     }
100 :     }
101 :     UNPROTECT(1);
102 :     return val;
103 :     }
104 :    
105 :     SEXP sscMatrix_inverse_factor(SEXP A)
106 :     {
107 :     return mat_from_taucs(taucs_ccs_factor_xxt(
108 :     csc_taucs_ptr(A, TAUCS_DOUBLE | TAUCS_LOWER | TAUCS_SYMMETRIC)));
109 :     }
110 :    
111 :     SEXP ssc_transpose(SEXP x)
112 :     {
113 :     SEXP
114 :     ans = PROTECT(NEW_OBJECT(MAKE_CLASS("sscMatrix"))),
115 :     islot = GET_SLOT(x, Matrix_iSym);
116 :     int nnz = length(islot),
117 :     *adims = INTEGER(GET_SLOT(ans, Matrix_DimSym)),
118 :     *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym));
119 :    
120 :     adims[0] = xdims[1]; adims[1] = xdims[0];
121 :     if (toupper(CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0]) == 'U')
122 :     SET_SLOT(ans, Matrix_uploSym, ScalarString(mkChar("L")));
123 :     SET_SLOT(ans, Matrix_pSym, allocVector(INTSXP, xdims[0] + 1));
124 :     SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nnz));
125 :     SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nnz));
126 :     csc_components_transpose(xdims[0], xdims[1], nnz,
127 :     INTEGER(GET_SLOT(x, Matrix_pSym)),
128 :     INTEGER(islot),
129 :     REAL(GET_SLOT(x, Matrix_xSym)),
130 :     INTEGER(GET_SLOT(ans, Matrix_pSym)),
131 :     INTEGER(GET_SLOT(ans, Matrix_iSym)),
132 :     REAL(GET_SLOT(ans, Matrix_xSym)));
133 :     UNPROTECT(1);
134 :     return ans;
135 :     }
136 : bates 70
137 :     SEXP sscMatrix_to_triplet(SEXP x)
138 :     {
139 :     SEXP
140 :     ans = PROTECT(NEW_OBJECT(MAKE_CLASS("tripletMatrix"))),
141 :     islot = GET_SLOT(x, Matrix_iSym),
142 :     pslot = GET_SLOT(x, Matrix_pSym);
143 :     int *ai, *aj, *iv = INTEGER(islot),
144 :     j, jj, nnz = length(islot), nout,
145 :     n = length(pslot) - 1,
146 :     *p = INTEGER(pslot), pos;
147 :     double *ax, *xv = REAL(GET_SLOT(x, Matrix_xSym));
148 :    
149 :     /* increment output count by number of off-diagonals */
150 :     nout = nnz;
151 :     for (j = 0; j < n; j++) {
152 :     int p2 = p[j+1];
153 :     for (jj = p[j]; jj < p2; jj++) {
154 :     if (iv[jj] != j) nout++;
155 :     }
156 :     }
157 :     SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(x, Matrix_DimSym)));
158 :     SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nout));
159 :     ai = INTEGER(GET_SLOT(ans, Matrix_iSym));
160 :     SET_SLOT(ans, Matrix_jSym, allocVector(INTSXP, nout));
161 :     aj = INTEGER(GET_SLOT(ans, Matrix_jSym));
162 :     SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nout));
163 :     ax = REAL(GET_SLOT(ans, Matrix_xSym));
164 :     pos = 0;
165 :     for (j = 0; j < n; j++) {
166 :     int p2 = p[j+1];
167 :     for (jj = p[j]; jj < p2; jj++) {
168 :     int ii = iv[jj];
169 :     double xx = xv[jj];
170 :    
171 :     ai[pos] = ii; aj[pos] = j; ax[pos] = xx; pos++;
172 :     if (ii != j) {
173 :     aj[pos] = ii; ai[pos] = j; ax[pos] = xx; pos++;
174 :     }
175 :     }
176 :     }
177 :     UNPROTECT(1);
178 :     return ans;
179 :     }

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