SCM

SCM Repository

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

Annotation of /pkg/src/Pattern.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 922 - (view) (download) (as text)

1 : bates 912 #include "Pattern.h"
2 :    
3 :     cholmod_sparse *factor_to_pattern(SEXP fact)
4 :     {
5 :     int *vals = INTEGER(fact), i, n = LENGTH(fact),
6 :     nlev = LENGTH(getAttrib(fact, R_LevelsSymbol));
7 :     cholmod_sparse *A;
8 :    
9 :     if (!isFactor(fact))
10 :     error("`fact' must be a factor");
11 :     A = cholmod_allocate_sparse((size_t) nlev, (size_t) n, (size_t) n,
12 :     TRUE, TRUE, 0, CHOLMOD_PATTERN, &c);
13 :     for (i = 0; i <= n; i++) ((int *) (A->p))[i] = i;
14 :     for (i = 0; i < n; i++) ((int *) (A->i))[i] = vals[i] - 1;
15 :     if (cholmod_check_sparse(A, &c) < 0)
16 :     error("sparse Pattern matrix has invalid structure");
17 :     return A;
18 :     }
19 :    
20 :     SEXP factor_prod(SEXP f1, SEXP f2)
21 :     {
22 :     SEXP ans = PROTECT(allocVector(VECSXP, 3));
23 : bates 922 int *dims, n = LENGTH(f1), i, nl1, nl2, nmax, super;
24 : bates 912 cholmod_triplet *A;
25 :     cholmod_sparse *B;
26 :     cholmod_factor *F;
27 :    
28 :     if (!isFactor(f1) || !isFactor(f2) || LENGTH(f2) != n)
29 :     error("f1 and f2 must be factors of the same length");
30 :     nl1 = LENGTH(getAttrib(f1, R_LevelsSymbol));
31 :     nl2 = LENGTH(getAttrib(f2, R_LevelsSymbol));
32 :    
33 :     A = cholmod_allocate_triplet((size_t) nl1, (size_t) nl2, (size_t) n,
34 :     0, CHOLMOD_PATTERN, &c);
35 :     for (i = 0; i < n; i++) {
36 :     ((int *)(A->i))[i] = INTEGER(f1)[i] - 1;
37 :     ((int *)(A->j))[i] = INTEGER(f2)[i] - 1;
38 :     }
39 :     A->nnz = n;
40 :     nmax = nl1 * nl2;
41 :     if (n < nmax) nmax = n;
42 :     B = cholmod_triplet_to_sparse(A, nmax, &c);
43 :     cholmod_free_triplet(&A, &c);
44 :     /* force a simplicial factorization */
45 :     super = c.supernodal;
46 :     c.supernodal = CHOLMOD_SIMPLICIAL;
47 :     F = cholmod_analyze(B, &c);
48 :     c.supernodal = super;
49 :    
50 :     SET_VECTOR_ELT(ans, 0, allocVector(INTSXP, 2));
51 :     dims = INTEGER(VECTOR_ELT(ans, 0));
52 :     dims[0] = nl1; dims[1] = nl2;
53 :     SET_VECTOR_ELT(ans, 1, allocVector(INTSXP, nl1));
54 :     Memcpy(INTEGER(VECTOR_ELT(ans, 1)), (int *) F->Perm, nl1);
55 :     SET_VECTOR_ELT(ans, 2, allocVector(INTSXP, nl1));
56 :     Memcpy(INTEGER(VECTOR_ELT(ans, 2)), (int *) F->ColCount, nl1);
57 :     Rprintf("Ordering used: %d\n", F->ordering);
58 :     Rprintf("is_super: %d\n", F->is_super);
59 :    
60 :     cholmod_free_sparse(&B, &c);
61 :     cholmod_free_factor(&F, &c);
62 :    
63 :     UNPROTECT(1);
64 :     return ans;
65 :     }
66 :    

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