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 733 - (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 :     *
69 :     * @param side LFT or RGT
70 :     * @param transa TRN or NTR
71 :     * @param m number of rows in B and C
72 :     * @param n number of columns in B and C
73 :     * @param Parent Parent array of A
74 :     * @param bi array of row indices of B
75 :     * @param bp array of column pointers of B
76 :     * @param CIP pointer to the row indices for C
77 :     * @param cp array of column pointers for C
78 :     *
79 :     * @return Pointer to the updated row indices for C. The column
80 :     * pointers may also be overwritten.
81 :     */
82 :     SEXP
83 :     lCholClgCsm(enum CBLAS_SIDE side, enum CBLAS_TRANSPOSE transa, int m,
84 :     int n, const int Parent[], const int bi[], const int bp[],
85 :     SEXP CIP, int cp[])
86 :     {
87 :     int *ci = INTEGER(CIP), cnz, extra, i, j, pari, pos;
88 :    
89 :     extra = 0;
90 :     if (side == LFT) {
91 :     if (transa == TRN) {
92 :     error(_("code not yet written"));
93 :     return R_NilValue;
94 :     } else {
95 :     int *Tci, *Tj, *Ti, ntot;
96 :     for (j = 0; j < n; j++) {
97 :     int ii, ii2 = bp[j + 1];
98 :     for (ii = bp[j]; ii < ii2; ii++)
99 :     for (pari = bi[ii]; pari >= 0; pari = Parent[pari])
100 :     if (check_csc_index(cp, ci, pari, j, -1) < 0) extra++;
101 :     }
102 :     /* This is not quite right. C may contain more elements than in the solution. */
103 :     if (!extra) return CIP;
104 :    
105 :     cnz = cp[n];
106 :     ntot = cnz + extra;
107 :     Ti = Memcpy(Calloc(ntot, int), ci, cnz);
108 :     Tj = expand_cmprPt(n, cp, Calloc(ntot, int));
109 :     Tci = Calloc(ntot, int);
110 :    
111 :     pos = cnz;
112 :     for (j = 0; j < n; j++) {
113 :     int ii, ii2 = bp[j + 1];
114 :     for (ii = bp[j]; ii < ii2; ii++)
115 :     for (pari = bi[ii]; pari >= 0; pari = Parent[pari])
116 :     if (check_csc_index(cp, ci, pari, j, -1) < 0) {
117 :     Ti[pos] = pari;
118 :     Tj[pos] = j;
119 :     pos++;
120 :     }
121 :     }
122 :     triplet_to_col(m, n, ntot, Ti, Tj, (double *) NULL,
123 :     cp, Tci, (double *) NULL);
124 :    
125 :     cnz = cp[n];
126 :     CIP = PROTECT(allocVector(INTSXP, cnz));
127 :     Memcpy(INTEGER(CIP), Tci, cnz);
128 :    
129 :     Free(Tci); Free(Ti); Free(Tj);
130 :     UNPROTECT(1);
131 :     return CIP;
132 :     }
133 :     } else {
134 :     if (transa == TRN) {
135 :     int *Tci, *Tj, *Ti, ntot;
136 :     for (j = 0; j < n; j++) {
137 :     int ii, ii2 = bp[j + 1];
138 :     for (pari = j; pari >= 0; pari = Parent[pari]) {
139 :     for (ii = bp[j]; ii < ii2; ii++)
140 :     if (check_csc_index(cp, ci, ii, pari, -1) < 0) extra++;
141 :     }
142 :     /* This is not quite right. C may contain more elements than in the solution. */
143 :     if (!extra) return CIP;
144 :    
145 :     cnz = cp[n];
146 :     ntot = cnz + extra;
147 :     Ti = Memcpy(Calloc(ntot, int), ci, cnz);
148 :     Tj = expand_cmprPt(n, cp, Calloc(ntot, int));
149 :     Tci = Calloc(ntot, int);
150 :    
151 :     pos = cnz;
152 :     for (j = 0; j < n; j++) {
153 :     int ii, ii2 = bp[j + 1];
154 :     for (pari = j; pari >= 0; pari = Parent[pari]) {
155 :     for (ii = bp[j]; ii < ii2; ii++)
156 :     if (check_csc_index(cp, ci, ii, pari, -1) < 0) {
157 :     Ti[pos] = ii;
158 :     Tj[pos] = pari;
159 :     pos++;
160 :     }
161 :     }
162 :     triplet_to_col(m, n, ntot, Ti, Tj, (double *) NULL,
163 :     cp, Tci, (double *) NULL);
164 :    
165 :     cnz = cp[n];
166 :     CIP = PROTECT(allocVector(INTSXP, cnz));
167 :     Memcpy(INTEGER(CIP), Tci, cnz);
168 :    
169 :     Free(Tci); Free(Ti); Free(Tj);
170 :     UNPROTECT(1);
171 :     return CIP;
172 :     } else {
173 :     error(_("code not yet written"));
174 :     return R_NilValue;
175 :     }
176 :     }
177 :     }
178 :    
179 :     SEXP lCholCMatrix_lgCMatrix_solve(SEXP a, SEXP b)
180 :     {
181 :     SEXP ans = PROTECT(duplicate(b));
182 :     int n = INTEGER(GET_SLOT(a, Matrix_DimSym))[0],
183 :     *bdims = INTEGER(GET_SLOT(b, Matrix_DimSym));
184 :    
185 :     if (n != bdims[0])
186 :     error(_("Dimensions of a (%d,%d) and b (%d,%d) do not conform"),
187 :     n, n, bdims[0], bdims[1]);
188 :     SET_SLOT(ans, Matrix_iSym,
189 :     lCholClgCsm(LFT, NTR, bdims[0], bdims[1],
190 :     INTEGER(GET_SLOT(a, Matrix_ParentSym)),
191 :     INTEGER(GET_SLOT(b, Matrix_iSym)),
192 :     INTEGER(GET_SLOT(b, Matrix_pSym)),
193 :     GET_SLOT(ans, Matrix_iSym),
194 :     INTEGER(GET_SLOT(ans, Matrix_pSym))));
195 :     UNPROTECT(1);
196 :     return ans;
197 :     }
198 :    
199 :    

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