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 10, Mon Mar 22 20:20:05 2004 UTC revision 209, Thu Jun 3 18:22:56 2004 UTC
# Line 35  Line 35 
35          SET_SLOT(val, Matrix_ipermSym, allocVector(INTSXP, tm->n));          SET_SLOT(val, Matrix_ipermSym, allocVector(INTSXP, tm->n));
36          perm = INTEGER(GET_SLOT(val, Matrix_permSym));          perm = INTEGER(GET_SLOT(val, Matrix_permSym));
37          iperm = INTEGER(GET_SLOT(val, Matrix_ipermSym));          iperm = INTEGER(GET_SLOT(val, Matrix_ipermSym));
38          ssc_metis_order(tm->n, (tm->colptr)[tm->n], tm->colptr,          ssc_metis_order(tm->n, tm->colptr, tm->rowind, perm, iperm);
                         tm->rowind, perm, iperm);  
39          tm = taucs_dccs_permute_symmetrically(tm, perm, iperm);          tm = taucs_dccs_permute_symmetrically(tm, perm, iperm);
40      }      }
41      if (!(L = taucs_ccs_factor_llt_mf(tm)))      if (!(L = taucs_ccs_factor_llt_mf(tm)))
# Line 134  Line 133 
133      UNPROTECT(1);      UNPROTECT(1);
134      return ans;      return ans;
135  }  }
136    
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    }
180    
181    SEXP sscMatrix_ldl_symbolic(SEXP x, SEXP doPerm)
182    {
183        SEXP Ax, Dims = GET_SLOT(x, Matrix_DimSym),
184            ans = PROTECT(allocVector(VECSXP, 3)), tsc;
185        int i, n = INTEGER(Dims)[0], nz, nza,
186            *Ap, *Ai, *Lp, *Li, *Parent,
187            doperm = asLogical(doPerm),
188            *Lnz = (int *) R_alloc(n, sizeof(int)),
189            *Flag = (int *) R_alloc(n, sizeof(int)),
190            *P = (int *) NULL, *Pinv = (int *) NULL;
191    
192    
193        if (toupper(CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0]) == 'L') {
194            x = PROTECT(ssc_transpose(x));
195        } else {
196            x = PROTECT(duplicate(x));
197        }
198        Ax = GET_SLOT(x, Matrix_xSym);
199        nza = length(Ax);
200        Ap = INTEGER(GET_SLOT(x, Matrix_pSym));
201        Ai = INTEGER(GET_SLOT(x, Matrix_iSym));
202        if (doperm) {
203            int *perm;
204            SET_VECTOR_ELT(ans, 2, allocVector(INTSXP, n));
205            perm = INTEGER(VECTOR_ELT(ans, 2));
206            ssc_metis_order(n, Ap, Ai, perm, Flag);
207            ssc_symbolic_permute(n, 1, Flag, Ap, Ai);
208        }
209        SET_VECTOR_ELT(ans, 0, allocVector(INTSXP, n));
210        Parent = INTEGER(VECTOR_ELT(ans, 0));
211        SET_VECTOR_ELT(ans, 1, NEW_OBJECT(MAKE_CLASS("tscMatrix")));
212        tsc = VECTOR_ELT(ans, 1);
213        SET_SLOT(tsc, Matrix_uploSym, ScalarString(mkChar("L")));
214        SET_SLOT(tsc, Matrix_diagSym, ScalarString(mkChar("U")));
215        SET_SLOT(tsc, Matrix_DimSym, Dims);
216        SET_SLOT(tsc, Matrix_pSym, allocVector(INTSXP, n + 1));
217        Lp = INTEGER(GET_SLOT(tsc, Matrix_pSym));
218        ldl_symbolic(n, Ap, Ai, Lp, Parent, Lnz, Flag, P, Pinv);
219        nz = Lp[n];
220        SET_SLOT(tsc, Matrix_iSym, allocVector(INTSXP, nz));
221        Li = INTEGER(GET_SLOT(tsc, Matrix_iSym));
222        SET_SLOT(tsc, Matrix_xSym, allocVector(REALSXP, nz));
223        for (i = 0; i < nza; i++) REAL(Ax)[i] = 0.00001;
224        for (i = 0; i < n; i++) REAL(Ax)[Ap[i+1]-1] = 10000.;
225        i = ldl_numeric(n, Ap, Ai, REAL(Ax), Lp, Parent, Lnz, Li,
226                        REAL(GET_SLOT(tsc, Matrix_xSym)),
227                        (double *) R_alloc(n, sizeof(double)), /* D */
228                        (double *) R_alloc(n, sizeof(double)), /* Y */
229                        (int *) R_alloc(n, sizeof(int)), /* Pattern */
230                        Flag, P, Pinv);
231        UNPROTECT(2);
232        return ans;
233    }
234    
235    SEXP sscMatrix_metis_perm(SEXP x)
236    {
237        SEXP pSlot = GET_SLOT(x, Matrix_pSym),
238            ans = PROTECT(allocVector(VECSXP, 2));
239        int n = length(pSlot) - 1;
240    
241        SET_VECTOR_ELT(ans, 0, allocVector(INTSXP, n));
242        SET_VECTOR_ELT(ans, 1, allocVector(INTSXP, n));
243        ssc_metis_order(n,
244                        INTEGER(pSlot),
245                        INTEGER(GET_SLOT(x, Matrix_iSym)),
246                        INTEGER(VECTOR_ELT(ans, 0)),
247                        INTEGER(VECTOR_ELT(ans, 1)));
248        UNPROTECT(1);
249        return ans;
250    }
251    
252    SEXP sscMatrix_metis_ldl_symbolic(SEXP x)
253    {
254        SEXP pSlot = GET_SLOT(x, Matrix_pSym),
255            ans = PROTECT(allocVector(VECSXP, 4));
256        int *Ai = INTEGER(GET_SLOT(x, Matrix_iSym)),
257            lo = toupper(CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0]) == 'L',
258            n = length(pSlot)-1;
259    
260    
261        if (lo) x = PROTECT(ssc_transpose(x));
262        SET_VECTOR_ELT(ans, 0, allocVector(INTSXP, n));
263        SET_VECTOR_ELT(ans, 1, allocVector(INTSXP, n));
264        SET_VECTOR_ELT(ans, 2, allocVector(INTSXP, n + 1));
265        SET_VECTOR_ELT(ans, 3, allocVector(INTSXP, n));
266        ssc_metis_order(n, INTEGER(pSlot), Ai,
267                        INTEGER(VECTOR_ELT(ans, 0)), /* P */
268                        INTEGER(VECTOR_ELT(ans, 1))); /* Pinv */
269        ldl_symbolic(n, INTEGER(pSlot), Ai,
270                     INTEGER(VECTOR_ELT(ans, 2)), /* Lp */
271                     INTEGER(VECTOR_ELT(ans, 3)), /* Parent */
272                     (int *) R_alloc(n, sizeof(int)), /* Lnz */
273                     (int *) R_alloc(n, sizeof(int)), /* Flag */
274                     INTEGER(VECTOR_ELT(ans, 0)), /* P */
275                     INTEGER(VECTOR_ELT(ans, 1))); /* Pinv */
276        UNPROTECT(lo ? 2 : 1);
277        return ans;
278    }

Legend:
Removed from v.10  
changed lines
  Added in v.209

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