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 3105 - (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 :     SEXP ans, Givens, Gcpy, nms, pivot, qraux, X;
174 :     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 :     setAttrib(ans, install("useLAPACK"), ScalarLogical(1));
243 : bates 1167 setAttrib(ans, install("rcond"), ScalarReal(rcond));
244 : bates 10 UNPROTECT(2);
245 :     return ans;
246 :     }
247 : bates 1059
248 :     SEXP dense_to_Csparse(SEXP x)
249 :     {
250 : mmaechler 3069 CHM_DN chxd = AS_CHM_xDN(PROTECT(mMatrix_as_geMatrix(x)));
251 : mmaechler 2661 /* cholmod_dense_to_sparse() in CHOLMOD/Core/ below does only work for
252 : maechler 1736 "REAL" 'xtypes', i.e. *not* for "nMatrix".
253 : mmaechler 3092 ===> need "_x" in above AS_CHM_xDN() call.
254 : maechler 1736
255 :     Also it cannot keep symmetric / triangular, hence the
256 :     as_geMatrix() above. Note that this is already a *waste* for
257 :     symmetric matrices; However, we could conceivably use an
258 : mmaechler 2661 enhanced cholmod_dense_to_sparse(), with an extra boolean
259 : maechler 1736 argument for symmetry.
260 : maechler 1725 */
261 : mmaechler 2661 CHM_SP chxs = cholmod_dense_to_sparse(chxd, 1, &c);
262 : maechler 2115 int Rkind = (chxd->xtype == CHOLMOD_REAL) ? Real_KIND2(x) : 0;
263 :     /* Note: when 'x' was integer Matrix, Real_KIND(x) = -1, but *_KIND2(.) = 0 */
264 : maechler 2166 R_CheckStack();
265 : bates 1059
266 : maechler 1960 UNPROTECT(1);
267 : maechler 1725 /* chm_sparse_to_SEXP() *could* deal with symmetric
268 :     * if chxs had such an stype; and we should be able to use uplo below */
269 :     return chm_sparse_to_SEXP(chxs, 1, 0/*TODO: uplo_P(x) if x has an uplo slot*/,
270 :     Rkind, "",
271 : bates 1369 isMatrix(x) ? getAttrib(x, R_DimNamesSymbol)
272 :     : GET_SLOT(x, Matrix_DimNamesSym));
273 : bates 1059 }
274 : bates 1453
275 :    
276 : maechler 2115 SEXP dense_band(SEXP x, SEXP k1P, SEXP k2P)
277 : bates 1453 /* Always returns a full matrix with entries outside the band zeroed
278 : maechler 2115 * Class of the value can be [dln]trMatrix or [dln]geMatrix
279 : bates 1453 */
280 :     {
281 : maechler 2112 int k1 = asInteger(k1P), k2 = asInteger(k2P);
282 : bates 1453
283 : maechler 2112 if (k1 > k2) {
284 : bates 1453 error(_("Lower band %d > upper band %d"), k1, k2);
285 : maechler 2115 return R_NilValue; /* -Wall */
286 : bates 1453 }
287 : maechler 2112 else {
288 : maechler 2115 SEXP ans = PROTECT(dup_mMatrix_as_geMatrix(x));
289 : maechler 2112 int *adims = INTEGER(GET_SLOT(ans, Matrix_DimSym)),
290 : maechler 2115 j, m = adims[0], n = adims[1],
291 : maechler 2112 sqr = (adims[0] == adims[1]),
292 :     tru = (k1 >= 0), trl = (k2 <= 0);
293 : maechler 2115 const char *cl = class_P(ans);
294 : mmaechler 3079 enum dense_enum M_type = ( (cl[0] == 'd') ? ddense :
295 :     ((cl[0] == 'l') ? ldense : ndense));
296 : maechler 2112
297 : maechler 2115
298 :     #define SET_ZERO_OUTSIDE \
299 :     for (j = 0; j < n; j++) { \
300 :     int i, i1 = j - k2, i2 = j + 1 - k1; \
301 : mmaechler 2346 if(i1 > m) i1 = m; \
302 :     if(i2 < 0) i2 = 0; \
303 : maechler 2115 for (i = 0; i < i1; i++) xx[i + j * m] = 0; \
304 :     for (i = i2; i < m; i++) xx[i + j * m] = 0; \
305 : maechler 2112 }
306 : maechler 2115
307 :     if(M_type == ddense) {
308 :     double *xx = REAL(GET_SLOT(ans, Matrix_xSym));
309 :     SET_ZERO_OUTSIDE
310 :     }
311 :     else { /* (M_type == ldense || M_type == ndense) */
312 :     int *xx = LOGICAL(GET_SLOT(ans, Matrix_xSym));
313 :     SET_ZERO_OUTSIDE
314 :     }
315 :    
316 :     if (!sqr || (!tru && !trl)) { /* return the *geMatrix */
317 : maechler 2112 UNPROTECT(1);
318 :     return ans;
319 :     }
320 : maechler 2115 else {
321 :     /* Copy ans to a *trMatrix object (must be square) */
322 :     SEXP aa= PROTECT(NEW_OBJECT(MAKE_CLASS(M_type == ddense? "dtrMatrix":
323 :     (M_type== ldense? "ltrMatrix"
324 :     : "ntrMatrix"))));
325 :     /* Because slots of ans are freshly allocated and ans will not be
326 :     * used, we use the slots themselves and don't duplicate */
327 :     SET_SLOT(aa, Matrix_xSym, GET_SLOT(ans, Matrix_xSym));
328 :     SET_SLOT(aa, Matrix_DimSym, GET_SLOT(ans, Matrix_DimSym));
329 :     SET_SLOT(aa, Matrix_DimNamesSym,GET_SLOT(ans, Matrix_DimNamesSym));
330 :     SET_SLOT(aa, Matrix_diagSym, mkString("N"));
331 :     SET_SLOT(aa, Matrix_uploSym, mkString(tru ? "U" : "L"));
332 :     UNPROTECT(2);
333 :     return aa;
334 :     }
335 : bates 1453 }
336 : maechler 2112 }
337 :    
338 : maechler 2115 SEXP dense_to_symmetric(SEXP x, SEXP uplo, SEXP symm_test)
339 :     /* Class of result will be [dln]syMatrix */
340 : maechler 2112 {
341 : mmaechler 3020 /*== FIXME: allow uplo = NA and then behave a bit like symmpart():
342 :     *== ----- would use the *dimnames* to determine U or L (??)
343 :     */
344 :    
345 : maechler 2115 int symm_tst = asLogical(symm_test);
346 :     SEXP dx = PROTECT(dup_mMatrix_as_geMatrix(x));
347 : mmaechler 3020 SEXP ans, dns, nms_dns;
348 : maechler 2115 const char *cl = class_P(dx);
349 :     /* same as in ..._geMatrix() above:*/
350 : mmaechler 3079 enum dense_enum M_type = ( (cl[0] == 'd') ? ddense :
351 :     ((cl[0] == 'l') ? ldense : ndense));
352 : mmaechler 2817 int *adims = INTEGER(GET_SLOT(dx, Matrix_DimSym)), n = adims[0];
353 :     if(n != adims[1]) {
354 :     UNPROTECT(1);
355 :     error(_("ddense_to_symmetric(): matrix is not square!"));
356 :     return R_NilValue; /* -Wall */
357 :     }
358 : maechler 2112
359 : maechler 2115 if(symm_tst) {
360 : mmaechler 2817 int i,j;
361 : mmaechler 3020 # define CHECK_SYMMETRIC \
362 : maechler 2115 for (j = 0; j < n; j++) \
363 :     for (i = 0; i < j; i++) \
364 :     if(xx[j * n + i] != xx[i * n + j]) { \
365 :     UNPROTECT(1); \
366 :     error(_("matrix is not symmetric [%d,%d]"), i+1, j+1); \
367 :     return R_NilValue; /* -Wall */ \
368 :     }
369 :     if(M_type == ddense) {
370 :    
371 :     double *xx = REAL(GET_SLOT(dx, Matrix_xSym));
372 :     CHECK_SYMMETRIC
373 :     }
374 :     else { /* (M_type == ldense || M_type == ndense) */
375 :    
376 :     int *xx = LOGICAL(GET_SLOT(dx, Matrix_xSym));
377 :     CHECK_SYMMETRIC
378 :     }
379 :     }
380 : mmaechler 3020 # undef CHECK_SYMMETRIC
381 : maechler 2115
382 : mmaechler 3020 ans = PROTECT(NEW_OBJECT(MAKE_CLASS( M_type == ddense ? "dsyMatrix" :
383 :     (M_type == ldense ? "lsyMatrix" :
384 :     "nsyMatrix"))));
385 :    
386 :    
387 :     // --- FIXME: Use MK_SYMMETRIC_DIMNAMES_AND_RETURN from below -- with "uplo" argument
388 :    
389 :     /* need _symmetric_ dimnames */
390 : maechler 2115 dns = GET_SLOT(dx, Matrix_DimNamesSym);
391 : maechler 2112 if(!equal_string_vectors(VECTOR_ELT(dns,0),
392 :     VECTOR_ELT(dns,1))) {
393 :     if(*CHAR(asChar(uplo)) == 'U')
394 :     SET_VECTOR_ELT(dns,0, VECTOR_ELT(dns,1));
395 :     else
396 :     SET_VECTOR_ELT(dns,1, VECTOR_ELT(dns,0));
397 :     }
398 : mmaechler 3020 if(!isNull(nms_dns = getAttrib(dns, R_NamesSymbol)) &&
399 :     !R_compute_identical(STRING_ELT(nms_dns, 0),
400 : mmaechler 3055 STRING_ELT(nms_dns, 1), 16)) { // names(dimnames(.)) :
401 : mmaechler 3020 if(*CHAR(asChar(uplo)) == 'U')
402 :     SET_STRING_ELT(nms_dns, 0, STRING_ELT(nms_dns,1));
403 :     else
404 :     SET_STRING_ELT(nms_dns, 1, STRING_ELT(nms_dns,0));
405 :     setAttrib(dns, R_NamesSymbol, nms_dns);
406 :     }
407 : maechler 2112
408 : maechler 2115 /* Copy dx to ans;
409 :     * Because slots of dx are freshly allocated and dx will not be
410 :     * used, we use the slots themselves and don't duplicate */
411 :     SET_SLOT(ans, Matrix_xSym, GET_SLOT(dx, Matrix_xSym));
412 :     SET_SLOT(ans, Matrix_DimSym, GET_SLOT(dx, Matrix_DimSym));
413 : maechler 2112 SET_SLOT(ans, Matrix_DimNamesSym, dns);
414 :     SET_SLOT(ans, Matrix_uploSym, ScalarString(asChar(uplo)));
415 :    
416 : bates 1453 UNPROTECT(2);
417 : maechler 2112 return ans;
418 : bates 1453 }
419 : maechler 2112
420 :     SEXP ddense_symmpart(SEXP x)
421 :     /* Class of the value will be dsyMatrix */
422 :     {
423 : mmaechler 3020 SEXP dx = dup_mMatrix_as_dgeMatrix(x);
424 : maechler 2115 int *adims = INTEGER(GET_SLOT(dx, Matrix_DimSym)), n = adims[0];
425 : maechler 2112
426 :     if(n != adims[1]) {
427 :     error(_("matrix is not square! (symmetric part)"));
428 : maechler 2115 return R_NilValue; /* -Wall */
429 : maechler 2112 } else {
430 : mmaechler 3020 PROTECT(dx);
431 :     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dsyMatrix"))), dns, nms_dns;
432 : maechler 2115 double *xx = REAL(GET_SLOT(dx, Matrix_xSym));
433 : maechler 2112
434 :     /* only need to assign the *upper* triangle (uplo = "U");
435 :     * noting that diagonal remains unchanged */
436 : mmaechler 3020 for (int j = 0; j < n; j++) {
437 :     for (int i = 0; i < j; i++) {
438 : maechler 2115 xx[j * n + i] = (xx[j * n + i] + xx[i * n + j]) / 2.;
439 : maechler 2112 }
440 :     }
441 :    
442 : mmaechler 3055 // FIXME?: Compare and synchronize with symmetric_DimNames() in ./Mutils.c
443 : mmaechler 3020 # define MK_SYMMETRIC_DIMNAMES_AND_RETURN \
444 :     \
445 :     dns = GET_SLOT(dx, Matrix_DimNamesSym); \
446 :     int J = 1; \
447 :     if(!equal_string_vectors(VECTOR_ELT(dns,0), \
448 :     VECTOR_ELT(dns,1))) { \
449 :     /* _symmetric_ dimnames: behave as symmDN(*, col=TRUE) */ \
450 :     if(isNull(VECTOR_ELT(dns, J))) \
451 :     J = !J; \
452 :     SET_VECTOR_ELT(dns, !J, VECTOR_ELT(dns, J)); \
453 : mmaechler 3055 } \
454 :     /* names(dimnames(.)): */ \
455 : mmaechler 3020 if(!isNull(nms_dns = getAttrib(dns, R_NamesSymbol)) && \
456 :     !R_compute_identical(STRING_ELT(nms_dns, 0), \
457 : mmaechler 3055 STRING_ELT(nms_dns, 1), 16)) { \
458 : mmaechler 3020 SET_STRING_ELT(nms_dns, !J, STRING_ELT(nms_dns, J)); \
459 :     setAttrib(dns, R_NamesSymbol, nms_dns); \
460 :     } \
461 :     \
462 :     /* Copy dx to ans; \
463 :     * Because slots of dx are freshly allocated and dx will not be \
464 :     * used, we use the slots themselves and don't duplicate */ \
465 :     \
466 :     SET_SLOT(ans, Matrix_xSym, GET_SLOT(dx, Matrix_xSym)); \
467 :     SET_SLOT(ans, Matrix_DimSym, GET_SLOT(dx, Matrix_DimSym)); \
468 :     SET_SLOT(ans, Matrix_DimNamesSym, dns); \
469 :     SET_SLOT(ans, Matrix_uploSym, mkString("U")); \
470 :     \
471 :     UNPROTECT(2); \
472 :     return ans
473 : maechler 2112
474 : mmaechler 3020 MK_SYMMETRIC_DIMNAMES_AND_RETURN;
475 : maechler 2112 }
476 :     }
477 :    
478 :     SEXP ddense_skewpart(SEXP x)
479 :     /* Class of the value will be dgeMatrix */
480 :     {
481 : mmaechler 3020 SEXP dx = dup_mMatrix_as_dgeMatrix(x);
482 : maechler 2115 int *adims = INTEGER(GET_SLOT(dx, Matrix_DimSym)), n = adims[0];
483 : maechler 2112
484 :     if(n != adims[1]) {
485 :     error(_("matrix is not square! (skew-symmetric part)"));
486 : maechler 2115 return R_NilValue; /* -Wall */
487 : maechler 2112 } else {
488 : mmaechler 3020 PROTECT(dx);
489 :     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix"))), dns, nms_dns;
490 : maechler 2115 double *xx = REAL(GET_SLOT(dx, Matrix_xSym));
491 : maechler 2112
492 : mmaechler 3020 for (int j = 0; j < n; j++) {
493 : maechler 2115 xx[j * n + j] = 0.;
494 : mmaechler 3020 for (int i = 0; i < j; i++) {
495 : maechler 2115 double s = (xx[j * n + i] - xx[i * n + j]) / 2.;
496 :     xx[j * n + i] = s;
497 :     xx[i * n + j] = -s;
498 : maechler 2112 }
499 :     }
500 :    
501 : mmaechler 3020 MK_SYMMETRIC_DIMNAMES_AND_RETURN;
502 : maechler 2112 }
503 :     }

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