SCM

SCM Repository

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

Annotation of /pkg/src/dense.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2166 - (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 :     xvals = (double *) R_alloc(n * p, sizeof(double));
153 :     Memcpy(xvals, REAL(X), n * p);
154 :     ans = PROTECT(duplicate(y));
155 :     lwork = -1;
156 :     F77_CALL(dgels)("N", &n, &p, &k, xvals, &n, REAL(ans), &n,
157 :     &tmp, &lwork, &info);
158 :     if (info)
159 : bates 582 error(_("First call to Lapack routine dgels returned error code %d"),
160 : bates 10 info);
161 :     lwork = (int) tmp;
162 :     work = (double *) R_alloc(lwork, sizeof(double));
163 :     F77_CALL(dgels)("N", &n, &p, &k, xvals, &n, REAL(ans), &n,
164 :     work, &lwork, &info);
165 :     if (info)
166 : bates 582 error(_("Second call to Lapack routine dgels returned error code %d"),
167 : bates 10 info);
168 :     UNPROTECT(1);
169 :     return ans;
170 :     }
171 : maechler 1548
172 : bates 10 SEXP lapack_qr(SEXP Xin, SEXP tl)
173 :     {
174 :     SEXP ans, Givens, Gcpy, nms, pivot, qraux, X;
175 :     int i, n, nGivens = 0, p, trsz, *Xdims, rank;
176 :     double rcond = 0., tol = asReal(tl), *work;
177 :    
178 :     if (!(isReal(Xin) & isMatrix(Xin)))
179 : bates 582 error(_("X must be a real (numeric) matrix"));
180 :     if (tol < 0.) error(_("tol, given as %g, must be non-negative"), tol);
181 :     if (tol > 1.) error(_("tol, given as %g, must be <= 1"), tol);
182 : bates 10 ans = PROTECT(allocVector(VECSXP,5));
183 :     SET_VECTOR_ELT(ans, 0, X = duplicate(Xin));
184 :     Xdims = INTEGER(coerceVector(getAttrib(X, R_DimSymbol), INTSXP));
185 :     n = Xdims[0]; p = Xdims[1];
186 :     SET_VECTOR_ELT(ans, 2, qraux = allocVector(REALSXP, (n < p) ? n : p));
187 :     SET_VECTOR_ELT(ans, 3, pivot = allocVector(INTSXP, p));
188 :     for (i = 0; i < p; i++) INTEGER(pivot)[i] = i + 1;
189 :     trsz = (n < p) ? n : p; /* size of triangular part of decomposition */
190 :     rank = trsz;
191 :     Givens = PROTECT(allocVector(VECSXP, rank - 1));
192 :     setAttrib(ans, R_NamesSymbol, nms = allocVector(STRSXP, 5));
193 :     SET_STRING_ELT(nms, 0, mkChar("qr"));
194 :     SET_STRING_ELT(nms, 1, mkChar("rank"));
195 :     SET_STRING_ELT(nms, 2, mkChar("qraux"));
196 :     SET_STRING_ELT(nms, 3, mkChar("pivot"));
197 :     SET_STRING_ELT(nms, 4, mkChar("Givens"));
198 :     if (n > 0 && p > 0) {
199 :     int info, *iwork, lwork;
200 :     double *xpt = REAL(X), tmp;
201 :    
202 :     lwork = -1;
203 :     F77_CALL(dgeqrf)(&n, &p, xpt, &n, REAL(qraux), &tmp, &lwork, &info);
204 :     if (info)
205 : bates 582 error(_("First call to dgeqrf returned error code %d"), info);
206 : bates 10 lwork = (int) tmp;
207 :     work = (double *) R_alloc((lwork < 3*trsz) ? 3*trsz : lwork,
208 :     sizeof(double));
209 :     F77_CALL(dgeqrf)(&n, &p, xpt, &n, REAL(qraux), work, &lwork, &info);
210 :     if (info)
211 : bates 582 error(_("Second call to dgeqrf returned error code %d"), info);
212 : bates 357 iwork = (int *) R_alloc(trsz, sizeof(int));
213 : bates 10 F77_CALL(dtrcon)("1", "U", "N", &rank, xpt, &n, &rcond,
214 :     work, iwork, &info);
215 :     if (info)
216 : bates 582 error(_("Lapack routine dtrcon returned error code %d"), info);
217 : bates 10 while (rcond < tol) { /* check diagonal elements */
218 :     double minabs = (xpt[0] < 0.) ? -xpt[0]: xpt[0];
219 :     int jmin = 0;
220 :     for (i = 1; i < rank; i++) {
221 :     double el = xpt[i*(n+1)];
222 :     el = (el < 0.) ? -el: el;
223 :     if (el < minabs) {
224 :     jmin = i;
225 :     minabs = el;
226 :     }
227 :     }
228 :     if (jmin < (rank - 1)) {
229 :     SET_VECTOR_ELT(Givens, nGivens, getGivens(xpt, n, jmin, rank));
230 :     nGivens++;
231 : maechler 1548 }
232 : bates 10 rank--;
233 :     F77_CALL(dtrcon)("1", "U", "N", &rank, xpt, &n, &rcond,
234 :     work, iwork, &info);
235 :     if (info)
236 : bates 582 error(_("Lapack routine dtrcon returned error code %d"), info);
237 : bates 10 }
238 :     }
239 :     SET_VECTOR_ELT(ans, 4, Gcpy = allocVector(VECSXP, nGivens));
240 :     for (i = 0; i < nGivens; i++)
241 :     SET_VECTOR_ELT(Gcpy, i, VECTOR_ELT(Givens, i));
242 :     SET_VECTOR_ELT(ans, 1, ScalarInteger(rank));
243 :     setAttrib(ans, install("useLAPACK"), ScalarLogical(1));
244 : bates 1167 setAttrib(ans, install("rcond"), ScalarReal(rcond));
245 : bates 10 UNPROTECT(2);
246 :     return ans;
247 :     }
248 : bates 1059
249 :     SEXP dense_to_Csparse(SEXP x)
250 :     {
251 : maechler 1960 CHM_DN chxd = AS_CHM_DN(PROTECT(mMatrix_as_geMatrix(x)));
252 : maechler 1736 /* cholmod_dense_to_sparse() in CHOLMOD/Core/ below does only work for
253 :     "REAL" 'xtypes', i.e. *not* for "nMatrix".
254 :     ===> need "_x" in above call.
255 :    
256 :     Also it cannot keep symmetric / triangular, hence the
257 :     as_geMatrix() above. Note that this is already a *waste* for
258 :     symmetric matrices; However, we could conceivably use an
259 :     enhanced cholmod_dense_to_sparse(), with an extra boolean
260 :     argument for symmetry.
261 : maechler 1725 */
262 : maechler 1960 CHM_SP chxs = cholmod_dense_to_sparse(chxd, 1, &c);
263 : maechler 2115 int Rkind = (chxd->xtype == CHOLMOD_REAL) ? Real_KIND2(x) : 0;
264 :     /* Note: when 'x' was integer Matrix, Real_KIND(x) = -1, but *_KIND2(.) = 0 */
265 : maechler 2166 R_CheckStack();
266 : bates 1059
267 : maechler 1960 UNPROTECT(1);
268 : maechler 1725 /* chm_sparse_to_SEXP() *could* deal with symmetric
269 :     * if chxs had such an stype; and we should be able to use uplo below */
270 :     return chm_sparse_to_SEXP(chxs, 1, 0/*TODO: uplo_P(x) if x has an uplo slot*/,
271 :     Rkind, "",
272 : bates 1369 isMatrix(x) ? getAttrib(x, R_DimNamesSymbol)
273 :     : GET_SLOT(x, Matrix_DimNamesSym));
274 : bates 1059 }
275 : bates 1453
276 :    
277 : maechler 2115 SEXP dense_band(SEXP x, SEXP k1P, SEXP k2P)
278 : bates 1453 /* Always returns a full matrix with entries outside the band zeroed
279 : maechler 2115 * Class of the value can be [dln]trMatrix or [dln]geMatrix
280 : bates 1453 */
281 :     {
282 : maechler 2112 int k1 = asInteger(k1P), k2 = asInteger(k2P);
283 : bates 1453
284 : maechler 2112 if (k1 > k2) {
285 : bates 1453 error(_("Lower band %d > upper band %d"), k1, k2);
286 : maechler 2115 return R_NilValue; /* -Wall */
287 : bates 1453 }
288 : maechler 2112 else {
289 : maechler 2115 SEXP ans = PROTECT(dup_mMatrix_as_geMatrix(x));
290 : maechler 2112 int *adims = INTEGER(GET_SLOT(ans, Matrix_DimSym)),
291 : maechler 2115 j, m = adims[0], n = adims[1],
292 : maechler 2112 sqr = (adims[0] == adims[1]),
293 :     tru = (k1 >= 0), trl = (k2 <= 0);
294 : maechler 2115 const char *cl = class_P(ans);
295 :     enum dense_enum { ddense, ldense, ndense
296 :     } M_type = ( (cl[0] == 'd') ? ddense :
297 :     ((cl[0] == 'l') ? ldense : ndense));
298 : maechler 2112
299 : maechler 2115
300 :     #define SET_ZERO_OUTSIDE \
301 :     for (j = 0; j < n; j++) { \
302 :     int i, i1 = j - k2, i2 = j + 1 - k1; \
303 :     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 : bates 2114 return R_NilValue; /* -Wall */
337 : maechler 2112 }
338 :    
339 : maechler 2115 SEXP dense_to_symmetric(SEXP x, SEXP uplo, SEXP symm_test)
340 :     /* Class of result will be [dln]syMatrix */
341 : maechler 2112 {
342 : maechler 2115 int symm_tst = asLogical(symm_test);
343 :     SEXP dx = PROTECT(dup_mMatrix_as_geMatrix(x));
344 :     SEXP ans, dns;
345 :     const char *cl = class_P(dx);
346 :     /* same as in ..._geMatrix() above:*/
347 :     enum dense_enum { ddense, ldense, ndense
348 :     } M_type = ( (cl[0] == 'd') ? ddense :
349 :     ((cl[0] == 'l') ? ldense : ndense));
350 : maechler 2112
351 : maechler 2115 if(symm_tst) {
352 :     int *adims = INTEGER(GET_SLOT(dx, Matrix_DimSym)), n = adims[0], i,j;
353 :     if(n != adims[1]) {
354 :     UNPROTECT(1);
355 :     error(_("ddense_to_symmetric(): matrix is not square!"));
356 :     return R_NilValue; /* -Wall */
357 :     }
358 :     #define CHECK_SYMMETRIC \
359 :     for (j = 0; j < n; j++) \
360 :     for (i = 0; i < j; i++) \
361 :     if(xx[j * n + i] != xx[i * n + j]) { \
362 :     UNPROTECT(1); \
363 :     error(_("matrix is not symmetric [%d,%d]"), i+1, j+1); \
364 :     return R_NilValue; /* -Wall */ \
365 :     }
366 :     if(M_type == ddense) {
367 :    
368 :     double *xx = REAL(GET_SLOT(dx, Matrix_xSym));
369 :     CHECK_SYMMETRIC
370 :     }
371 :     else { /* (M_type == ldense || M_type == ndense) */
372 :    
373 :     int *xx = LOGICAL(GET_SLOT(dx, Matrix_xSym));
374 :     CHECK_SYMMETRIC
375 :     }
376 :     }
377 :    
378 :     dns = GET_SLOT(dx, Matrix_DimNamesSym);
379 : maechler 2112 if(!equal_string_vectors(VECTOR_ELT(dns,0),
380 :     VECTOR_ELT(dns,1))) {
381 :     /* need _symmetric_ dimnames */
382 :     if(*CHAR(asChar(uplo)) == 'U')
383 :     SET_VECTOR_ELT(dns,0, VECTOR_ELT(dns,1));
384 :     else
385 :     SET_VECTOR_ELT(dns,1, VECTOR_ELT(dns,0));
386 :     }
387 :    
388 : maechler 2115 /* Copy dx to ans;
389 :     * Because slots of dx are freshly allocated and dx will not be
390 :     * used, we use the slots themselves and don't duplicate */
391 :     ans = PROTECT(NEW_OBJECT(MAKE_CLASS( M_type == ddense ? "dsyMatrix" :
392 :     (M_type == ldense ? "lsyMatrix" :
393 :     "nsyMatrix"))));
394 :     SET_SLOT(ans, Matrix_xSym, GET_SLOT(dx, Matrix_xSym));
395 :     SET_SLOT(ans, Matrix_DimSym, GET_SLOT(dx, Matrix_DimSym));
396 : maechler 2112 SET_SLOT(ans, Matrix_DimNamesSym, dns);
397 :     SET_SLOT(ans, Matrix_uploSym, ScalarString(asChar(uplo)));
398 :    
399 : bates 1453 UNPROTECT(2);
400 : maechler 2112 return ans;
401 : bates 1453 }
402 : maechler 2112
403 :     SEXP ddense_symmpart(SEXP x)
404 :     /* Class of the value will be dsyMatrix */
405 :     {
406 : maechler 2115 SEXP dx = PROTECT(dup_mMatrix_as_dgeMatrix(x));
407 :     int *adims = INTEGER(GET_SLOT(dx, Matrix_DimSym)), n = adims[0];
408 : maechler 2112
409 :     if(n != adims[1]) {
410 :     UNPROTECT(1);
411 :     error(_("matrix is not square! (symmetric part)"));
412 : maechler 2115 return R_NilValue; /* -Wall */
413 : maechler 2112 } else {
414 :     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dsyMatrix"))), dns;
415 : maechler 2115 double *xx = REAL(GET_SLOT(dx, Matrix_xSym));
416 : maechler 2112 int i,j;
417 :    
418 :     /* only need to assign the *upper* triangle (uplo = "U");
419 :     * noting that diagonal remains unchanged */
420 :     for (j = 0; j < n; j++) {
421 :     for (i = 0; i < j; i++) {
422 : maechler 2115 xx[j * n + i] = (xx[j * n + i] + xx[i * n + j]) / 2.;
423 : maechler 2112 }
424 :     }
425 :    
426 : maechler 2115 dns = GET_SLOT(dx, Matrix_DimNamesSym);
427 : maechler 2112 if(!equal_string_vectors(VECTOR_ELT(dns,0),
428 :     VECTOR_ELT(dns,1))) {
429 :     /* need _symmetric_ dimnames */
430 :     SET_VECTOR_ELT(dns,0, VECTOR_ELT(dns,1));/* ==> uplo = "U" */
431 :     }
432 :    
433 : maechler 2115 /* Copy dx to ans;
434 :     * Because slots of dx are freshly allocated and dx will not be
435 :     * used, we use the slots themselves and don't duplicate */
436 :    
437 :     SET_SLOT(ans, Matrix_xSym, GET_SLOT(dx, Matrix_xSym));
438 :     SET_SLOT(ans, Matrix_DimSym, GET_SLOT(dx, Matrix_DimSym));
439 : maechler 2112 SET_SLOT(ans, Matrix_DimNamesSym, dns);
440 :     SET_SLOT(ans, Matrix_uploSym, mkString("U"));
441 :    
442 :     UNPROTECT(2);
443 :     return ans;
444 :     }
445 : bates 2114 return R_NilValue; /* -Wall */
446 : maechler 2112 }
447 :    
448 :     SEXP ddense_skewpart(SEXP x)
449 :     /* Class of the value will be dgeMatrix */
450 :     {
451 : maechler 2115 SEXP dx = PROTECT(dup_mMatrix_as_dgeMatrix(x));
452 :     int *adims = INTEGER(GET_SLOT(dx, Matrix_DimSym)), n = adims[0];
453 : maechler 2112
454 :     if(n != adims[1]) {
455 :     UNPROTECT(1);
456 :     error(_("matrix is not square! (skew-symmetric part)"));
457 : maechler 2115 return R_NilValue; /* -Wall */
458 : maechler 2112 } else {
459 :     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix"))), dns;
460 : maechler 2115 double *xx = REAL(GET_SLOT(dx, Matrix_xSym));
461 : maechler 2112 int i,j;
462 :    
463 :     for (j = 0; j < n; j++) {
464 : maechler 2115 xx[j * n + j] = 0.;
465 : maechler 2112 for (i = 0; i < j; i++) {
466 : maechler 2115 double s = (xx[j * n + i] - xx[i * n + j]) / 2.;
467 :     xx[j * n + i] = s;
468 :     xx[i * n + j] = -s;
469 : maechler 2112 }
470 :     }
471 :    
472 : maechler 2115 dns = GET_SLOT(dx, Matrix_DimNamesSym);
473 : maechler 2112 if(!equal_string_vectors(VECTOR_ELT(dns,0),
474 :     VECTOR_ELT(dns,1))) {
475 :     /* need _symmetric_ dimnames */
476 :     SET_VECTOR_ELT(dns,0, VECTOR_ELT(dns,1));/* uplo = "U" */
477 :     }
478 :    
479 : maechler 2115 /* Copy dx to ans;
480 :     * Because slots of dx are freshly allocated and dx will not be
481 :     * used, we use the slots themselves and don't duplicate */
482 :    
483 :     SET_SLOT(ans, Matrix_xSym, GET_SLOT(dx, Matrix_xSym));
484 :     SET_SLOT(ans, Matrix_DimSym, GET_SLOT(dx, Matrix_DimSym));
485 : maechler 2112 SET_SLOT(ans, Matrix_DimNamesSym, dns);
486 :     SET_SLOT(ans, Matrix_uploSym, mkString("U"));
487 :    
488 :     UNPROTECT(2);
489 :     return ans;
490 :     }
491 : bates 2114 return R_NilValue; /* -Wall */
492 : maechler 2112 }

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