SCM

SCM Repository

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

Annotation of /pkg/Matrix/src/dense.c

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : bates 10 #include "dense.h"
2 : maechler 1725 #include "Mutils.h"
3 : bates 1059 #include "chm_common.h"
4 : bates 10
5 : maechler 1548 /**
6 : bates 10 * Perform a left cyclic shift of columns j to k in the upper triangular
7 :     * matrix x, then restore it to upper triangular form with Givens rotations.
8 :     * The algorithm is based on the Fortran routine DCHEX from Linpack.
9 :     *
10 :     * The lower triangle of x is not modified.
11 : maechler 1548 *
12 : bates 10 * @param x Matrix stored in column-major order
13 :     * @param ldx leading dimension of x
14 :     * @param j column number (0-based) that will be shifted to position k
15 :     * @param k last column number (0-based) to be shifted
16 :     * @param cosines cosines of the Givens rotations
17 :     * @param sines sines of the Givens rotations
18 : maechler 1548 *
19 : bates 10 * @return 0 for success
20 :     */
21 :     static
22 :     int left_cyclic(double x[], int ldx, int j, int k,
23 :     double cosines[], double sines[])
24 :     {
25 :     double *lastcol;
26 :     int i, jj;
27 :    
28 :     if (j >= k)
29 : bates 582 error(_("incorrect left cyclic shift, j (%d) >= k (%d)"), j, k);
30 : bates 10 if (j < 0)
31 : bates 582 error(_("incorrect left cyclic shift, j (%d) < 0"), j, k);
32 : bates 10 if (ldx < k)
33 : bates 582 error(_("incorrect left cyclic shift, k (%d) > ldx (%d)"), k, ldx);
34 : bates 10 lastcol = (double*) R_alloc(k+1, sizeof(double));
35 :     /* keep a copy of column j */
36 :     for(i = 0; i <= j; i++) lastcol[i] = x[i + j*ldx];
37 :     /* For safety, zero the rest */
38 :     for(i = j+1; i <= k; i++) lastcol[i] = 0.;
39 :     for(jj = j+1; jj <= k; jj++) { /* columns to be shifted */
40 :     int diagind = jj*(ldx+1), ind = (jj-j) - 1;
41 :     double tmp = x[diagind], cc, ss;
42 :     /* Calculate the Givens rotation. */
43 :     /* This modified the super-diagonal element */
44 :     F77_CALL(drotg)(x + diagind-1, &tmp, cosines + ind, sines + ind);
45 :     cc = cosines[ind]; ss = sines[ind];
46 :     /* Copy column jj+1 to column jj. */
47 :     for(i = 0; i < jj; i++) x[i + (jj-1)*ldx] = x[i+jj*ldx];
48 :     /* Apply rotation to columns up to k */
49 :     for(i = jj; i < k; i++) {
50 :     tmp = cc*x[(jj-1)+i*ldx] + ss*x[jj+i*ldx];
51 :     x[jj+i*ldx] = cc*x[jj+i*ldx] - ss*x[(jj-1)+i*ldx];
52 :     x[(jj-1)+i*ldx] = tmp;
53 :     }
54 :     /* Apply rotation to lastcol */
55 :     lastcol[jj] = -ss*lastcol[jj-1]; lastcol[jj-1] *= cc;
56 :     }
57 :     /* Copy lastcol to column k */
58 :     for(i = 0; i <= k; i++) x[i+k*ldx] = lastcol[i];
59 :     return 0;
60 :     }
61 :    
62 :     static
63 :     SEXP getGivens(double x[], int ldx, int jmin, int rank)
64 :     {
65 :     int shiftlen = (rank - jmin) - 1;
66 :     SEXP ans = PROTECT(allocVector(VECSXP, 4)), nms, cosines, sines;
67 :    
68 :     SET_VECTOR_ELT(ans, 0, ScalarInteger(jmin));
69 :     SET_VECTOR_ELT(ans, 1, ScalarInteger(rank));
70 :     SET_VECTOR_ELT(ans, 2, cosines = allocVector(REALSXP, shiftlen));
71 :     SET_VECTOR_ELT(ans, 3, sines = allocVector(REALSXP, shiftlen));
72 :     setAttrib(ans, R_NamesSymbol, nms = allocVector(STRSXP, 4));
73 :     SET_STRING_ELT(nms, 0, mkChar("jmin"));
74 :     SET_STRING_ELT(nms, 1, mkChar("rank"));
75 :     SET_STRING_ELT(nms, 2, mkChar("cosines"));
76 :     SET_STRING_ELT(nms, 3, mkChar("sines"));
77 :     if (left_cyclic(x, ldx, jmin, rank - 1, REAL(cosines), REAL(sines)))
78 : bates 582 error(_("Unknown error in getGivens"));
79 : bates 10 UNPROTECT(1);
80 :     return ans;
81 :     }
82 :    
83 :     SEXP checkGivens(SEXP X, SEXP jmin, SEXP rank)
84 :     {
85 :     SEXP ans = PROTECT(allocVector(VECSXP, 2)),
86 :     Xcp = PROTECT(duplicate(X));
87 :     int *Xdims;
88 :    
89 :     if (!(isReal(X) & isMatrix(X)))
90 : bates 582 error(_("X must be a numeric (double precision) matrix"));
91 : bates 10 Xdims = INTEGER(coerceVector(getAttrib(X, R_DimSymbol), INTSXP));
92 :     SET_VECTOR_ELT(ans, 1, getGivens(REAL(Xcp), Xdims[0],
93 :     asInteger(jmin), asInteger(rank)));
94 :     SET_VECTOR_ELT(ans, 0, Xcp);
95 :     UNPROTECT(2);
96 :     return ans;
97 :     }
98 :    
99 :     SEXP lsq_dense_Chol(SEXP X, SEXP y)
100 :     {
101 :     SEXP ans;
102 :     int info, n, p, k, *Xdims, *ydims;
103 :     double *xpx, d_one = 1., d_zero = 0.;
104 :    
105 :     if (!(isReal(X) & isMatrix(X)))
106 : bates 582 error(_("X must be a numeric (double precision) matrix"));
107 : bates 10 Xdims = INTEGER(coerceVector(getAttrib(X, R_DimSymbol), INTSXP));
108 :     n = Xdims[0];
109 :     p = Xdims[1];
110 :     if (!(isReal(y) & isMatrix(y)))
111 : bates 582 error(_("y must be a numeric (double precision) matrix"));
112 : bates 10 ydims = INTEGER(coerceVector(getAttrib(y, R_DimSymbol), INTSXP));
113 :     if (ydims[0] != n)
114 : bates 582 error(_(
115 :     "number of rows in y (%d) does not match number of rows in X (%d)"),
116 : bates 10 ydims[0], n);
117 :     k = ydims[1];
118 :     if (k < 1 || p < 1) return allocMatrix(REALSXP, p, k);
119 :     ans = PROTECT(allocMatrix(REALSXP, p, k));
120 :     F77_CALL(dgemm)("T", "N", &p, &k, &n, &d_one, REAL(X), &n, REAL(y), &n,
121 :     &d_zero, REAL(ans), &p);
122 :     xpx = (double *) R_alloc(p * p, sizeof(double));
123 :     F77_CALL(dsyrk)("U", "T", &p, &n, &d_one, REAL(X), &n, &d_zero,
124 :     xpx, &p);
125 :     F77_CALL(dposv)("U", &p, &k, xpx, &p, REAL(ans), &p, &info);
126 : bates 582 if (info) error(_("Lapack routine dposv returned error code %d"), info);
127 : bates 10 UNPROTECT(1);
128 :     return ans;
129 :     }
130 : maechler 1548
131 :    
132 : bates 10 SEXP lsq_dense_QR(SEXP X, SEXP y)
133 :     {
134 :     SEXP ans;
135 :     int info, n, p, k, *Xdims, *ydims, lwork;
136 :     double *work, tmp, *xvals;
137 :    
138 :     if (!(isReal(X) & isMatrix(X)))
139 : bates 582 error(_("X must be a numeric (double precision) matrix"));
140 : bates 10 Xdims = INTEGER(coerceVector(getAttrib(X, R_DimSymbol), INTSXP));
141 :     n = Xdims[0];
142 :     p = Xdims[1];
143 :     if (!(isReal(y) & isMatrix(y)))
144 : bates 582 error(_("y must be a numeric (double precision) matrix"));
145 : bates 10 ydims = INTEGER(coerceVector(getAttrib(y, R_DimSymbol), INTSXP));
146 :     if (ydims[0] != n)
147 : bates 582 error(_(
148 :     "number of rows in y (%d) does not match number of rows in X (%d)"),
149 : bates 10 ydims[0], n);
150 :     k = ydims[1];
151 :     if (k < 1 || p < 1) return allocMatrix(REALSXP, p, k);
152 : mmaechler 3105 xvals = (double *) Memcpy(R_alloc(n * p, sizeof(double)), REAL(X), n * p);
153 : bates 10 ans = PROTECT(duplicate(y));
154 :     lwork = -1;
155 :     F77_CALL(dgels)("N", &n, &p, &k, xvals, &n, REAL(ans), &n,
156 :     &tmp, &lwork, &info);
157 :     if (info)
158 : bates 582 error(_("First call to Lapack routine dgels returned error code %d"),
159 : bates 10 info);
160 :     lwork = (int) tmp;
161 :     work = (double *) R_alloc(lwork, sizeof(double));
162 :     F77_CALL(dgels)("N", &n, &p, &k, xvals, &n, REAL(ans), &n,
163 :     work, &lwork, &info);
164 :     if (info)
165 : bates 582 error(_("Second call to Lapack routine dgels returned error code %d"),
166 : bates 10 info);
167 :     UNPROTECT(1);
168 :     return ans;
169 :     }
170 : maechler 1548
171 : bates 10 SEXP lapack_qr(SEXP Xin, SEXP tl)
172 :     {
173 : mmaechler 3213 SEXP ans, Givens, Gcpy, nms, pivot, qraux, X, sym;
174 : bates 10 int i, n, nGivens = 0, p, trsz, *Xdims, rank;
175 :     double rcond = 0., tol = asReal(tl), *work;
176 :    
177 :     if (!(isReal(Xin) & isMatrix(Xin)))
178 : bates 582 error(_("X must be a real (numeric) matrix"));
179 :     if (tol < 0.) error(_("tol, given as %g, must be non-negative"), tol);
180 :     if (tol > 1.) error(_("tol, given as %g, must be <= 1"), tol);
181 : bates 10 ans = PROTECT(allocVector(VECSXP,5));
182 :     SET_VECTOR_ELT(ans, 0, X = duplicate(Xin));
183 :     Xdims = INTEGER(coerceVector(getAttrib(X, R_DimSymbol), INTSXP));
184 :     n = Xdims[0]; p = Xdims[1];
185 :     SET_VECTOR_ELT(ans, 2, qraux = allocVector(REALSXP, (n < p) ? n : p));
186 :     SET_VECTOR_ELT(ans, 3, pivot = allocVector(INTSXP, p));
187 :     for (i = 0; i < p; i++) INTEGER(pivot)[i] = i + 1;
188 :     trsz = (n < p) ? n : p; /* size of triangular part of decomposition */
189 :     rank = trsz;
190 :     Givens = PROTECT(allocVector(VECSXP, rank - 1));
191 :     setAttrib(ans, R_NamesSymbol, nms = allocVector(STRSXP, 5));
192 :     SET_STRING_ELT(nms, 0, mkChar("qr"));
193 :     SET_STRING_ELT(nms, 1, mkChar("rank"));
194 :     SET_STRING_ELT(nms, 2, mkChar("qraux"));
195 :     SET_STRING_ELT(nms, 3, mkChar("pivot"));
196 :     SET_STRING_ELT(nms, 4, mkChar("Givens"));
197 :     if (n > 0 && p > 0) {
198 :     int info, *iwork, lwork;
199 :     double *xpt = REAL(X), tmp;
200 :    
201 :     lwork = -1;
202 :     F77_CALL(dgeqrf)(&n, &p, xpt, &n, REAL(qraux), &tmp, &lwork, &info);
203 :     if (info)
204 : bates 582 error(_("First call to dgeqrf returned error code %d"), info);
205 : bates 10 lwork = (int) tmp;
206 :     work = (double *) R_alloc((lwork < 3*trsz) ? 3*trsz : lwork,
207 :     sizeof(double));
208 :     F77_CALL(dgeqrf)(&n, &p, xpt, &n, REAL(qraux), work, &lwork, &info);
209 :     if (info)
210 : bates 582 error(_("Second call to dgeqrf returned error code %d"), info);
211 : bates 357 iwork = (int *) R_alloc(trsz, sizeof(int));
212 : bates 10 F77_CALL(dtrcon)("1", "U", "N", &rank, xpt, &n, &rcond,
213 :     work, iwork, &info);
214 :     if (info)
215 : bates 582 error(_("Lapack routine dtrcon returned error code %d"), info);
216 : bates 10 while (rcond < tol) { /* check diagonal elements */
217 :     double minabs = (xpt[0] < 0.) ? -xpt[0]: xpt[0];
218 :     int jmin = 0;
219 :     for (i = 1; i < rank; i++) {
220 :     double el = xpt[i*(n+1)];
221 :     el = (el < 0.) ? -el: el;
222 :     if (el < minabs) {
223 :     jmin = i;
224 :     minabs = el;
225 :     }
226 :     }
227 :     if (jmin < (rank - 1)) {
228 :     SET_VECTOR_ELT(Givens, nGivens, getGivens(xpt, n, jmin, rank));
229 :     nGivens++;
230 : maechler 1548 }
231 : bates 10 rank--;
232 :     F77_CALL(dtrcon)("1", "U", "N", &rank, xpt, &n, &rcond,
233 :     work, iwork, &info);
234 :     if (info)
235 : bates 582 error(_("Lapack routine dtrcon returned error code %d"), info);
236 : bates 10 }
237 :     }
238 :     SET_VECTOR_ELT(ans, 4, Gcpy = allocVector(VECSXP, nGivens));
239 :     for (i = 0; i < nGivens; i++)
240 :     SET_VECTOR_ELT(Gcpy, i, VECTOR_ELT(Givens, i));
241 :     SET_VECTOR_ELT(ans, 1, ScalarInteger(rank));
242 : mmaechler 3213 sym = PROTECT(install("useLAPACK")); setAttrib(ans, sym, ScalarLogical(1)); UNPROTECT(1);
243 :     sym = PROTECT(install("rcond")); setAttrib(ans, sym, ScalarReal(rcond));UNPROTECT(1);
244 : bates 10 UNPROTECT(2);
245 :     return ans;
246 :     }
247 : bates 1059
248 :     SEXP dense_to_Csparse(SEXP x)
249 :     {
250 : mmaechler 3257 SEXP ge_x = PROTECT(mMatrix_as_geMatrix(x)),
251 :     Dim = GET_SLOT(ge_x, Matrix_DimSym);
252 :     int *dims = INTEGER(Dim);
253 :     Rboolean longi = (dims[0] * (double)dims[1] > INT_MAX);
254 :     // int itype = longi ? CHOLMOD_LONG : CHOLMOD_INT;
255 :     CHM_DN chxd = AS_CHM_xDN(ge_x); // cholmod_dense (has no itype)
256 :     CHM_SP chxs;
257 : mmaechler 2661 /* cholmod_dense_to_sparse() in CHOLMOD/Core/ below does only work for
258 : maechler 1736 "REAL" 'xtypes', i.e. *not* for "nMatrix".
259 : mmaechler 3092 ===> need "_x" in above AS_CHM_xDN() call.
260 : maechler 1736
261 :     Also it cannot keep symmetric / triangular, hence the
262 :     as_geMatrix() above. Note that this is already a *waste* for
263 :     symmetric matrices; However, we could conceivably use an
264 : mmaechler 2661 enhanced cholmod_dense_to_sparse(), with an extra boolean
265 : maechler 1736 argument for symmetry.
266 : maechler 1725 */
267 : mmaechler 3257 #define DLONG
268 :     /* You can try defining DLONG -- then just get a seg.fault :
269 :     * I think it is because of this in ./CHOLMOD/Include/cholmod_core.h :
270 :     *
271 :     The itype of all parameters for all CHOLMOD routines must match.
272 :     --- ^^^^^ ------------------------------------------------------
273 :     but then as_cholmod_dense should *not* make a difference: cholmod_dense has *no* itype (????)
274 :     */
275 :     if(longi) { // calling cholmod_dense_to_sparse() gives wrong matrix
276 :     #ifdef DLONG
277 :     chxs = cholmod_l_dense_to_sparse(chxd, 1, &cl);
278 :     // in gdb, I found that 'chxs' seems "basically empty": all
279 :     // p chxs->foo give ''Cannot access memory at address 0x....''
280 :     // for now rather give error:
281 :     if(cl.status)
282 :     error(_("dense_to_Csparse(<LARGE>): cholmod_l_dense_to_sparse failure status=%d"),
283 :     cl.status);
284 :     #else
285 :     error(_("Matrix dimension %d x %d (= %g) too large [FIXME calling cholmod_l_dense_to_sparse]"),
286 :     m,n, m * (double)n);
287 :     #endif
288 :     } else { // fits, using integer (instead of long int) 'itype'
289 :     chxs = cholmod_dense_to_sparse(chxd, 1, &c);
290 :     }
291 : mmaechler 3167
292 : maechler 2115 int Rkind = (chxd->xtype == CHOLMOD_REAL) ? Real_KIND2(x) : 0;
293 :     /* Note: when 'x' was integer Matrix, Real_KIND(x) = -1, but *_KIND2(.) = 0 */
294 : maechler 2166 R_CheckStack();
295 : bates 1059
296 : maechler 1960 UNPROTECT(1);
297 : maechler 1725 /* chm_sparse_to_SEXP() *could* deal with symmetric
298 :     * if chxs had such an stype; and we should be able to use uplo below */
299 :     return chm_sparse_to_SEXP(chxs, 1, 0/*TODO: uplo_P(x) if x has an uplo slot*/,
300 :     Rkind, "",
301 : bates 1369 isMatrix(x) ? getAttrib(x, R_DimNamesSymbol)
302 :     : GET_SLOT(x, Matrix_DimNamesSym));
303 : bates 1059 }
304 : bates 1453
305 :    
306 : maechler 2115 SEXP dense_band(SEXP x, SEXP k1P, SEXP k2P)
307 : bates 1453 /* Always returns a full matrix with entries outside the band zeroed
308 : maechler 2115 * Class of the value can be [dln]trMatrix or [dln]geMatrix
309 : bates 1453 */
310 :     {
311 : maechler 2112 int k1 = asInteger(k1P), k2 = asInteger(k2P);
312 : bates 1453
313 : maechler 2112 if (k1 > k2) {
314 : bates 1453 error(_("Lower band %d > upper band %d"), k1, k2);
315 : maechler 2115 return R_NilValue; /* -Wall */
316 : bates 1453 }
317 : maechler 2112 else {
318 : maechler 2115 SEXP ans = PROTECT(dup_mMatrix_as_geMatrix(x));
319 : maechler 2112 int *adims = INTEGER(GET_SLOT(ans, Matrix_DimSym)),
320 : maechler 2115 j, m = adims[0], n = adims[1],
321 : maechler 2112 sqr = (adims[0] == adims[1]),
322 :     tru = (k1 >= 0), trl = (k2 <= 0);
323 : maechler 2115 const char *cl = class_P(ans);
324 : mmaechler 3079 enum dense_enum M_type = ( (cl[0] == 'd') ? ddense :
325 :     ((cl[0] == 'l') ? ldense : ndense));
326 : maechler 2112
327 : maechler 2115
328 :     #define SET_ZERO_OUTSIDE \
329 :     for (j = 0; j < n; j++) { \
330 :     int i, i1 = j - k2, i2 = j + 1 - k1; \
331 : mmaechler 2346 if(i1 > m) i1 = m; \
332 :     if(i2 < 0) i2 = 0; \
333 : maechler 2115 for (i = 0; i < i1; i++) xx[i + j * m] = 0; \
334 :     for (i = i2; i < m; i++) xx[i + j * m] = 0; \
335 : maechler 2112 }
336 : maechler 2115
337 :     if(M_type == ddense) {
338 :     double *xx = REAL(GET_SLOT(ans, Matrix_xSym));
339 :     SET_ZERO_OUTSIDE
340 :     }
341 :     else { /* (M_type == ldense || M_type == ndense) */
342 :     int *xx = LOGICAL(GET_SLOT(ans, Matrix_xSym));
343 :     SET_ZERO_OUTSIDE
344 :     }
345 :    
346 :     if (!sqr || (!tru && !trl)) { /* return the *geMatrix */
347 : maechler 2112 UNPROTECT(1);
348 :     return ans;
349 :     }
350 : maechler 2115 else {
351 :     /* Copy ans to a *trMatrix object (must be square) */
352 : mmaechler 3273 SEXP aa= PROTECT(NEW_OBJECT_OF_CLASS(M_type == ddense? "dtrMatrix":
353 :     (M_type== ldense? "ltrMatrix": "ntrMatrix")));
354 : maechler 2115 /* Because slots of ans are freshly allocated and ans will not be
355 :     * used, we use the slots themselves and don't duplicate */
356 :     SET_SLOT(aa, Matrix_xSym, GET_SLOT(ans, Matrix_xSym));
357 :     SET_SLOT(aa, Matrix_DimSym, GET_SLOT(ans, Matrix_DimSym));
358 :     SET_SLOT(aa, Matrix_DimNamesSym,GET_SLOT(ans, Matrix_DimNamesSym));
359 :     SET_SLOT(aa, Matrix_diagSym, mkString("N"));
360 :     SET_SLOT(aa, Matrix_uploSym, mkString(tru ? "U" : "L"));
361 :     UNPROTECT(2);
362 :     return aa;
363 :     }
364 : bates 1453 }
365 : maechler 2112 }
366 :    
367 : maechler 2115 SEXP dense_to_symmetric(SEXP x, SEXP uplo, SEXP symm_test)
368 :     /* Class of result will be [dln]syMatrix */
369 : maechler 2112 {
370 : mmaechler 3020 /*== FIXME: allow uplo = NA and then behave a bit like symmpart():
371 :     *== ----- would use the *dimnames* to determine U or L (??)
372 :     */
373 :    
374 : maechler 2115 int symm_tst = asLogical(symm_test);
375 :     SEXP dx = PROTECT(dup_mMatrix_as_geMatrix(x));
376 : mmaechler 3020 SEXP ans, dns, nms_dns;
377 : maechler 2115 const char *cl = class_P(dx);
378 :     /* same as in ..._geMatrix() above:*/
379 : mmaechler 3079 enum dense_enum M_type = ( (cl[0] == 'd') ? ddense :
380 :     ((cl[0] == 'l') ? ldense : ndense));
381 : mmaechler 2817 int *adims = INTEGER(GET_SLOT(dx, Matrix_DimSym)), n = adims[0];
382 :     if(n != adims[1]) {
383 :     UNPROTECT(1);
384 :     error(_("ddense_to_symmetric(): matrix is not square!"));
385 :     return R_NilValue; /* -Wall */
386 :     }
387 : maechler 2112
388 : maechler 2115 if(symm_tst) {
389 : mmaechler 2817 int i,j;
390 : mmaechler 3020 # define CHECK_SYMMETRIC \
391 : maechler 2115 for (j = 0; j < n; j++) \
392 :     for (i = 0; i < j; i++) \
393 :     if(xx[j * n + i] != xx[i * n + j]) { \
394 :     UNPROTECT(1); \
395 :     error(_("matrix is not symmetric [%d,%d]"), i+1, j+1); \
396 :     return R_NilValue; /* -Wall */ \
397 :     }
398 :     if(M_type == ddense) {
399 :    
400 :     double *xx = REAL(GET_SLOT(dx, Matrix_xSym));
401 :     CHECK_SYMMETRIC
402 :     }
403 :     else { /* (M_type == ldense || M_type == ndense) */
404 :    
405 :     int *xx = LOGICAL(GET_SLOT(dx, Matrix_xSym));
406 :     CHECK_SYMMETRIC
407 :     }
408 :     }
409 : mmaechler 3020 # undef CHECK_SYMMETRIC
410 : maechler 2115
411 : mmaechler 3273 ans = PROTECT(NEW_OBJECT_OF_CLASS(M_type == ddense ? "dsyMatrix" :
412 :     (M_type == ldense ? "lsyMatrix" : "nsyMatrix")));
413 : mmaechler 3020
414 :     // --- FIXME: Use MK_SYMMETRIC_DIMNAMES_AND_RETURN from below -- with "uplo" argument
415 :    
416 :     /* need _symmetric_ dimnames */
417 : maechler 2115 dns = GET_SLOT(dx, Matrix_DimNamesSym);
418 : maechler 2112 if(!equal_string_vectors(VECTOR_ELT(dns,0),
419 :     VECTOR_ELT(dns,1))) {
420 :     if(*CHAR(asChar(uplo)) == 'U')
421 :     SET_VECTOR_ELT(dns,0, VECTOR_ELT(dns,1));
422 :     else
423 :     SET_VECTOR_ELT(dns,1, VECTOR_ELT(dns,0));
424 :     }
425 : mmaechler 3213 nms_dns = PROTECT(getAttrib(dns, R_NamesSymbol));
426 :     if(!isNull(nms_dns) &&
427 : mmaechler 3020 !R_compute_identical(STRING_ELT(nms_dns, 0),
428 : mmaechler 3055 STRING_ELT(nms_dns, 1), 16)) { // names(dimnames(.)) :
429 : mmaechler 3020 if(*CHAR(asChar(uplo)) == 'U')
430 :     SET_STRING_ELT(nms_dns, 0, STRING_ELT(nms_dns,1));
431 :     else
432 :     SET_STRING_ELT(nms_dns, 1, STRING_ELT(nms_dns,0));
433 :     setAttrib(dns, R_NamesSymbol, nms_dns);
434 :     }
435 : maechler 2112
436 : maechler 2115 /* Copy dx to ans;
437 :     * Because slots of dx are freshly allocated and dx will not be
438 :     * used, we use the slots themselves and don't duplicate */
439 :     SET_SLOT(ans, Matrix_xSym, GET_SLOT(dx, Matrix_xSym));
440 :     SET_SLOT(ans, Matrix_DimSym, GET_SLOT(dx, Matrix_DimSym));
441 : maechler 2112 SET_SLOT(ans, Matrix_DimNamesSym, dns);
442 :     SET_SLOT(ans, Matrix_uploSym, ScalarString(asChar(uplo)));
443 :    
444 : mmaechler 3213 UNPROTECT(3);
445 : maechler 2112 return ans;
446 : bates 1453 }
447 : maechler 2112
448 :     SEXP ddense_symmpart(SEXP x)
449 :     /* Class of the value will be dsyMatrix */
450 :     {
451 : mmaechler 3020 SEXP dx = dup_mMatrix_as_dgeMatrix(x);
452 : maechler 2115 int *adims = INTEGER(GET_SLOT(dx, Matrix_DimSym)), n = adims[0];
453 : maechler 2112
454 :     if(n != adims[1]) {
455 :     error(_("matrix is not square! (symmetric part)"));
456 : maechler 2115 return R_NilValue; /* -Wall */
457 : maechler 2112 } else {
458 : mmaechler 3020 PROTECT(dx);
459 : mmaechler 3273 SEXP ans = PROTECT(NEW_OBJECT_OF_CLASS("dsyMatrix")), dns, nms_dns;
460 : maechler 2115 double *xx = REAL(GET_SLOT(dx, Matrix_xSym));
461 : maechler 2112
462 :     /* only need to assign the *upper* triangle (uplo = "U");
463 :     * noting that diagonal remains unchanged */
464 : mmaechler 3020 for (int j = 0; j < n; j++) {
465 :     for (int i = 0; i < j; i++) {
466 : maechler 2115 xx[j * n + i] = (xx[j * n + i] + xx[i * n + j]) / 2.;
467 : maechler 2112 }
468 :     }
469 :    
470 : mmaechler 3055 // FIXME?: Compare and synchronize with symmetric_DimNames() in ./Mutils.c
471 : mmaechler 3020 # define MK_SYMMETRIC_DIMNAMES_AND_RETURN \
472 :     \
473 :     dns = GET_SLOT(dx, Matrix_DimNamesSym); \
474 :     int J = 1; \
475 :     if(!equal_string_vectors(VECTOR_ELT(dns,0), \
476 :     VECTOR_ELT(dns,1))) { \
477 :     /* _symmetric_ dimnames: behave as symmDN(*, col=TRUE) */ \
478 :     if(isNull(VECTOR_ELT(dns, J))) \
479 :     J = !J; \
480 :     SET_VECTOR_ELT(dns, !J, VECTOR_ELT(dns, J)); \
481 : mmaechler 3055 } \
482 :     /* names(dimnames(.)): */ \
483 : mmaechler 3213 nms_dns = PROTECT(getAttrib(dns, R_NamesSymbol)); \
484 :     if(!isNull(nms_dns) && \
485 : mmaechler 3020 !R_compute_identical(STRING_ELT(nms_dns, 0), \
486 : mmaechler 3055 STRING_ELT(nms_dns, 1), 16)) { \
487 : mmaechler 3020 SET_STRING_ELT(nms_dns, !J, STRING_ELT(nms_dns, J)); \
488 :     setAttrib(dns, R_NamesSymbol, nms_dns); \
489 :     } \
490 :     \
491 :     /* Copy dx to ans; \
492 :     * Because slots of dx are freshly allocated and dx will not be \
493 :     * used, we use the slots themselves and don't duplicate */ \
494 :     \
495 :     SET_SLOT(ans, Matrix_xSym, GET_SLOT(dx, Matrix_xSym)); \
496 :     SET_SLOT(ans, Matrix_DimSym, GET_SLOT(dx, Matrix_DimSym)); \
497 :     SET_SLOT(ans, Matrix_DimNamesSym, dns); \
498 :     SET_SLOT(ans, Matrix_uploSym, mkString("U")); \
499 :     \
500 : mmaechler 3213 UNPROTECT(3); \
501 : mmaechler 3020 return ans
502 : maechler 2112
503 : mmaechler 3020 MK_SYMMETRIC_DIMNAMES_AND_RETURN;
504 : maechler 2112 }
505 :     }
506 :    
507 :     SEXP ddense_skewpart(SEXP x)
508 :     /* Class of the value will be dgeMatrix */
509 :     {
510 : mmaechler 3020 SEXP dx = dup_mMatrix_as_dgeMatrix(x);
511 : maechler 2115 int *adims = INTEGER(GET_SLOT(dx, Matrix_DimSym)), n = adims[0];
512 : maechler 2112
513 :     if(n != adims[1]) {
514 :     error(_("matrix is not square! (skew-symmetric part)"));
515 : maechler 2115 return R_NilValue; /* -Wall */
516 : maechler 2112 } else {
517 : mmaechler 3020 PROTECT(dx);
518 : mmaechler 3273 SEXP ans = PROTECT(NEW_OBJECT_OF_CLASS("dgeMatrix")), dns, nms_dns;
519 : maechler 2115 double *xx = REAL(GET_SLOT(dx, Matrix_xSym));
520 : maechler 2112
521 : mmaechler 3020 for (int j = 0; j < n; j++) {
522 : maechler 2115 xx[j * n + j] = 0.;
523 : mmaechler 3020 for (int i = 0; i < j; i++) {
524 : maechler 2115 double s = (xx[j * n + i] - xx[i * n + j]) / 2.;
525 :     xx[j * n + i] = s;
526 :     xx[i * n + j] = -s;
527 : maechler 2112 }
528 :     }
529 :    
530 : mmaechler 3020 MK_SYMMETRIC_DIMNAMES_AND_RETURN;
531 : maechler 2112 }
532 :     }

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