SCM

SCM Repository

[matrix] Annotation of /branches/Matrix-mer2/src/lCholCMatrix.c
ViewVC logotype

Annotation of /branches/Matrix-mer2/src/lCholCMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : bates 733 /* LDL' factorization of a logical,
2 :     * sparse column-oriented matrix */
3 :     #include "lCholCMatrix.h"
4 :    
5 :     SEXP lCholCMatrix_validate(SEXP x)
6 :     {
7 :     int i, n = INTEGER(GET_SLOT(x, Matrix_DimSym))[0];
8 :     SEXP perm = GET_SLOT(x, Matrix_permSym),
9 :     Parent = GET_SLOT(x, Matrix_ParentSym);
10 :    
11 :     if (length(perm) != n)
12 :     return mkString(_("slot perm must have length n"));
13 :     if (length(Parent) != n)
14 :     return mkString(_("slot Parent must have length n"));
15 :     if (!R_ldl_valid_perm(n, INTEGER(perm)))
16 :     return mkString(_("slot perm is not a valid 0-based permutation"));
17 :     for (i = 0; i < n; i++) {
18 :     int pari = INTEGER(Parent)[i];
19 :     if (pari < -1 || pari > (n - 1))
20 :     return mkString(_("an element of the Parent array is not in range [-1,n-1]"));
21 :     }
22 :    
23 :     return ScalarLogical(1);
24 :     }
25 :    
26 :     /**
27 :     * Create the structure of the inverse of L from the LDL' factorization.
28 :     *
29 :     * @param x Pointer to a lCholCMatrix object
30 :     *
31 :     * @return An ltCMatrix object representing L^{-1}
32 :     */
33 :     SEXP lCholCMatrix_solve(SEXP x)
34 :     {
35 :     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("ltCMatrix")));
36 :     SEXP Parent = GET_SLOT(x, Matrix_ParentSym);
37 :     int i, n = length(Parent), pari, pos;
38 :     int *ai, *ap, *nzc = Calloc(n, int), ntot;
39 :    
40 :     ntot = 0;
41 :     for (i = n - 1; i >= 0; i--) { /* count non-zeros in each column */
42 :     pari = INTEGER(Parent)[i];
43 :     ntot += (nzc[i] = (pari >= 0) ? 1 + nzc[pari] : 0);
44 :     }
45 :    
46 :     slot_dup(ans, x, Matrix_DimSym);
47 :     slot_dup(ans, x, Matrix_DimNamesSym);
48 :     slot_dup(ans, x, Matrix_uploSym);
49 :     slot_dup(ans, x, Matrix_diagSym);
50 :     ai = INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, ntot));
51 :     ap = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, n + 1));
52 :     ap[0] = pos = 0;
53 :     for (i = 0; i < n; i++) {
54 :     ap[i + 1] = ap[i] + nzc[i];
55 :     for (pari = INTEGER(Parent)[i]; pari >= 0; pari = INTEGER(Parent)[pari])
56 :     ai[pos++] = pari;
57 :     }
58 :    
59 :     Free(nzc);
60 :     UNPROTECT(1);
61 :     return ans;
62 :     }
63 :    
64 :     /**
65 :     * Solve one of the matrix equations op(A)*C = B, or * C*op(A) = B
66 :     * where A is an lCholCMatrix object and B and C are lgCMatrix
67 :     * objects.
68 : bates 746 *
69 :     * An early check is done to see if A is the identity, in which case
70 :     * BIP is returned and bp is unmodified.
71 : bates 733 *
72 :     * @param side LFT or RGT
73 :     * @param transa TRN or NTR
74 :     * @param m number of rows in B and C
75 :     * @param n number of columns in B and C
76 :     * @param Parent Parent array of A
77 : bates 746 * @param BIP pointer to the row indices for B
78 :     * @param bp array of column pointers for B (may be overwritten)
79 : bates 733 *
80 :     * @return Pointer to the updated row indices for C. The column
81 : bates 746 * pointers bp are overwritten with cp.
82 : bates 733 */
83 :     SEXP
84 :     lCholClgCsm(enum CBLAS_SIDE side, enum CBLAS_TRANSPOSE transa, int m,
85 : bates 746 int n, const int Parent[], SEXP BIP, int bp[])
86 : bates 733 {
87 : bates 746 int *bi = INTEGER(BIP), bnz, extra, ident, j, pari, pos;
88 :     int nca = (transa == TRN) ? n : m;
89 : bates 733
90 : bates 746 ident = 1;
91 :     for (j = 0; j < nca; j++)
92 :     if (Parent[j] >= 0) {
93 :     ident = 0;
94 :     break;
95 :     }
96 :     if (ident) return BIP;
97 :    
98 : bates 733 extra = 0;
99 :     if (side == LFT) {
100 :     if (transa == TRN) {
101 :     error(_("code not yet written"));
102 :     return R_NilValue;
103 :     } else {
104 :     int *Tci, *Tj, *Ti, ntot;
105 :     for (j = 0; j < n; j++) {
106 :     int ii, ii2 = bp[j + 1];
107 :     for (ii = bp[j]; ii < ii2; ii++)
108 : bates 746 for (pari = Parent[bi[ii]]; pari >= 0;
109 :     pari = Parent[pari]) extra++;
110 : bates 733 }
111 :    
112 : bates 746 bnz = bp[n];
113 :     ntot = bnz + extra;
114 :     Ti = Memcpy(Calloc(ntot, int), bi, bnz);
115 :     Tj = expand_cmprPt(n, bp, Calloc(ntot, int));
116 : bates 733 Tci = Calloc(ntot, int);
117 :    
118 : bates 746 pos = bnz;
119 : bates 733 for (j = 0; j < n; j++) {
120 :     int ii, ii2 = bp[j + 1];
121 :     for (ii = bp[j]; ii < ii2; ii++)
122 : bates 746 for (pari = Parent[bi[ii]]; pari >= 0;
123 :     pari = Parent[pari]) {
124 : bates 733 Ti[pos] = pari;
125 :     Tj[pos] = j;
126 :     pos++;
127 :     }
128 :     }
129 :     triplet_to_col(m, n, ntot, Ti, Tj, (double *) NULL,
130 : bates 746 bp, Tci, (double *) NULL);
131 : bates 733
132 : bates 746 bnz = bp[n];
133 :     BIP = PROTECT(allocVector(INTSXP, bnz));
134 :     Memcpy(INTEGER(BIP), Tci, bnz);
135 : bates 733
136 :     Free(Tci); Free(Ti); Free(Tj);
137 :     UNPROTECT(1);
138 : bates 746 return BIP;
139 : bates 733 }
140 :     } else {
141 :     if (transa == TRN) {
142 :     int *Tci, *Tj, *Ti, ntot;
143 :     for (j = 0; j < n; j++) {
144 :     int ii, ii2 = bp[j + 1];
145 : bates 746 for (pari = Parent[j]; pari >= 0; pari = Parent[pari])
146 :     for (ii = bp[j]; ii < ii2; ii++) extra++;
147 : bates 733 }
148 :    
149 : bates 746 bnz = bp[n];
150 :     ntot = bnz + extra;
151 :     Ti = Memcpy(Calloc(ntot, int), bi, bnz);
152 :     Tj = expand_cmprPt(n, bp, Calloc(ntot, int));
153 : bates 733 Tci = Calloc(ntot, int);
154 :    
155 : bates 746 pos = bnz;
156 : bates 733 for (j = 0; j < n; j++) {
157 :     int ii, ii2 = bp[j + 1];
158 : bates 746 for (pari = Parent[j]; pari >= 0; pari = Parent[pari]) {
159 :     for (ii = bp[j]; ii < ii2; ii++) {
160 :     Ti[pos] = bi[ii];
161 :     Tj[pos] = pari;
162 :     pos++;
163 :     }
164 :     }
165 : bates 733 }
166 :     triplet_to_col(m, n, ntot, Ti, Tj, (double *) NULL,
167 : bates 746 bp, Tci, (double *) NULL);
168 :     bnz = bp[n];
169 :     BIP = PROTECT(allocVector(INTSXP, bnz));
170 :     Memcpy(INTEGER(BIP), Tci, bnz);
171 : bates 733
172 :     Free(Tci); Free(Ti); Free(Tj);
173 :     UNPROTECT(1);
174 : bates 746 return BIP;
175 : bates 733 } else {
176 :     error(_("code not yet written"));
177 :     return R_NilValue;
178 :     }
179 :     }
180 :     }
181 :    
182 :     SEXP lCholCMatrix_lgCMatrix_solve(SEXP a, SEXP b)
183 :     {
184 :     SEXP ans = PROTECT(duplicate(b));
185 :     int n = INTEGER(GET_SLOT(a, Matrix_DimSym))[0],
186 :     *bdims = INTEGER(GET_SLOT(b, Matrix_DimSym));
187 :    
188 :     if (n != bdims[0])
189 :     error(_("Dimensions of a (%d,%d) and b (%d,%d) do not conform"),
190 :     n, n, bdims[0], bdims[1]);
191 :     SET_SLOT(ans, Matrix_iSym,
192 :     lCholClgCsm(LFT, NTR, bdims[0], bdims[1],
193 :     INTEGER(GET_SLOT(a, Matrix_ParentSym)),
194 :     GET_SLOT(ans, Matrix_iSym),
195 :     INTEGER(GET_SLOT(ans, Matrix_pSym))));
196 :     UNPROTECT(1);
197 :     return ans;
198 :     }
199 :    
200 :    

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