Matrix  $Rev: 3071 $ at $LastChangedDate: 2015-03-26 15:35:47 +0100 (Thu, 26 Mar 2015) $
dtrMatrix.c
Go to the documentation of this file.
1 /* double (precision) TRiangular Matrices */
2 
3 #include "dtrMatrix.h"
4 
6 {
7  SEXP val = GET_SLOT(obj, Matrix_DimSym);
8 
9  if (LENGTH(val) < 2)
10  return mkString(_("'Dim' slot has length less than two"));
11  if (INTEGER(val)[0] != INTEGER(val)[1])
12  return mkString(_("Matrix is not square"));
13  if (isString(val = check_scalar_string(GET_SLOT(obj, Matrix_uploSym),
14  "LU", "uplo"))) return val;
15  if (isString(val = check_scalar_string(GET_SLOT(obj, Matrix_diagSym),
16  "NU", "diag"))) return val;
17  return ScalarLogical(1);
18 }
19 
20 static
21 double get_norm(SEXP obj, const char *typstr)
22 {
23  char typnm[] = {'\0', '\0'};
24  int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym));
25  double *work = (double *) NULL;
26 
27  typnm[0] = La_norm_type(typstr);
28  if (*typnm == 'I') {
29  work = (double *) R_alloc(dims[0], sizeof(double));
30  }
31  return F77_CALL(dlantr)(typnm, uplo_P(obj), diag_P(obj), dims, dims+1,
32  REAL(GET_SLOT(obj, Matrix_xSym)), dims, work);
33 }
34 
35 
36 SEXP dtrMatrix_norm(SEXP obj, SEXP type)
37 {
38  return ScalarReal(get_norm(obj, CHAR(asChar(type))));
39 }
40 
41 SEXP dtrMatrix_rcond(SEXP obj, SEXP type)
42 {
43  char typnm[] = {'\0', '\0'};
44  int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym)), info;
45  double rcond;
46 
47  typnm[0] = La_rcond_type(CHAR(asChar(type)));
48  F77_CALL(dtrcon)(typnm, uplo_P(obj), diag_P(obj), dims,
49  REAL(GET_SLOT(obj, Matrix_xSym)), dims, &rcond,
50  (double *) R_alloc(3*dims[0], sizeof(double)),
51  (int *) R_alloc(dims[0], sizeof(int)), &info);
52  return ScalarReal(rcond);
53 }
54 
55 SEXP dtrMatrix_solve(SEXP a)
56 {
57  SEXP val = PROTECT(duplicate(a));
58  int info, *Dim = INTEGER(GET_SLOT(val, Matrix_DimSym));
59  F77_CALL(dtrtri)(uplo_P(val), diag_P(val), Dim,
60  REAL(GET_SLOT(val, Matrix_xSym)), Dim, &info);
61  UNPROTECT(1);
62  return val;
63 }
64 
65 SEXP dtrMatrix_chol2inv(SEXP a)
66 {
67  SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dpoMatrix")));
68  int info, n;
69 
70  slot_dup(val, a, Matrix_DimSym);
71  slot_dup(val, a, Matrix_uploSym);
72  slot_dup(val, a, Matrix_diagSym);
74  slot_dup(val, a, Matrix_xSym);
75  n = *INTEGER(GET_SLOT(val, Matrix_DimSym));
76  F77_CALL(dpotri)(uplo_P(val), &n,
77  REAL(GET_SLOT(val, Matrix_xSym)), &n, &info);
78  UNPROTECT(1);
79  return val;
80 }
81 
82 SEXP dtrMatrix_matrix_solve(SEXP a, SEXP b)
83 {
84  SEXP ans = PROTECT(dup_mMatrix_as_dgeMatrix(b));
85  int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
86  *bdims = INTEGER(GET_SLOT(ans, Matrix_DimSym));
87  int n = bdims[0], nrhs = bdims[1];
88  double one = 1.0;
89 
90  if (adims[0] != n || n != adims[1])
91  error(_("Dimensions of system to be solved are inconsistent"));
92  F77_CALL(dtrsm)("L", uplo_P(a), "N", diag_P(a),
93  &n, &nrhs, &one, REAL(GET_SLOT(a, Matrix_xSym)), &n,
94  REAL(GET_SLOT(ans, Matrix_xSym)), &n);
95  UNPROTECT(1);
96  return ans;
97 }
98 
99 // to be used for all three: '%*%', crossprod() and tcrossprod()
110 SEXP dtrMatrix_matrix_mm(SEXP a, SEXP b, SEXP right, SEXP trans)
111 {
112  /* called from "%*%", crossprod() and tcrossprod() in ../R/products.R
113  *
114  * Because 'a' must be square, the size of the answer 'val',
115  * is the same as the size of 'b' */
116  SEXP val = PROTECT(dup_mMatrix_as_dgeMatrix(b));
117  int rt = asLogical(right); /* if(rt), compute b %*% op(a), else op(a) %*% b */
118  int tr = asLogical(trans);/* if true, use t(a) */
119  int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
120  *bdims = INTEGER(GET_SLOT(val, Matrix_DimSym));
121  int m = bdims[0], n = bdims[1];
122  double one = 1.;
123 
124  if (adims[0] != adims[1])
125  error(_("dtrMatrix must be square"));
126  if ((rt && adims[0] != n) || (!rt && adims[1] != m))
127  error(_("Matrices are not conformable for multiplication"));
128  if (m >= 1 && n >= 1) {
129  // Level 3 BLAS - DTRMM() --> see call further below
130  F77_CALL(dtrmm)(rt ? "R" : "L", uplo_P(a),
131  /*trans_A = */ tr ? "T" : "N",
132  diag_P(a), &m, &n, &one,
133  REAL(GET_SLOT(a, Matrix_xSym)), adims,
134  REAL(GET_SLOT(val, Matrix_xSym)), &m);
135  }
136 
137  SEXP
138  dn_a = GET_SLOT( a, Matrix_DimNamesSym),
139  dn = GET_SLOT(val, Matrix_DimNamesSym);
140  /* matrix product a %*% b, t(a) %*% b, b %*% a, or b %*% t(a)
141  * (right, trans) = (F, F) (F, T) (T, F) (T, T)
142  * set:from_a = 0:0 0:1 1:1 1:0
143  */
144  SET_VECTOR_ELT(dn, rt ? 1 : 0, VECTOR_ELT(dn_a, (rt+tr) % 2));
145 
146  UNPROTECT(1);
147  return val;
148 }
149 
160 SEXP dtrMatrix_dtrMatrix_mm(SEXP a, SEXP b, SEXP right, SEXP trans)
161 {
162  /* called from "%*%" : (x,y, FALSE,FALSE),
163  crossprod() : (x,y, FALSE, TRUE) , and
164  tcrossprod(): (y,x, TRUE , TRUE)
165  * -
166  * TWO cases : (1) result is triangular <=> uplo's "match" (i.e., non-equal iff trans)
167  * === (2) result is "general"
168  */
169  SEXP val,/* = in case (2): PROTECT(dup_mMatrix_as_dgeMatrix(b)); */
170  d_a = GET_SLOT(a, Matrix_DimSym),
171  uplo_a = GET_SLOT(a, Matrix_uploSym), diag_a = GET_SLOT(a, Matrix_diagSym),
172  uplo_b = GET_SLOT(b, Matrix_uploSym), diag_b = GET_SLOT(b, Matrix_diagSym);
173  int rt = asLogical(right);
174  int tr = asLogical(trans);
175  int *adims = INTEGER(d_a), n = adims[0];
176  double *valx = (double *) NULL /*Wall*/;
177  const char
178  *uplo_a_ch = CHAR(STRING_ELT(uplo_a, 0)), /* = uplo_P(a) */
179  *diag_a_ch = CHAR(STRING_ELT(diag_a, 0)), /* = diag_P(a) */
180  *uplo_b_ch = CHAR(STRING_ELT(uplo_b, 0)), /* = uplo_P(b) */
181  *diag_b_ch = CHAR(STRING_ELT(diag_b, 0)); /* = diag_P(b) */
182  Rboolean same_uplo = (*uplo_a_ch == *uplo_b_ch),
183  matching_uplo = tr ? (!same_uplo) : same_uplo,
184  uDiag_b = /* -Wall: */ FALSE;
185 
186  if (INTEGER(GET_SLOT(b, Matrix_DimSym))[0] != n)
187  /* validity checking already "assures" square matrices ... */
188  error(_("\"dtrMatrix\" objects in '%*%' must have matching (square) dimension"));
189  if(matching_uplo) {
190  /* ==> result is triangular -- "dtrMatrix" !
191  * val := dup_mMatrix_as_dtrMatrix(b) : */
192  int sz = n * n;
193  val = PROTECT(NEW_OBJECT(MAKE_CLASS("dtrMatrix")));
194  SET_SLOT(val, Matrix_uploSym, duplicate(uplo_b));
195  SET_SLOT(val, Matrix_DimSym, duplicate(d_a));
196  SET_DimNames(val, b);
197  valx = REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, sz));
198  Memcpy(valx, REAL(GET_SLOT(b, Matrix_xSym)), sz);
199  if((uDiag_b = (*diag_b_ch == 'U'))) {
200  /* unit-diagonal b - may contain garbage in diagonal */
201  for (int i = 0; i < n; i++)
202  valx[i * (n+1)] = 1.;
203  }
204  } else { /* different "uplo" ==> result is "dgeMatrix" ! */
205  val = PROTECT(dup_mMatrix_as_dgeMatrix(b));
206  SEXP
207  dn_a = GET_SLOT( a , Matrix_DimNamesSym),
208  dn = GET_SLOT(val, Matrix_DimNamesSym);
209  /* matrix product a %*% b, t(a) %*% b, b %*% a, or b %*% t(a)
210  * (right, trans) = (F, F) (F, T) (T, F) (T, T)
211  * set:from_a = 0:0 0:1 1:1 1:0
212  */
213  SET_VECTOR_ELT(dn, rt ? 1 : 0, VECTOR_ELT(dn_a, (rt+tr) % 2));
214  }
215  if (n >= 1) {
216  double alpha = 1.;
217  /* Level 3 BLAS - DTRMM(): Compute one of the matrix multiplication operations
218  * B := alpha*op( A )*B ["L"], or B := alpha*B*op( A ) ["R"],
219  * where trans_A determines op(A):= A "N"one or
220  * op(A):= t(A) "T"ransposed */
221  F77_CALL(dtrmm)(rt ? "R" : "L", uplo_a_ch,
222  /*trans_A = */ tr ? "T" : "N", diag_a_ch, &n, &n, &alpha,
223  REAL(GET_SLOT(a, Matrix_xSym)), adims,
224  REAL(GET_SLOT(val, Matrix_xSym)), &n);
225  }
226  if(matching_uplo) {
227  make_d_matrix_triangular(valx, tr ? b : a); /* set "other triangle" to 0 */
228  if(*diag_a_ch == 'U' && uDiag_b) /* result remains uni-diagonal */
229  SET_SLOT(val, Matrix_diagSym, duplicate(diag_a));
230  }
231  UNPROTECT(1);
232  return val;
233 }
234 
235 
236 SEXP dtrMatrix_as_matrix(SEXP from, SEXP keep_dimnames)
237 {
238  int *Dim = INTEGER(GET_SLOT(from, Matrix_DimSym));
239  int m = Dim[0], n = Dim[1];
240  SEXP val = PROTECT(allocMatrix(REALSXP, m, n));
241  make_d_matrix_triangular(Memcpy(REAL(val),
242  REAL(GET_SLOT(from, Matrix_xSym)), m * n),
243  from);
244  if(asLogical(keep_dimnames))
245  setAttrib(val, R_DimNamesSymbol, GET_SLOT(from, Matrix_DimNamesSym));
246  UNPROTECT(1);
247  return val;
248 }
249 
250 #define GET_trMatrix_Diag(_C_TYPE_, _SEXPTYPE_, _SEXP_, _ONE_) \
251  int i, n = INTEGER(GET_SLOT(x, Matrix_DimSym))[0]; \
252  SEXP x_x = GET_SLOT(x, Matrix_xSym); \
253  \
254  SEXP ret = PROTECT(allocVector(_SEXPTYPE_, n)); \
255  _C_TYPE_ *rv = _SEXP_(ret), \
256  *xv = _SEXP_(x_x); \
257  \
258  if ('U' == diag_P(x)[0]) { \
259  for (i = 0; i < n; i++) rv[i] = _ONE_; \
260  } else { \
261  for (i = 0; i < n; i++) rv[i] = xv[i * (n + 1)]; \
262  } \
263  UNPROTECT(1); \
264  return ret
265 
266 
267 SEXP dtrMatrix_getDiag(SEXP x) {
268  GET_trMatrix_Diag(double, REALSXP, REAL, 1.);
269 }
270 
271 SEXP ltrMatrix_getDiag(SEXP x) {
272  GET_trMatrix_Diag( int, LGLSXP, LOGICAL, 1);
273 }
274 
275 #define SET_trMatrix_Diag(_C_TYPE_, _SEXP_) \
276  if ('U' == diag_P(x)[0]) \
277  error(_("cannot set diag() as long as 'diag = \"U\"'")); \
278  /* careful to recycle RHS value: */ \
279  int n = INTEGER(GET_SLOT(x, Matrix_DimSym))[0]; \
280  int l_d = LENGTH(d); Rboolean d_full = (l_d == n); \
281  if (!d_full && l_d != 1) \
282  error(_("replacement diagonal has wrong length")); \
283  SEXP ret = PROTECT(duplicate(x)), \
284  r_x = GET_SLOT(ret, Matrix_xSym); \
285  _C_TYPE_ *dv = _SEXP_(d), \
286  *rv = _SEXP_(r_x); \
287  \
288  if(d_full) for (int i = 0; i < n; i++) \
289  rv[i * (n + 1)] = dv[i]; \
290  else for (int i = 0; i < n; i++) \
291  rv[i * (n + 1)] = *dv; \
292  \
293  UNPROTECT(1); \
294  return ret
295 
296 SEXP dtrMatrix_setDiag(SEXP x, SEXP d) {
297  SET_trMatrix_Diag(double, REAL);
298 }
299 
300 SEXP ltrMatrix_setDiag(SEXP x, SEXP d) {
301  SET_trMatrix_Diag( int, LOGICAL);
302 }
303 
304 SEXP dtrMatrix_addDiag(SEXP x, SEXP d) {
305  int n = INTEGER(GET_SLOT(x, Matrix_DimSym))[0];
306  SEXP ret = PROTECT(duplicate(x)),
307  r_x = GET_SLOT(ret, Matrix_xSym);
308  double *dv = REAL(d), *rv = REAL(r_x);
309 
310  if ('U' == diag_P(x)[0])
311  error(_("cannot add diag() as long as 'diag = \"U\"'"));
312  for (int i = 0; i < n; i++) rv[i * (n + 1)] += dv[i];
313 
314  UNPROTECT(1);
315  return ret;
316 }
317 
318 
319 SEXP dtrMatrix_as_dtpMatrix(SEXP from)
320 {
321  SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dtpMatrix"))),
322  uplo = GET_SLOT(from, Matrix_uploSym),
323  diag = GET_SLOT(from, Matrix_diagSym),
324  dimP = GET_SLOT(from, Matrix_DimSym);
325  int n = *INTEGER(dimP);
326 
327  SET_SLOT(val, Matrix_DimSym, duplicate(dimP));
328  SET_SLOT(val, Matrix_diagSym, duplicate(diag));
329  SET_SLOT(val, Matrix_uploSym, duplicate(uplo));
330  full_to_packed_double(
331  REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, (n*(n+1))/2)),
332  REAL(GET_SLOT(from, Matrix_xSym)), n,
333  *CHAR(STRING_ELT(uplo, 0)) == 'U' ? UPP : LOW,
334  *CHAR(STRING_ELT(diag, 0)) == 'U' ? UNT : NUN);
335  SET_SLOT(val, Matrix_DimNamesSym,
336  duplicate(GET_SLOT(from, Matrix_DimNamesSym)));
337  UNPROTECT(1);
338  return val;
339 }
SEXP Matrix_DimSym
Definition: Syms.h:2
static double get_norm(SEXP obj, const char *typstr)
Definition: dtrMatrix.c:21
SEXP ltrMatrix_getDiag(SEXP x)
Definition: dtrMatrix.c:271
SEXP dtrMatrix_as_matrix(SEXP from, SEXP keep_dimnames)
Definition: dtrMatrix.c:236
SEXP Matrix_xSym
Definition: Syms.h:2
#define slot_dup(dest, src, sym)
Definition: Mutils.h:149
SEXP Matrix_DimNamesSym
Definition: Syms.h:2
SEXP dtrMatrix_chol2inv(SEXP a)
Definition: dtrMatrix.c:65
SEXP Matrix_uploSym
Definition: Syms.h:2
SEXP dtrMatrix_norm(SEXP obj, SEXP type)
Definition: dtrMatrix.c:36
SEXP dup_mMatrix_as_dgeMatrix(SEXP A)
Definition: Mutils.c:852
void make_d_matrix_triangular(double *to, SEXP from) void make_i_matrix_triangular(int *to
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
SEXP dtrMatrix_getDiag(SEXP x)
Definition: dtrMatrix.c:267
SEXP ltrMatrix_setDiag(SEXP x, SEXP d)
Definition: dtrMatrix.c:300
char La_rcond_type(const char *typstr)
Definition: Mutils.c:27
SEXP dtrMatrix_matrix_solve(SEXP a, SEXP b)
Definition: dtrMatrix.c:82
SEXP dtrMatrix_rcond(SEXP obj, SEXP type)
Definition: dtrMatrix.c:41
#define UPP
Definition: Mutils.h:81
SEXP dtrMatrix_as_dtpMatrix(SEXP from)
Definition: dtrMatrix.c:319
SEXP dtrMatrix_setDiag(SEXP x, SEXP d)
Definition: dtrMatrix.c:296
#define uplo_P(_x_)
Definition: Mutils.h:174
SEXP dtrMatrix_addDiag(SEXP x, SEXP d)
Definition: dtrMatrix.c:304
#define _(String)
Definition: Mutils.h:32
SEXP dtrMatrix_matrix_mm(SEXP a, SEXP b, SEXP right, SEXP trans)
Matrix products dense triangular Matrices o
Definition: dtrMatrix.c:110
SEXP dtrMatrix_dtrMatrix_mm(SEXP a, SEXP b, SEXP right, SEXP trans)
Matrix products of dense triangular Matrices.
Definition: dtrMatrix.c:160
SEXP triangularMatrix_validate(SEXP obj)
Definition: dtrMatrix.c:5
#define SET_trMatrix_Diag(_C_TYPE_, _SEXP_)
Definition: dtrMatrix.c:275
char La_norm_type(const char *typstr)
Definition: Mutils.c:8
#define LOW
Definition: Mutils.h:82
#define GET_trMatrix_Diag(_C_TYPE_, _SEXPTYPE_, _SEXP_, _ONE_)
Definition: dtrMatrix.c:250
#define diag_P(_x_)
Definition: Mutils.h:175
#define UNT
Definition: Mutils.h:84
SEXP Matrix_diagSym
Definition: Syms.h:2
#define NUN
Definition: Mutils.h:83
static R_INLINE void SET_DimNames(SEXP dest, SEXP src)
Definition: Mutils.h:161
SEXP dtrMatrix_solve(SEXP a)
Definition: dtrMatrix.c:55
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