Matrix  $Rev: 3071 $ at $LastChangedDate: 2015-03-26 15:35:47 +0100 (Thu, 26 Mar 2015) $
Mutils.c
Go to the documentation of this file.
1 #include <limits.h>
2 
3 #include "Mutils.h"
4 #include <R_ext/Lapack.h>
5 
6 // La_norm_type() & La_rcond_type() have been in R_ext/Lapack.h
7 // but have still not been available to package writers ...
8 char La_norm_type(const char *typstr)
9 {
10  char typup;
11 
12  if (strlen(typstr) != 1)
13  error(
14  _("argument type[1]='%s' must be a one-letter character string"),
15  typstr);
16  typup = toupper(*typstr);
17  if (typup == '1')
18  typup = 'O'; /* alias */
19  else if (typup == 'E')
20  typup = 'F';
21  else if (typup != 'M' && typup != 'O' && typup != 'I' && typup != 'F')
22  error(_("argument type[1]='%s' must be one of 'M','1','O','I','F' or 'E'"),
23  typstr);
24  return typup;
25 }
26 
27 char La_rcond_type(const char *typstr)
28 {
29  char typup;
30 
31  if (strlen(typstr) != 1)
32  error(
33  _("argument type[1]='%s' must be a one-letter character string"),
34  typstr);
35  typup = toupper(*typstr);
36  if (typup == '1')
37  typup = 'O'; /* alias */
38  else if (typup != 'O' && typup != 'I')
39  error(_("argument type[1]='%s' must be one of '1','O', or 'I'"),
40  typstr);
41  return typup;
42 }
43 
44 double get_double_by_name(SEXP obj, char *nm)
45 {
46  SEXP nms = getAttrib(obj, R_NamesSymbol);
47  int i, len = length(obj);
48 
49  if ((!isReal(obj)) || (length(obj) > 0 && nms == R_NilValue))
50  error(_("object must be a named, numeric vector"));
51  for (i = 0; i < len; i++) {
52  if (!strcmp(nm, CHAR(STRING_ELT(nms, i)))) {
53  return REAL(obj)[i];
54  }
55  }
56  return R_NaReal;
57 }
58 
59 SEXP
60 set_double_by_name(SEXP obj, double val, char *nm)
61 {
62  SEXP nms = getAttrib(obj, R_NamesSymbol);
63  int i, len = length(obj);
64 
65  if ((!isReal(obj)) || (length(obj) > 0 && nms == R_NilValue))
66  error(_("object must be a named, numeric vector"));
67  for (i = 0; i < len; i++) {
68  if (!strcmp(nm, CHAR(STRING_ELT(nms, i)))) {
69  REAL(obj)[i] = val;
70  return obj;
71  }
72  }
73  {
74  SEXP nx = PROTECT(allocVector(REALSXP, len + 1)),
75  nnms = allocVector(STRSXP, len + 1);
76 
77  setAttrib(nx, R_NamesSymbol, nnms);
78  for (i = 0; i < len; i++) {
79  REAL(nx)[i] = REAL(obj)[i];
80  SET_STRING_ELT(nnms, i, duplicate(STRING_ELT(nms, i)));
81  }
82  REAL(nx)[len] = val;
83  SET_STRING_ELT(nnms, len, mkChar(nm));
84  UNPROTECT(1);
85  return nx;
86  }
87 }
88 
89 SEXP as_det_obj(double val, int log, int sign)
90 {
91  SEXP det = PROTECT(allocVector(VECSXP, 2)),
92  nms = PROTECT(allocVector(STRSXP, 2)),
93  vv = PROTECT(ScalarReal(val));
94 
95  setAttrib(det, R_NamesSymbol, nms);
96  SET_STRING_ELT(nms, 0, mkChar("modulus"));
97  SET_STRING_ELT(nms, 1, mkChar("sign"));
98  setAttrib(vv, install("logarithm"), ScalarLogical(log));
99  SET_VECTOR_ELT(det, 0, vv);
100  SET_VECTOR_ELT(det, 1, ScalarInteger(sign));
101  setAttrib(det, R_ClassSymbol, mkString("det"));
102  UNPROTECT(3);
103  return det;
104 }
105 
106 SEXP get_factors(SEXP obj, char *nm)
107 {
108  SEXP fac = GET_SLOT(obj, Matrix_factorSym),
109  nms = getAttrib(fac, R_NamesSymbol);
110  int i, len = length(fac);
111 
112  if ((!isNewList(fac)) || (length(fac) > 0 && nms == R_NilValue))
113  error(_("'factors' slot must be a named list"));
114  for (i = 0; i < len; i++) {
115  if (!strcmp(nm, CHAR(STRING_ELT(nms, i)))) {
116  return VECTOR_ELT(fac, i);
117  }
118  }
119  return R_NilValue;
120 }
121 
130 SEXP set_factors(SEXP obj, SEXP val, char *nm)
131 {
132  SEXP fac = GET_SLOT(obj, Matrix_factorSym),
133  nms = getAttrib(fac, R_NamesSymbol);
134  int i, len = length(fac);
135 
136  if ((!isNewList(fac)) || (length(fac) > 0 && nms == R_NilValue))
137  error(_("'factors' slot must be a named list"));
138  PROTECT(val); /* set_factors(..) may be called as "finalizer" after UNPROTECT()*/
139  // if there's already a 'nm' factor, we replace it and return:
140  for (i = 0; i < len; i++) {
141  if (!strcmp(nm, CHAR(STRING_ELT(nms, i)))) {
142  SET_VECTOR_ELT(fac, i, duplicate(val));
143  UNPROTECT(1);
144  return val;
145  }
146  }
147  // Otherwise: create a new 'nm' entry in the 'factors' list:
148  // create a list of length (len + 1),
149  SEXP nfac = PROTECT(allocVector(VECSXP, len + 1)),
150  nnms = PROTECT(allocVector(STRSXP, len + 1));
151  setAttrib(nfac, R_NamesSymbol, nnms);
152  // copy all the previous entries,
153  for (i = 0; i < len; i++) {
154  SET_VECTOR_ELT(nfac, i, VECTOR_ELT(fac, i));
155  SET_STRING_ELT(nnms, i, duplicate(STRING_ELT(nms, i)));
156  }
157  // and add the new entry at the end.
158  SET_VECTOR_ELT(nfac, len, duplicate(val));
159  SET_STRING_ELT(nnms, len, mkChar(nm));
160  SET_SLOT(obj, Matrix_factorSym, nfac);
161  UNPROTECT(3);
162  return VECTOR_ELT(nfac, len);
163 }
164 
165 // R interface [for updating the '@ factors' slot of a function *argument* [CARE!]
166 SEXP R_set_factors(SEXP obj, SEXP val, SEXP name, SEXP warn)
167 {
168  Rboolean do_warn = asLogical(warn);
169  if(R_has_slot(obj, Matrix_factorSym))
170  return set_factors(obj, val, (char *)CHAR(asChar(name)));
171  else {
172  if(do_warn) warning(_("Matrix object has no 'factors' slot"));
173  return val;
174  }
175 }
176 
177 #if 0 /* unused */
178 /* useful for all the ..CMatrix classes (and ..R by [0] <-> [1]); but unused */
179 SEXP CMatrix_set_Dim(SEXP x, int nrow)
180 {
181  int *dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
182 
183  dims[0] = nrow;
184  dims[1] = length(GET_SLOT(x, Matrix_pSym)) - 1;
185  return x;
186 }
187 #endif /* unused */
188 
189 /* Fill in the "trivial remainder" in n*m array ;
190  * typically the 'x' slot of a "dtrMatrix", such that
191  * it should be usable for double/logical/int/complex : */
192 #define MAKE_TRIANGULAR_BODY(_TO_, _FROM_, _ZERO_, _ONE_) \
193 { \
194  int i, j, *dims = INTEGER(GET_SLOT(_FROM_, Matrix_DimSym)); \
195  int n = dims[0], m = dims[1]; \
196  \
197  if (*uplo_P(_FROM_) == 'U') { \
198  for (j = 0; j < n; j++) \
199  for (i = j+1; i < m; i++) \
200  _TO_[i + j*m] = _ZERO_; \
201  } else { \
202  for (j = 1; j < n; j++) \
203  for (i = 0; i < j && i < m; i++) \
204  _TO_[i + j*m] = _ZERO_; \
205  } \
206  if (*diag_P(_FROM_) == 'U') { \
207  j = (n < m) ? n : m; \
208  for (i = 0; i < j; i++) \
209  _TO_[i * (m + 1)] = _ONE_; \
210  } \
211 }
212 
213 void make_d_matrix_triangular(double *to, SEXP from)
214  MAKE_TRIANGULAR_BODY(to, from, 0., 1.)
215 void make_i_matrix_triangular(int *to, SEXP from)
216  MAKE_TRIANGULAR_BODY(to, from, 0, 1)
217 
218 
219 /* Should work for double/logical/int/complex : */
220 #define MAKE_SYMMETRIC_BODY(_TO_, _FROM_) \
221 { \
222  int i, j, n = INTEGER(GET_SLOT(_FROM_, Matrix_DimSym))[0]; \
223  \
224  if (*uplo_P(_FROM_) == 'U') { \
225  for (j = 0; j < n; j++) \
226  for (i = j+1; i < n; i++) \
227  _TO_[i + j*n] = _TO_[j + i*n]; \
228  } else { \
229  for (j = 1; j < n; j++) \
230  for (i = 0; i < j && i < n; i++) \
231  _TO_[i + j*n] = _TO_[j + i*n]; \
232  } \
233 }
234 
235 void make_d_matrix_symmetric(double *to, SEXP from)
236  MAKE_SYMMETRIC_BODY(to, from)
237 
238 void make_i_matrix_symmetric(int *to, SEXP from)
239  MAKE_SYMMETRIC_BODY(to, from)
240 
241 
242 #define Matrix_Error_Bufsiz 4096
243 
254 SEXP check_scalar_string(SEXP sP, char *vals, char *nm)
255 {
256  SEXP val = ScalarLogical(1);
257  char *buf;
258  /* only allocate when needed: in good case, none is needed */
259 #define SPRINTF buf = Alloca(Matrix_Error_Bufsiz, char); R_CheckStack(); sprintf
260 
261  if (length(sP) != 1) {
262  SPRINTF(buf, _("'%s' slot must have length 1"), nm);
263  } else {
264  const char *str = CHAR(STRING_ELT(sP, 0));
265  if (strlen(str) != 1) {
266  SPRINTF(buf, _("'%s' must have string length 1"), nm);
267  } else {
268  int i, len;
269  for (i = 0, len = strlen(vals); i < len; i++) {
270  if (str[0] == vals[i])
271  return R_NilValue;
272  }
273  SPRINTF(buf, _("'%s' must be in '%s'"), nm, vals);
274  }
275  }
276  /* 'error' returns : */
277  val = mkString(buf);
278  return val;
279 #undef SPRINTF
280 }
281 
282 /* FIXME? Something like this should be part of the R API ?
283  * But then, R has the more general R_compute_identical()
284  * in src/main/identical.c: Rboolean R_compute_identical(SEXP x, SEXP y);
285 */
286 Rboolean equal_string_vectors(SEXP s1, SEXP s2)
287 {
288  Rboolean n1 = isNull(s1), n2 = isNull(s2);
289  if (n1 || n2) // if one is NULL : "equal" if both are
290  return (n1 == n2) ? TRUE : FALSE;
291  else if (TYPEOF(s1) != STRSXP || TYPEOF(s2) != STRSXP) {
292  error(_("'s1' and 's2' must be \"character\" vectors"));
293  return FALSE; /* -Wall */
294  } else {
295  int n = LENGTH(s1), i;
296  if (n != LENGTH(s2))
297  return FALSE;
298  for(i = 0; i < n; i++) {
299  /* note that compute_identical() code for STRSXP
300  is careful about NA's which we don't need */
301  if(strcmp(CHAR(STRING_ELT(s1, i)),
302  CHAR(STRING_ELT(s2, i))))
303  return FALSE;
304  }
305  return TRUE; /* they *are* equal */
306  }
307 }
308 
309 
311 {
312  int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym));
313  if ((((double) dims[0]) * dims[1]) != XLENGTH(GET_SLOT(obj, Matrix_xSym)))
314  return mkString(_("length of x slot != prod(Dim)"));
315  return ScalarLogical(1);
316 }
317 
318 SEXP dim_validate(SEXP Dim, const char* name) {
319  if (length(Dim) != 2)
320  return mkString(_("Dim slot must have length 2"));
321  /* if (TYPEOF(Dim) != INTSXP && TYPEOF(Dim) != REALSXP) */
322  /* return mkString(_("Dim slot is not numeric")); */
323  if (TYPEOF(Dim) != INTSXP)
324  // TODO?: coerce REALSXP to INTSXP in the "double" case ???
325  return mkString(_("Dim slot is not integer"));
326  int
327  m = INTEGER(Dim)[0],
328  n = INTEGER(Dim)[1];
329  if (m < 0 || n < 0)
330  return mkString(dngettext(name,
331  "Negative value in Dim",
332  "Negative values in Dim",
333  (m*n > 0) ? 2 : 1));
334  return ScalarLogical(1);
335 }
336 // to be called from R :
337 SEXP Dim_validate(SEXP obj, SEXP name) {
338  return dim_validate(GET_SLOT(obj, Matrix_DimSym),
339  CHAR(STRING_ELT(name, 0)));
340 }
341 
342 // obj must have @Dim and @Dimnames. Assume 'Dim' is already checked.
343 SEXP dimNames_validate(SEXP obj)
344 {
345  int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym));
346  SEXP dmNms = GET_SLOT(obj, Matrix_DimNamesSym);
347  if(!isNewList(dmNms))
348  return mkString(_("Dimnames slot is not a list"));
349  if(length(dmNms) != 2)
350  return mkString(_("Dimnames slot is a list, but not of length 2"));
351  char buf[99];
352  for(int j=0; j < 2; j++) { // x@Dimnames[j] must be NULL or character(length(x@Dim[j]))
353  if(!isNull(VECTOR_ELT(dmNms, j))) {
354  if(TYPEOF(VECTOR_ELT(dmNms, j)) != STRSXP) {
355  sprintf(buf, _("Dimnames[%d] is not a character vector"), j+1);
356  return mkString(buf);
357  }
358  if(LENGTH(VECTOR_ELT(dmNms, j)) != 0 && // character(0) allowed here
359  LENGTH(VECTOR_ELT(dmNms, j)) != dims[j]) {
360  sprintf(buf, _("length(Dimnames[%d]) differs from Dim[%d] which is %d"),
361  j+1, j+1, dims[j]);
362  return mkString(buf);
363  }
364  }
365  }
366  return ScalarLogical(1);
367 }
368 
369 
370 
371 #define PACKED_TO_FULL(TYPE) \
372 TYPE *packed_to_full_ ## TYPE(TYPE *dest, const TYPE *src, \
373  int n, enum CBLAS_UPLO uplo) \
374 { \
375  int i, j, pos = 0; \
376  \
377  AZERO(dest, n*n); \
378  for (j = 0; j < n; j++) { \
379  switch(uplo) { \
380  case UPP: \
381  for (i = 0; i <= j; i++) dest[i + j * n] = src[pos++]; \
382  break; \
383  case LOW: \
384  for (i = j; i < n; i++) dest[i + j * n] = src[pos++]; \
385  break; \
386  default: \
387  error(_("'uplo' must be UPP or LOW")); \
388  } \
389  } \
390  return dest; \
391 }
392 
393 PACKED_TO_FULL(double)
394 PACKED_TO_FULL(int)
395 
396 #define FULL_TO_PACKED(TYPE) \
397 TYPE *full_to_packed_ ## TYPE(TYPE *dest, const TYPE *src, int n, \
398  enum CBLAS_UPLO uplo, enum CBLAS_DIAG diag) \
399 { \
400  int i, j, pos = 0; \
401  \
402  for (j = 0; j < n; j++) { \
403  switch(uplo) { \
404  case UPP: \
405  for (i = 0; i <= j; i++) \
406  dest[pos++] = (i == j && diag== UNT) ? 1 : src[i + j*n]; \
407  break; \
408  case LOW: \
409  for (i = j; i < n; i++) \
410  dest[pos++] = (i == j && diag== UNT) ? 1 : src[i + j*n]; \
411  break; \
412  default: \
413  error(_("'uplo' must be UPP or LOW")); \
414  } \
415  } \
416  return dest; \
417 }
418 
419 FULL_TO_PACKED(double)
420 FULL_TO_PACKED(int)
421 
422 
423 
433 void d_packed_getDiag(double *dest, SEXP x, int n)
434 {
435  double *xx = REAL(GET_SLOT(x, Matrix_xSym));
436 
437 #define END_packed_getDiag \
438  int j, pos = 0; \
439  \
440  if (*uplo_P(x) == 'U') { \
441  for(pos= 0, j=0; j < n; pos += 1+(++j)) dest[j] = xx[pos]; \
442  } else { \
443  for(pos= 0, j=0; j < n; pos += (n - j), j++) dest[j] = xx[pos]; \
444  } \
445  return
446 
448 }
449 
450 void l_packed_getDiag(int *dest, SEXP x, int n)
451 {
452  int *xx = LOGICAL(GET_SLOT(x, Matrix_xSym));
453 
455 }
456 
457 #undef END_packed_getDiag
458 
459 //----------------------------------------------------------------------
460 
461 SEXP d_packed_setDiag(double *diag, int l_d, SEXP x, int n)
462 {
463 #define SET_packed_setDiag \
464  SEXP ret = PROTECT(duplicate(x)), \
465  r_x = GET_SLOT(ret, Matrix_xSym); \
466  Rboolean d_full = (l_d == n); \
467  if (!d_full && l_d != 1) \
468  error(_("replacement diagonal has wrong length"))
469 
470 #define END_packed_setDiag \
471  int j, pos = 0; \
472  \
473  if (*uplo_P(x) == 'U') { \
474  if(d_full) \
475  for(pos= 0, j=0; j < n; pos += 1+(++j)) xx[pos] = diag[j]; \
476  else /* l_d == 1 */ \
477  for(pos= 0, j=0; j < n; pos += 1+(++j)) xx[pos] = *diag; \
478  } else { \
479  if(d_full) \
480  for(pos= 0, j=0; j < n; pos += (n - j), j++) xx[pos] = diag[j]; \
481  else /* l_d == 1 */ \
482  for(pos= 0, j=0; j < n; pos += (n - j), j++) xx[pos] = *diag; \
483  } \
484  UNPROTECT(1); \
485  return ret
486 
487  SET_packed_setDiag; double *xx = REAL(r_x);
489 }
490 
491 SEXP l_packed_setDiag(int *diag, int l_d, SEXP x, int n)
492 {
493  SET_packed_setDiag; int *xx = LOGICAL(r_x);
495 }
496 
497 #define tr_END_packed_setDiag \
498  if (*diag_P(x) == 'U') { /* uni-triangular */ \
499  /* after setting, typically is not uni-triangular anymore: */ \
500  SET_STRING_ELT(GET_SLOT(ret, Matrix_diagSym), 0, mkChar("N")); \
501  } \
502  END_packed_setDiag
503 
504 
505 SEXP tr_d_packed_setDiag(double *diag, int l_d, SEXP x, int n)
506 {
507  SET_packed_setDiag; double *xx = REAL(r_x);
509 }
510 
511 SEXP tr_l_packed_setDiag(int *diag, int l_d, SEXP x, int n)
512 {
513  SET_packed_setDiag; int *xx = LOGICAL(r_x);
515 }
516 
517 
518 #undef SET_packed_setDiag
519 #undef END_packed_setDiag
520 #undef tr_END_packed_setDiag
521 //----------------------------------------------------------------------
522 
523 SEXP d_packed_addDiag(double *diag, int l_d, SEXP x, int n)
524 {
525  SEXP ret = PROTECT(duplicate(x)),
526  r_x = GET_SLOT(ret, Matrix_xSym);
527  double *xx = REAL(r_x);
528  int j, pos = 0;
529 
530  if (*uplo_P(x) == 'U') {
531  for(pos= 0, j=0; j < n; pos += 1+(++j)) xx[pos] += diag[j];
532  } else {
533  for(pos= 0, j=0; j < n; pos += (n - j), j++) xx[pos] += diag[j];
534  }
535  UNPROTECT(1);
536  return ret;
537 }
538 
539 SEXP tr_d_packed_addDiag(double *diag, int l_d, SEXP x, int n)
540 {
541  SEXP ret = PROTECT(d_packed_addDiag(diag, l_d, x, n));
542  if (*diag_P(x) == 'U') /* uni-triangular */
543  SET_STRING_ELT(GET_SLOT(ret, Matrix_diagSym), 0, mkChar("N"));
544  UNPROTECT(1);
545  return ret;
546 }
547 
548 
549 //----------------------------------------------------------------------
550 
551 void tr_d_packed_getDiag(double *dest, SEXP x, int n)
552 {
553  if (*diag_P(x) == 'U') {
554  for (int j = 0; j < n; j++) dest[j] = 1.;
555  } else {
556  d_packed_getDiag(dest, x, n);
557  }
558  return;
559 }
560 
561 void tr_l_packed_getDiag( int *dest, SEXP x, int n)
562 {
563  if (*diag_P(x) == 'U')
564  for (int j = 0; j < n; j++) dest[j] = 1;
565  else
566  l_packed_getDiag(dest, x, n);
567  return;
568 }
569 
570 
572 {
573  int n = length(pP) - 1;
574  int *p = INTEGER(pP);
575  SEXP ans = PROTECT(allocVector(INTSXP, p[n]));
576 
577  expand_cmprPt(n, p, INTEGER(ans));
578  UNPROTECT(1);
579  return ans;
580 }
581 
582 
591 SEXP
592 Matrix_getElement(SEXP list, char *nm) {
593  SEXP names = getAttrib(list, R_NamesSymbol);
594  int i;
595 
596  for (i = 0; i < LENGTH(names); i++)
597  if (!strcmp(CHAR(STRING_ELT(names, i)), nm))
598  return(VECTOR_ELT(list, i));
599  return R_NilValue;
600 }
601 
610 static double *
611 install_diagonal(double *dest, SEXP A)
612 {
613  int nc = INTEGER(GET_SLOT(A, Matrix_DimSym))[0];
614  int i, ncp1 = nc + 1, unit = *diag_P(A) == 'U';
615  double *ax = REAL(GET_SLOT(A, Matrix_xSym));
616 
617  AZERO(dest, nc * nc);
618  for (i = 0; i < nc; i++)
619  dest[i * ncp1] = (unit) ? 1. : ax[i];
620  return dest;
621 }
622 
623 static int *
624 install_diagonal_int(int *dest, SEXP A)
625 {
626  int nc = INTEGER(GET_SLOT(A, Matrix_DimSym))[0];
627  int i, ncp1 = nc + 1, unit = *diag_P(A) == 'U';
628  int *ax = INTEGER(GET_SLOT(A, Matrix_xSym));
629 
630  AZERO(dest, nc * nc);
631  for (i = 0; i < nc; i++)
632  dest[i * ncp1] = (unit) ? 1 : ax[i];
633  return dest;
634 }
635 
636 
649 {
650  /* NOTA BENE: If you enlarge this list, do change '14' and '6' below !
651  * --------- ddiMatrix & ldiMatrix are no longer ddense or ldense on the R level,
652  * --- --- but are still dealt with here: */
653  static const char *valid[] = {
654  "_NOT_A_CLASS_",// *_CLASSES defined in ./Mutils.h :
655  MATRIX_VALID_ddense, /* 14 */
656  MATRIX_VALID_ldense, /* 6 */
657  MATRIX_VALID_ndense, /* 5 */
658  ""};
659  SEXP ans, ad = R_NilValue, an = R_NilValue; /* -Wall */
660  int sz, ctype = R_check_class_etc(A, valid),
661  nprot = 1;
662  enum dense_enum M_type = ddense /* -Wall */;
663 
664  if (ctype > 0) { /* a [nld]denseMatrix or [dl]diMatrix object */
665  ad = GET_SLOT(A, Matrix_DimSym);
666  an = GET_SLOT(A, Matrix_DimNamesSym);
667  M_type = (ctype <= 14) ? ddense :
668  ((ctype <= 14+6) ? ldense : ndense);
669  }
670  else if (ctype < 0) { /* not a (recognized) classed matrix */
671 
672  if (isReal(A))
673  M_type = ddense;
674  else if (isInteger(A)) {
675  A = PROTECT(coerceVector(A, REALSXP));
676  nprot++;
677  M_type = ddense;
678  }
679  else if (isLogical(A))
680  M_type = ldense;
681  else
682  error(_("invalid class '%s' to dup_mMatrix_as_geMatrix"),
683  class_P(A));
684 
685 #define DUP_MMATRIX_NON_CLASS(transpose_if_vec) \
686  if (isMatrix(A)) { /* "matrix" */ \
687  ad = getAttrib(A, R_DimSymbol); \
688  an = getAttrib(A, R_DimNamesSymbol); \
689  } else {/* maybe "numeric" (incl integer,logical) --> (n x 1) */ \
690  int* dd = INTEGER(ad = PROTECT(allocVector(INTSXP, 2))); \
691  nprot++; \
692  if(transpose_if_vec) { \
693  dd[0] = 1; \
694  dd[1] = LENGTH(A); \
695  } else { \
696  dd[0] = LENGTH(A); \
697  dd[1] = 1; \
698  } \
699  SEXP nms = getAttrib(A, R_NamesSymbol); \
700  if(nms != R_NilValue) { \
701  an = PROTECT(allocVector(VECSXP, 2)); \
702  nprot++; \
703  SET_VECTOR_ELT(an, (transpose_if_vec)? 1 : 0, nms); \
704  /* not needed: SET_VECTOR_ELT(an, 1, R_NilValue); */ \
705  } /* else nms = NULL ==> an remains NULL */ \
706  } \
707  ctype = 0
708 
709  DUP_MMATRIX_NON_CLASS(FALSE);
710  }
711 
712  ans = PROTECT(NEW_OBJECT(MAKE_CLASS(M_type == ddense ? "dgeMatrix" :
713  (M_type == ldense ? "lgeMatrix" :
714  "ngeMatrix"))));
715 #define DUP_MMATRIX_SET_1 \
716  SET_SLOT(ans, Matrix_DimSym, duplicate(ad)); \
717  SET_SLOT(ans, Matrix_DimNamesSym, (LENGTH(an) == 2) ? \
718  duplicate(an): allocVector(VECSXP, 2)); \
719  sz = INTEGER(ad)[0] * INTEGER(ad)[1]
720 
722 
723  if(M_type == ddense) { /* ddense -> dge */
724 
725  double *ansx;
726 
727 #define DUP_MMATRIX_ddense_CASES \
728  ansx = REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, sz)); \
729  switch(ctype) { \
730  case 0: /* unclassed real matrix */ \
731  Memcpy(ansx, REAL(A), sz); \
732  break; \
733  case 1: /* dgeMatrix */ \
734  Memcpy(ansx, REAL(GET_SLOT(A, Matrix_xSym)), sz); \
735  break; \
736  case 2: /* dtrMatrix and subclasses */ \
737  case 9: case 10: case 11: /* --- Cholesky, LDL, BunchKaufman */ \
738  Memcpy(ansx, REAL(GET_SLOT(A, Matrix_xSym)), sz); \
739  make_d_matrix_triangular(ansx, A); \
740  break; \
741  case 3: /* dsyMatrix */ \
742  case 4: /* dpoMatrix + subclass */ \
743  case 14: /* --- corMatrix */ \
744  Memcpy(ansx, REAL(GET_SLOT(A, Matrix_xSym)), sz); \
745  make_d_matrix_symmetric(ansx, A); \
746  break; \
747  case 5: /* ddiMatrix */ \
748  install_diagonal(ansx, A); \
749  break; \
750  case 6: /* dtpMatrix + subclasses */ \
751  case 12: case 13: /* --- pCholesky, pBunchKaufman */ \
752  packed_to_full_double(ansx, REAL(GET_SLOT(A, Matrix_xSym)), \
753  INTEGER(ad)[0], \
754  *uplo_P(A) == 'U' ? UPP : LOW); \
755  make_d_matrix_triangular(ansx, A); \
756  break; \
757  case 7: /* dspMatrix */ \
758  case 8: /* dppMatrix */ \
759  packed_to_full_double(ansx, REAL(GET_SLOT(A, Matrix_xSym)), \
760  INTEGER(ad)[0], \
761  *uplo_P(A) == 'U' ? UPP : LOW); \
762  make_d_matrix_symmetric(ansx, A); \
763  break; \
764  } /* switch(ctype) */
765 
767  }
768  else { /* M_type == ldense || M_type = ndense */
769  /* ldense -> lge */
770  /* ndense -> nge */
771  int *ansx = LOGICAL(ALLOC_SLOT(ans, Matrix_xSym, LGLSXP, sz));
772 
773  switch(ctype) {
774  case 0: /* unclassed logical matrix */
775  Memcpy(ansx, LOGICAL(A), sz);
776  break;
777 
778  case 1+14: /* lgeMatrix */
779  case 1+14+6: /* ngeMatrix */
780  Memcpy(ansx, LOGICAL(GET_SLOT(A, Matrix_xSym)), sz);
781  break;
782  case 2+14: /* ltrMatrix */
783  case 2+14+6: /* ntrMatrix */
784  Memcpy(ansx, LOGICAL(GET_SLOT(A, Matrix_xSym)), sz);
785  make_i_matrix_triangular(ansx, A);
786  break;
787  case 3+14: /* lsyMatrix */
788  case 3+14+6: /* nsyMatrix */
789  Memcpy(ansx, LOGICAL(GET_SLOT(A, Matrix_xSym)), sz);
790  make_i_matrix_symmetric(ansx, A);
791  break;
792  case 4+14: /* ldiMatrix */
793  // case 4+14+6: /* ndiMatrix _DOES NOT EXIST_ */
794  install_diagonal_int(ansx, A);
795  break;
796  case 5+14: /* ltpMatrix */
797  case 4+14+6: /* ntpMatrix */
798  packed_to_full_int(ansx, LOGICAL(GET_SLOT(A, Matrix_xSym)),
799  INTEGER(ad)[0],
800  *uplo_P(A) == 'U' ? UPP : LOW);
801  make_i_matrix_triangular(ansx, A);
802  break;
803  case 6+14: /* lspMatrix */
804  case 5+14+6: /* nspMatrix */
805  packed_to_full_int(ansx, LOGICAL(GET_SLOT(A, Matrix_xSym)),
806  INTEGER(ad)[0],
807  *uplo_P(A) == 'U' ? UPP : LOW);
808  make_i_matrix_symmetric(ansx, A);
809  break;
810 
811  default:
812  error(_("unexpected ctype = %d in dup_mMatrix_as_geMatrix"), ctype);
813  } /* switch(ctype) */
814 
815  } /* if(M_type == .) */
816 
817  UNPROTECT(nprot);
818  return ans;
819 }
820 
821 SEXP dup_mMatrix_as_dgeMatrix2(SEXP A, Rboolean tr_if_vec)
822 {
823  SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix"))),
824  ad = R_NilValue , an = R_NilValue; /* -Wall */
825  static const char *valid[] = {"_NOT_A_CLASS_", MATRIX_VALID_ddense, ""};
826  int ctype = R_check_class_etc(A, valid),
827  nprot = 1, sz;
828  double *ansx;
829 
830  if (ctype > 0) { /* a ddenseMatrix object */
831  ad = GET_SLOT(A, Matrix_DimSym);
832  an = GET_SLOT(A, Matrix_DimNamesSym);
833  }
834  else if (ctype < 0) { /* not a (recognized) classed matrix */
835  if (!isReal(A)) {
836  if (isInteger(A) || isLogical(A)) {
837  A = PROTECT(coerceVector(A, REALSXP));
838  nprot++;
839  } else
840  error(_("invalid class '%s' to dup_mMatrix_as_dgeMatrix"),
841  class_P(A));
842  }
843  DUP_MMATRIX_NON_CLASS(tr_if_vec);
844  }
845 
848  UNPROTECT(nprot);
849  return ans;
850 }
851 
853  return dup_mMatrix_as_dgeMatrix2(A, FALSE);
854 }
855 
856 SEXP new_dgeMatrix(int nrow, int ncol)
857 {
858  SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix"))),
859  ad = PROTECT(allocVector(INTSXP, 2));
860 
861  INTEGER(ad)[0] = nrow;
862  INTEGER(ad)[1] = ncol;
863  SET_SLOT(ans, Matrix_DimSym, ad);
864  SET_SLOT(ans, Matrix_DimNamesSym, allocVector(VECSXP, 2));
865  ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nrow * ncol);
866 
867  UNPROTECT(2);
868  return ans;
869 }
870 
880 SEXP m_encodeInd(SEXP ij, SEXP di, SEXP orig_1, SEXP chk_bnds)
881 {
882  SEXP ans;
883  int *ij_di = NULL, n, nprot=1;
884  Rboolean check_bounds = asLogical(chk_bnds), one_ind = asLogical(orig_1);
885 
886  if(TYPEOF(di) != INTSXP) {di = PROTECT(coerceVector(di, INTSXP)); nprot++; }
887  if(TYPEOF(ij) != INTSXP) {ij = PROTECT(coerceVector(ij, INTSXP)); nprot++; }
888  if(!isMatrix(ij) ||
889  (ij_di = INTEGER(getAttrib(ij, R_DimSymbol)))[1] != 2)
890  error(_("Argument ij must be 2-column integer matrix"));
891  n = ij_di[0];
892  int *Di = INTEGER(di), *IJ = INTEGER(ij),
893  *j_ = IJ+n;/* pointer offset! */
894 
895  if((Di[0] * (double) Di[1]) >= 1 + (double)INT_MAX) { /* need double */
896  ans = PROTECT(allocVector(REALSXP, n));
897  double *ii = REAL(ans), nr = (double) Di[0];
898 #define do_ii_FILL(_i_, _j_) \
899  int i; \
900  if(check_bounds) { \
901  for(i=0; i < n; i++) { \
902  if(_i_[i] == NA_INTEGER || _j_[i] == NA_INTEGER) \
903  ii[i] = NA_INTEGER; \
904  else { \
905  register int i_i, j_i; \
906  if(one_ind) { i_i = _i_[i]-1; j_i = _j_[i]-1; } \
907  else { i_i = _i_[i] ; j_i = _j_[i] ; } \
908  if(i_i < 0 || i_i >= Di[0]) \
909  error(_("subscript 'i' out of bounds in M[ij]")); \
910  if(j_i < 0 || j_i >= Di[1]) \
911  error(_("subscript 'j' out of bounds in M[ij]")); \
912  ii[i] = i_i + j_i * nr; \
913  } \
914  } \
915  } else { \
916  for(i=0; i < n; i++) \
917  ii[i] = (_i_[i] == NA_INTEGER || _j_[i] == NA_INTEGER) \
918  ? NA_INTEGER \
919  : (one_ind ? ((_i_[i]-1) + (_j_[i]-1)*nr) \
920  : _i_[i] + _j_[i] *nr); \
921  }
922 
923  do_ii_FILL(IJ, j_);
924  } else {
925  ans = PROTECT(allocVector(INTSXP, n));
926  int *ii = INTEGER(ans), nr = Di[0];
927 
928  do_ii_FILL(IJ, j_);
929  }
930  UNPROTECT(nprot);
931  return ans;
932 }
933 
945 SEXP m_encodeInd2(SEXP i, SEXP j, SEXP di, SEXP orig_1, SEXP chk_bnds)
946 {
947  SEXP ans;
948  int n = LENGTH(i), nprot = 1;
949  Rboolean check_bounds = asLogical(chk_bnds), one_ind = asLogical(orig_1);
950 
951  if(TYPEOF(di)!= INTSXP) {di = PROTECT(coerceVector(di,INTSXP)); nprot++; }
952  if(TYPEOF(i) != INTSXP) { i = PROTECT(coerceVector(i, INTSXP)); nprot++; }
953  if(TYPEOF(j) != INTSXP) { j = PROTECT(coerceVector(j, INTSXP)); nprot++; }
954  if(LENGTH(j) != n)
955  error(_("i and j must be integer vectors of the same length"));
956  int *Di = INTEGER(di), *i_ = INTEGER(i), *j_ = INTEGER(j);
957 
958  if((Di[0] * (double) Di[1]) >= 1 + (double)INT_MAX) { /* need double */
959  ans = PROTECT(allocVector(REALSXP, n));
960  double *ii = REAL(ans), nr = (double) Di[0];
961 
962  do_ii_FILL(i_, j_);
963  } else {
964  ans = PROTECT(allocVector(INTSXP, n));
965  int *ii = INTEGER(ans), nr = Di[0];
966 
967  do_ii_FILL(i_, j_);
968  }
969  UNPROTECT(nprot);
970  return ans;
971 }
972 #undef do_ii_FILL
973 
974 // Almost "Cut n Paste" from ...R../src/main/array.c do_matrix() :
975 // used in ../R/Matrix.R as
976 //
977 // .External(Mmatrix,
978 // data, nrow, ncol, byrow, dimnames,
979 // missing(nrow), missing(ncol))
980 SEXP Mmatrix(SEXP args)
981 {
982  SEXP vals, ans, snr, snc, dimnames;
983  int nr = 1, nc = 1, byrow, miss_nr, miss_nc;
984  R_xlen_t lendat;
985 
986  args = CDR(args); /* skip 'name' */
987  vals = CAR(args); args = CDR(args);
988  /* Supposedly as.vector() gave a vector type, but we check */
989  switch(TYPEOF(vals)) {
990  case LGLSXP:
991  case INTSXP:
992  case REALSXP:
993  case CPLXSXP:
994  case STRSXP:
995  case RAWSXP:
996  case EXPRSXP:
997  case VECSXP:
998  break;
999  default:
1000  error(_("'data' must be of a vector type"));
1001  }
1002  lendat = XLENGTH(vals);
1003  snr = CAR(args); args = CDR(args);
1004  snc = CAR(args); args = CDR(args);
1005  byrow = asLogical(CAR(args)); args = CDR(args);
1006  if (byrow == NA_INTEGER)
1007  error(_("invalid '%s' argument"), "byrow");
1008  dimnames = CAR(args);
1009  args = CDR(args);
1010  miss_nr = asLogical(CAR(args)); args = CDR(args);
1011  miss_nc = asLogical(CAR(args));
1012 
1013  if (!miss_nr) {
1014  if (!isNumeric(snr)) error(_("non-numeric matrix extent"));
1015  nr = asInteger(snr);
1016  if (nr == NA_INTEGER)
1017  error(_("invalid 'nrow' value (too large or NA)"));
1018  if (nr < 0)
1019  error(_("invalid 'nrow' value (< 0)"));
1020  }
1021  if (!miss_nc) {
1022  if (!isNumeric(snc)) error(_("non-numeric matrix extent"));
1023  nc = asInteger(snc);
1024  if (nc == NA_INTEGER)
1025  error(_("invalid 'ncol' value (too large or NA)"));
1026  if (nc < 0)
1027  error(_("invalid 'ncol' value (< 0)"));
1028  }
1029  if (miss_nr && miss_nc) {
1030  if (lendat > INT_MAX) error("data is too long");
1031  nr = (int) lendat;
1032  } else if (miss_nr) {
1033  if (lendat > (double) nc * INT_MAX) error("data is too long");
1034  nr = (int) ceil((double) lendat / (double) nc);
1035  } else if (miss_nc) {
1036  if (lendat > (double) nr * INT_MAX) error("data is too long");
1037  nc = (int) ceil((double) lendat / (double) nr);
1038  }
1039 
1040  if(lendat > 0) {
1041  R_xlen_t nrc = (R_xlen_t) nr * nc;
1042  if (lendat > 1 && nrc % lendat != 0) {
1043  if (((lendat > nr) && (lendat / nr) * nr != lendat) ||
1044  ((lendat < nr) && (nr / lendat) * lendat != nr))
1045  warning(_("data length [%d] is not a sub-multiple or multiple of the number of rows [%d]"), lendat, nr);
1046  else if (((lendat > nc) && (lendat / nc) * nc != lendat) ||
1047  ((lendat < nc) && (nc / lendat) * lendat != nc))
1048  warning(_("data length [%d] is not a sub-multiple or multiple of the number of columns [%d]"), lendat, nc);
1049  }
1050  else if ((lendat > 1) && (nrc == 0)){
1051  warning(_("data length exceeds size of matrix"));
1052  }
1053  }
1054 
1055 #ifndef LONG_VECTOR_SUPPORT
1056  if ((double)nr * (double)nc > INT_MAX)
1057  error(_("too many elements specified"));
1058 #endif
1059 
1060  PROTECT(ans = allocMatrix(TYPEOF(vals), nr, nc));
1061  if(lendat) {
1062  if (isVector(vals))
1063  copyMatrix(ans, vals, byrow);
1064  else
1065  copyListMatrix(ans, vals, byrow);
1066  } else if (isVector(vals)) { /* fill with NAs */
1067  R_xlen_t N = (R_xlen_t) nr * nc, i;
1068  switch(TYPEOF(vals)) {
1069  case STRSXP:
1070  for (i = 0; i < N; i++)
1071  SET_STRING_ELT(ans, i, NA_STRING);
1072  break;
1073  case LGLSXP:
1074  for (i = 0; i < N; i++)
1075  LOGICAL(ans)[i] = NA_LOGICAL;
1076  break;
1077  case INTSXP:
1078  for (i = 0; i < N; i++)
1079  INTEGER(ans)[i] = NA_INTEGER;
1080  break;
1081  case REALSXP:
1082  for (i = 0; i < N; i++)
1083  REAL(ans)[i] = NA_REAL;
1084  break;
1085  case CPLXSXP:
1086  {
1087  Rcomplex na_cmplx;
1088  na_cmplx.r = NA_REAL;
1089  na_cmplx.i = 0;
1090  for (i = 0; i < N; i++)
1091  COMPLEX(ans)[i] = na_cmplx;
1092  }
1093  break;
1094  case RAWSXP:
1095  // FIXME: N may overflow size_t !!
1096  memset(RAW(ans), 0, N);
1097  break;
1098  default:
1099  /* don't fill with anything */
1100  ;
1101  }
1102  }
1103  if(!isNull(dimnames)&& length(dimnames) > 0)
1104  ans = dimnamesgets(ans, dimnames);
1105  UNPROTECT(1);
1106  return ans;
1107 }
1108 SEXP R_rbind2_vector(SEXP a, SEXP b) {
1121  int *d_a = INTEGER(GET_SLOT(a, Matrix_DimSym)),
1122  *d_b = INTEGER(GET_SLOT(b, Matrix_DimSym)),
1123  n1 = d_a[0], m = d_a[1],
1124  n2 = d_b[0],
1125  nprot = 1;
1126  SEXP ans,
1127  a_x = GET_SLOT(a, Matrix_xSym),
1128  b_x = GET_SLOT(b, Matrix_xSym);
1129  if(d_b[1] != m)
1130  error(_("the number of columns differ in R_rbind2_vector: %d != %d"),
1131  m, d_b[1]);
1132  // Care: can have "ddenseMatrix" "ldenseMatrix" or "ndenseMatrix"
1133  if(TYPEOF(a_x) != TYPEOF(b_x)) { // choose the "common type"
1134  // Now know: either LGLSXP or REALSXP. FIXME for iMatrix, zMatrix,..
1135  if(TYPEOF(a_x) != REALSXP)
1136  a_x = PROTECT(duplicate(coerceVector(a_x, REALSXP)));
1137  else if(TYPEOF(b_x) != REALSXP)
1138  b_x = PROTECT(duplicate(coerceVector(b_x, REALSXP)));
1139  nprot++;
1140  }
1141 
1142  ans = PROTECT(allocVector(TYPEOF(a_x), m * (n1 + n2)));
1143  int i, ii = 0;
1144  switch(TYPEOF(a_x)) {
1145  case LGLSXP: {
1146  int
1147  *r = LOGICAL(ans),
1148  *ax= LOGICAL(a_x),
1149  *bx= LOGICAL(b_x);
1150 
1151 #define COPY_a_AND_b_j \
1152 /* FIXME faster: use Memcpy() : */ \
1153  for(int j=0; j < m; j++) { \
1154  for(i=0; i < n1; i++) r[ii++] = ax[j*n1 + i]; \
1155  for(i=0; i < n2; i++) r[ii++] = bx[j*n2 + i]; \
1156  }
1157 
1159  }
1160  case REALSXP: {
1161  double
1162  *r = REAL(ans),
1163  *ax= REAL(a_x),
1164  *bx= REAL(b_x);
1165 
1167  }
1168  } // switch
1169 
1170  UNPROTECT(nprot);
1171  return ans;
1172 }
1173 
1174 #define TRUE_ ScalarLogical(1)
1175 #define FALSE_ ScalarLogical(0)
1176 
1177 // Fast implementation of [ originally in ../R/Auxiliaries.R ]
1178 // all0 <- function(x) !any(is.na(x)) && all(!x) ## ~= allFalse
1179 // allFalse <- function(x) !any(x) && !any(is.na(x)) ## ~= all0
1180 SEXP R_all0(SEXP x) {
1181  if (!isVectorAtomic(x)) {
1182  if(length(x) == 0) return TRUE_;
1183  // Typically S4. TODO: Call the R code above, instead!
1184  error(_("Argument must be numeric-like atomic vector"));
1185  }
1186  R_xlen_t i, n = XLENGTH(x);
1187  if(n == 0) return TRUE_;
1188 
1189  switch(TYPEOF(x)) {
1190  case LGLSXP: {
1191  int *xx = LOGICAL(x);
1192  for(i=0; i < n; i++)
1193  if(xx[i] == NA_LOGICAL || xx[i] != 0) return FALSE_;
1194  return TRUE_;
1195  }
1196  case INTSXP: {
1197  int *xx = INTEGER(x);
1198  for(i=0; i < n; i++)
1199  if(xx[i] == NA_INTEGER || xx[i] != 0) return FALSE_;
1200  return TRUE_;
1201  }
1202  case REALSXP: {
1203  double *xx = REAL(x);
1204  for(i=0; i < n; i++)
1205  if(ISNAN(xx[i]) || xx[i] != 0.) return FALSE_;
1206  return TRUE_;
1207  }
1208  case RAWSXP: {
1209  unsigned char *xx = RAW(x);
1210  for(i=0; i < n; i++)
1211  if(xx[i] != 0) return FALSE_;
1212  return TRUE_;
1213  }
1214  }
1215  error(_("Argument must be numeric-like atomic vector"));
1216  return R_NilValue; // -Wall
1217 }
1218 
1219 // Fast implementation of [ originally in ../R/Auxiliaries.R ]
1220 // any0 <- function(x) isTRUE(any(x == 0)) ## ~= anyFalse
1221 // anyFalse <- function(x) isTRUE(any(!x)) ## ~= any0
1222 SEXP R_any0(SEXP x) {
1223  if (!isVectorAtomic(x)) {
1224  if(length(x) == 0) return FALSE_;
1225  // Typically S4. TODO: Call the R code above, instead!
1226  error(_("Argument must be numeric-like atomic vector"));
1227  }
1228  R_xlen_t i, n = XLENGTH(x);
1229  if(n == 0) return FALSE_;
1230 
1231  switch(TYPEOF(x)) {
1232  case LGLSXP: {
1233  int *xx = LOGICAL(x);
1234  for(i=0; i < n; i++) if(xx[i] == 0) return TRUE_;
1235  return FALSE_;
1236  }
1237  case INTSXP: {
1238  int *xx = INTEGER(x);
1239  for(i=0; i < n; i++) if(xx[i] == 0) return TRUE_;
1240  return FALSE_;
1241  }
1242  case REALSXP: {
1243  double *xx = REAL(x);
1244  for(i=0; i < n; i++) if(xx[i] == 0.) return TRUE_;
1245  return FALSE_;
1246  }
1247  case RAWSXP: {
1248  unsigned char *xx = RAW(x);
1249  for(i=0; i < n; i++) if(xx[i] == 0) return TRUE_;
1250  return FALSE_;
1251  }
1252  }
1253  error(_("Argument must be numeric-like atomic vector"));
1254  return R_NilValue; // -Wall
1255 }
1256 
1257 #undef TRUE_
1258 #undef FALSE_
1259 
1260 /* FIXME: Compare and synchronize with MK_SYMMETRIC_DIMNAMES.. in ./dense.c
1261  * ----- which *also* considers names(dimnames(.)) !!
1262 */
1263 
1269 SEXP symmetric_DimNames(SEXP dn) {
1270  Rboolean do_nm = FALSE;
1271 #define NON_TRIVIAL_DN \
1272  !isNull(VECTOR_ELT(dn, 0)) || \
1273  !isNull(VECTOR_ELT(dn, 1)) || \
1274  (do_nm = !isNull(getAttrib(dn, R_NamesSymbol)))
1275 
1276 #define SYMM_DIMNAMES \
1277  /* Fixup dimnames to be symmetric <==> \
1278  symmetricDimnames() in ../R/symmetricMatrix.R */ \
1279  PROTECT(dn = duplicate(dn)); \
1280  if (isNull(VECTOR_ELT(dn,0))) \
1281  SET_VECTOR_ELT(dn, 0, VECTOR_ELT(dn, 1)); \
1282  if (isNull(VECTOR_ELT(dn,1))) \
1283  SET_VECTOR_ELT(dn, 1, VECTOR_ELT(dn, 0)); \
1284  if(do_nm) { /* names(dimnames(.)): */ \
1285  SEXP nms_dn = getAttrib(dn, R_NamesSymbol); \
1286  if(!R_compute_identical(STRING_ELT(nms_dn, 0), \
1287  STRING_ELT(nms_dn, 1), 16)) { \
1288  PROTECT(nms_dn); \
1289  int J = LENGTH(STRING_ELT(nms_dn, 0)) == 0; /* 0/1 */ \
1290  SET_STRING_ELT(nms_dn, !J, STRING_ELT(nms_dn, J)); \
1291  setAttrib(dn, R_NamesSymbol, nms_dn); \
1292  UNPROTECT(1); \
1293  } \
1294  } \
1295  UNPROTECT(1)
1296 
1297  // Be fast (do nothing!) for the case where dimnames = list(NULL,NULL) :
1298  if (NON_TRIVIAL_DN) {
1299  SYMM_DIMNAMES;
1300  }
1301  return dn;
1302 }
1303 
1315 SEXP R_symmetric_Dimnames(SEXP x) {
1316  return symmetric_DimNames(GET_SLOT(x, Matrix_DimNamesSym));
1317 }
1318 
1319 
1328 void SET_DimNames_symm(SEXP dest, SEXP src) {
1329  SEXP dn = GET_SLOT(src, Matrix_DimNamesSym);
1330  Rboolean do_nm = FALSE;
1331  // Be fast (do nothing!) for the case where dimnames = list(NULL,NULL) :
1332  if (NON_TRIVIAL_DN) {
1333  SYMM_DIMNAMES;
1334  SET_SLOT(dest, Matrix_DimNamesSym, dn);
1335  }
1336  return;
1337 }
1338 
SEXP Matrix_DimSym
Definition: Syms.h:2
#define class_P(_x_)
Definition: Mutils.h:178
static double * install_diagonal(double *dest, SEXP A)
Zero a square matrix of size nc then copy a vector to the diagonal.
Definition: Mutils.c:611
SEXP d_packed_addDiag(double *diag, int l_d, SEXP x, int n)
Definition: Mutils.c:523
SEXP as_det_obj(double val, int log, int sign)
Definition: Mutils.c:89
void tr_l_packed_getDiag(int *dest, SEXP x, int n)
Definition: Mutils.c:561
SEXP Matrix_expand_pointers(SEXP pP)
Definition: Mutils.c:571
SEXP Matrix_xSym
Definition: Syms.h:2
#define FALSE_
Definition: Mutils.c:1175
#define tr_END_packed_setDiag
Definition: Mutils.c:497
#define END_packed_getDiag
#define MATRIX_VALID_ldense
Definition: Mutils.h:369
SEXP m_encodeInd(SEXP ij, SEXP di, SEXP orig_1, SEXP chk_bnds)
Encode Matrix index (i,j) |–> i + j * nrow {i,j : 0-origin}.
Definition: Mutils.c:880
SEXP R_all0(SEXP x)
Definition: Mutils.c:1180
SEXP R_any0(SEXP x)
Definition: Mutils.c:1222
SEXP Matrix_factorSym
Definition: Syms.h:2
#define XLENGTH(x)
Definition: Mutils.h:39
SEXP Mmatrix(SEXP args)
Definition: Mutils.c:980
SEXP Matrix_DimNamesSym
Definition: Syms.h:2
static int * install_diagonal_int(int *dest, SEXP A)
Definition: Mutils.c:624
void make_i_matrix_symmetric(int *to, SEXP from)
SEXP tr_d_packed_addDiag(double *diag, int l_d, SEXP x, int n)
Definition: Mutils.c:539
void l_packed_getDiag(int *dest, SEXP x, int n)
Definition: Mutils.c:450
#define FULL_TO_PACKED(TYPE)
Definition: Mutils.c:396
SEXP dup_mMatrix_as_dgeMatrix(SEXP A)
Definition: Mutils.c:852
SEXP tr_l_packed_setDiag(int *diag, int l_d, SEXP x, int n)
Definition: Mutils.c:511
SEXP dim_validate(SEXP Dim, const char *name)
Definition: Mutils.c:318
dense_enum
Definition: Mutils.h:181
#define END_packed_setDiag
#define SYMM_DIMNAMES
void make_d_matrix_triangular(double *to, SEXP from) void make_i_matrix_triangular(int *to
void make_i_matrix_triangular(int *x, SEXP from)
Rboolean equal_string_vectors(SEXP s1, SEXP s2)
Definition: Mutils.c:286
SEXP check_scalar_string(SEXP sP, char *vals, char *nm)
Check validity of 1-letter string from a set of possible values (typically used in S4 validity method...
Definition: Mutils.c:254
#define dngettext(pkg, String, StringP, N)
Definition: Mutils.h:34
SEXP set_double_by_name(SEXP obj, double val, char *nm)
Definition: Mutils.c:60
char La_rcond_type(const char *typstr)
Definition: Mutils.c:27
Definition: Mutils.h:181
#define do_ii_FILL(_i_, _j_)
#define SET_packed_setDiag
SEXP set_factors(SEXP obj, SEXP val, char *nm)
Caches 'val' in the 'factors' slot of obj, i.e.
Definition: Mutils.c:130
SEXP get_factors(SEXP obj, char *nm)
Definition: Mutils.c:106
#define MAKE_SYMMETRIC_BODY(_TO_, _FROM_)
Definition: Mutils.c:220
#define TRUE_
Definition: Mutils.c:1174
#define UPP
Definition: Mutils.h:81
#define NON_TRIVIAL_DN
#define MATRIX_VALID_ndense
Definition: Mutils.h:374
static R_INLINE int * expand_cmprPt(int ncol, const int mp[], int mj[])
Expand compressed pointers in the array mp into a full set of indices in the array mj...
Definition: Mutils.h:259
SEXP m_encodeInd2(SEXP i, SEXP j, SEXP di, SEXP orig_1, SEXP chk_bnds)
Encode Matrix index (i,j) |–> i + j * nrow {i,j : 0-origin}.
Definition: Mutils.c:945
#define AZERO(x, n)
Definition: Mutils.h:140
#define MATRIX_VALID_ddense
Definition: Mutils.h:360
#define uplo_P(_x_)
Definition: Mutils.h:174
SEXP dup_mMatrix_as_dgeMatrix2(SEXP A, Rboolean tr_if_vec)
Definition: Mutils.c:821
#define _(String)
Definition: Mutils.h:32
#define DUP_MMATRIX_SET_1
SEXP dup_mMatrix_as_geMatrix(SEXP A)
Duplicate a [dln]denseMatrix or a numeric matrix or even a vector as a [dln]geMatrix.
Definition: Mutils.c:648
void SET_DimNames_symm(SEXP dest, SEXP src)
Set 'Dimnames' slot of 'dest' from the one of 'src' when 'src' is a "symmetricMatrix" with possibly a...
Definition: Mutils.c:1328
#define SPRINTF
SEXP d_packed_setDiag(double *diag, int l_d, SEXP x, int n)
Definition: Mutils.c:461
#define MAKE_TRIANGULAR_BODY(_TO_, _FROM_, _ZERO_, _ONE_)
Definition: Mutils.c:192
char La_norm_type(const char *typstr)
Definition: Mutils.c:8
SEXP new_dgeMatrix(int nrow, int ncol)
Definition: Mutils.c:856
void tr_d_packed_getDiag(double *dest, SEXP x, int n)
Definition: Mutils.c:551
#define LOW
Definition: Mutils.h:82
#define COPY_a_AND_b_j
SEXP dimNames_validate(SEXP obj)
Definition: Mutils.c:343
SEXP Matrix_getElement(SEXP list, char *nm)
Return the element of a given name from a named list.
Definition: Mutils.c:592
void make_d_matrix_symmetric(double *to, SEXP from) void make_i_matrix_symmetric(int *to
#define diag_P(_x_)
Definition: Mutils.h:175
Definition: Mutils.h:181
#define DUP_MMATRIX_NON_CLASS(transpose_if_vec)
SEXP Matrix_pSym
Definition: Syms.h:2
#define DUP_MMATRIX_ddense_CASES
SEXP symmetric_DimNames(SEXP dn)
Produce symmetric 'Dimnames' from possibly asymmetric ones.
Definition: Mutils.c:1269
SEXP tr_d_packed_setDiag(double *diag, int l_d, SEXP x, int n)
Definition: Mutils.c:505
SEXP Matrix_diagSym
Definition: Syms.h:2
SEXP Dim_validate(SEXP obj, SEXP name)
Definition: Mutils.c:337
SEXP R_rbind2_vector(SEXP a, SEXP b)
From the two 'x' slots of two dense matrices a and b, compute the 'x' slot of rbind(a, b)
Definition: Mutils.c:1120
Definition: Mutils.h:181
SEXP R_set_factors(SEXP obj, SEXP val, SEXP name, SEXP warn)
Definition: Mutils.c:166
static R_INLINE SEXP ALLOC_SLOT(SEXP obj, SEXP nm, SEXPTYPE type, int length)
Allocate an SEXP of given type and length, assign it as slot nm in the object, and return the SEXP...
Definition: Mutils.h:240
SEXP R_symmetric_Dimnames(SEXP x)
Even if the Dimnames slot is list(NULL, ) etc, return symmetric dimnames: Get ...
Definition: Mutils.c:1315
void d_packed_getDiag(double *dest, SEXP x, int n)
Copy the diagonal elements of the packed denseMatrix x to dest.
Definition: Mutils.c:433
SEXP dense_nonpacked_validate(SEXP obj)
Definition: Mutils.c:310
double get_double_by_name(SEXP obj, char *nm)
Definition: Mutils.c:44
SEXP l_packed_setDiag(int *diag, int l_d, SEXP x, int n)
Definition: Mutils.c:491
#define PACKED_TO_FULL(TYPE)
Definition: Mutils.c:371