Matrix  $Rev: 3071 $ at $LastChangedDate: 2015-03-26 15:35:47 +0100 (Thu, 26 Mar 2015) $
dgTMatrix.c
Go to the documentation of this file.
1 #include <Rinternals.h>
2 /* for R_LEN... */
3 
4 #include "dgTMatrix.h"
5 
6 #include "chm_common.h"
7 #include "Tsparse.h"
8 
9 SEXP xTMatrix_validate(SEXP x)
10 {
11  /* Almost everything now in Tsparse_validate ( ./Tsparse.c )
12  * *but* the checking of the 'x' slot : */
13  if (LENGTH(GET_SLOT(x, Matrix_iSym)) !=
14  LENGTH(GET_SLOT(x, Matrix_xSym)))
15  return mkString(_("lengths of slots i and x must match"));
16  return ScalarLogical(1);
17 }
18 
19 static void
20 d_insert_triplets_in_array(int m, int n, int nnz,
21  const int xi[], const int xj[], const double xx[],
22  /* --> */ double vx[])
23 {
24  // For ( m*n ) > INT_MAX, we here assume that size_t is using 64-bit !
25  size_t m_ = (size_t) m, len = sizeof(double) * m_ * n;
26  if(len == sizeof(double) * (double)m_ *n)
27  memset(vx, 0, len);
28  else { // len did overflow -- this should call memset() several times:
29  size_t max_l = (1 << (sizeof(size_t)-1)); // = 2^(N-1)
30  max_l += ((long)max_l - 1); // = 2^(N-1) + 2^(N-1) - 1 = 2^N - 1
31  double dlen = ((double)m_) * n;
32  if(dlen > max_l)
33  error(_("too large matrix: %.0f"), dlen);
34  // else : m * n does fit -- call memset() several times:
35  dlen *= sizeof(double);
36  memset(vx, 0, max_l);
37  double len_done = 0.; // length in bytes // but also need length in double
38  while((len_done += max_l) < dlen) {
39  size_t dd = (dlen - len_done < max_l) ? (size_t)(dlen - len_done) : max_l;
40  memset(vx + (int)(len_done/sizeof(double)), 0, dd);
41  }
42  }
43  for (int i = 0; i < nnz; i++) {
44  vx[xi[i] + xj[i] * m_] += xx[i]; /* allow redundant entries in x */
45  }
46 }
47 
48 static void
49 l_insert_triplets_in_array(int m, int n, int nnz,
50  const int xi[], const int xj[], const int xx[],
51  /* --> */ int vx[])
52 {
53  // For ( m*n ) > INT_MAX, we here assume that size_t is using 64-bit !
54  size_t m_ = (size_t) m, len = sizeof(int) * m_ * n;
55  if(len == sizeof(int) * (double)m_ *n)
56  memset(vx, 0, len);
57  else { // len did overflow -- this should call memset() several times:
58  size_t max_l = (1 << (sizeof(size_t)-1)); // = 2^(N-1)
59  max_l += ((long)max_l - 1); // = 2^(N-1) + 2^(N-1) - 1 = 2^N - 1
60  double dlen = ((double)m_) * n;
61  if(dlen > max_l)
62  error(_("too large matrix: %.0f"), dlen);
63  // else : m * n does fit -- call memset() several times:
64  dlen *= sizeof(int);
65  memset(vx, 0, max_l);
66  double len_done = 0.; // length in bytes // but also need length in int
67  while((len_done += max_l) < dlen) {
68  size_t dd = (dlen - len_done < max_l) ? (size_t)(dlen - len_done) : max_l;
69  memset(vx + (int)(len_done/sizeof(int)), 0, dd);
70  }
71  }
72 
73  for (int i = 0; i < nnz; i++) {
74  size_t ind = xi[i] + xj[i] * m_;
75  if(vx[ind] == NA_LOGICAL) {
76  // do nothing: remains NA
77  } else if(xx[i] == NA_LOGICAL)
78  vx[ind] = NA_LOGICAL;
79  else // "or" :
80  vx[ind] |= xx[i];
81  }
82 }
83 
84 #define MAKE_gTMatrix_to_geMatrix(_t1_, _SEXPTYPE_, _SEXP_) \
85 SEXP _t1_ ## gTMatrix_to_ ## _t1_ ## geMatrix(SEXP x) \
86 { \
87  SEXP dd = GET_SLOT(x, Matrix_DimSym), \
88  islot = GET_SLOT(x, Matrix_iSym), \
89  ans = PROTECT(NEW_OBJECT(MAKE_CLASS(#_t1_ "geMatrix"))); \
90  \
91  int *dims = INTEGER(dd), \
92  m = dims[0], \
93  n = dims[1]; \
94  double len = m * (double)n; \
95  \
96  if (len > R_XLEN_T_MAX) \
97  error(_("Cannot coerce to too large *geMatrix with %.0f entries"), \
98  len); \
99  \
100  SET_SLOT(ans, Matrix_factorSym, allocVector(VECSXP, 0)); \
101  SET_SLOT(ans, Matrix_DimSym, duplicate(dd)); \
102  SET_DimNames(ans, x); \
103  SET_SLOT(ans, Matrix_xSym, allocVector(_SEXPTYPE_, (R_xlen_t)len)); \
104  _t1_ ## _insert_triplets_in_array(m, n, length(islot), \
105  INTEGER(islot), \
106  INTEGER(GET_SLOT(x, Matrix_jSym)),\
107  _SEXP_(GET_SLOT(x, Matrix_xSym)), \
108  _SEXP_(GET_SLOT(ans, Matrix_xSym))); \
109  UNPROTECT(1); \
110  return ans; \
111 }
112 
113 MAKE_gTMatrix_to_geMatrix(d, REALSXP, REAL)
114 
115 MAKE_gTMatrix_to_geMatrix(l, LGLSXP, LOGICAL)
116 
117 #undef MAKE_gTMatrix_to_geMatrix
118 
119 #define MAKE_gTMatrix_to_matrix(_t1_, _SEXPTYPE_, _SEXP_) \
120 SEXP _t1_ ## gTMatrix_to_matrix(SEXP x) \
121 { \
122  SEXP dd = GET_SLOT(x, Matrix_DimSym), \
123  dn = GET_SLOT(x, Matrix_DimNamesSym), \
124  islot = GET_SLOT(x, Matrix_iSym); \
125  int m = INTEGER(dd)[0], \
126  n = INTEGER(dd)[1]; \
127  SEXP ans = PROTECT(allocMatrix(_SEXPTYPE_, m, n)); \
128  if(VECTOR_ELT(dn, 0) != R_NilValue || VECTOR_ELT(dn, 1) != R_NilValue) \
129  /* matrix() with non-trivial dimnames */ \
130  setAttrib(ans, R_DimNamesSymbol, duplicate(dn)); \
131  _t1_ ## _insert_triplets_in_array(m, n, length(islot), \
132  INTEGER(islot), \
133  INTEGER(GET_SLOT(x, Matrix_jSym)),\
134  _SEXP_(GET_SLOT(x, Matrix_xSym)), \
135  _SEXP_(ans)); \
136  UNPROTECT(1); \
137  return ans; \
138 }
139 
140 MAKE_gTMatrix_to_matrix(d, REALSXP, REAL)
141 
142 MAKE_gTMatrix_to_matrix(l, LGLSXP, LOGICAL)
143 
144 #undef MAKE_gTMatrix_to_matrix
SEXP Matrix_xSym
Definition: Syms.h:2
#define MAKE_gTMatrix_to_matrix(_t1_, _SEXPTYPE_, _SEXP_)
#define MAKE_gTMatrix_to_geMatrix(_t1_, _SEXPTYPE_, _SEXP_)
Definition: dgTMatrix.c:84
static void l_insert_triplets_in_array(int m, int n, int nnz, const int xi[], const int xj[], const int xx[], int vx[])
Definition: dgTMatrix.c:49
static void d_insert_triplets_in_array(int m, int n, int nnz, const int xi[], const int xj[], const double xx[], double vx[])
Definition: dgTMatrix.c:20
#define _(String)
Definition: Mutils.h:32
SEXP Matrix_iSym
Definition: Syms.h:2
SEXP xTMatrix_validate(SEXP x)
Definition: dgTMatrix.c:9