Matrix  $Rev: 3071 $ at $LastChangedDate: 2015-03-26 15:35:47 +0100 (Thu, 26 Mar 2015) $
dgeMatrix.c
Go to the documentation of this file.
1 #include "dgeMatrix.h"
2 // -> Mutils.h etc
3 
4 SEXP dMatrix_validate(SEXP obj)
5 {
6  SEXP x = GET_SLOT(obj, Matrix_xSym),
7  Dim = GET_SLOT(obj, Matrix_DimSym);
8  if (!isReal(x))
9  return mkString(_("x slot must be numeric \"double\""));
10  SEXP val;
11  if (isString(val = dim_validate(Dim, "Matrix")))
12  return val;
13  return ScalarLogical(1);
14 }
15 
16 SEXP dgeMatrix_validate(SEXP obj)
17 {
18  SEXP val;
19  if (isString(val = dim_validate(GET_SLOT(obj, Matrix_DimSym), "dgeMatrix")))
20  return(val);
21  if (isString(val = dense_nonpacked_validate(obj)))
22  return(val);
23  SEXP fact = GET_SLOT(obj, Matrix_factorSym);
24  if (length(fact) > 0 && getAttrib(fact, R_NamesSymbol) == R_NilValue)
25  return mkString(_("factors slot must be named list"));
26  return ScalarLogical(1);
27 }
28 
29 static
30 double get_norm(SEXP obj, const char *typstr)
31 {
32  if(any_NA_in_x(obj))
33  return NA_REAL;
34  else {
35  char typnm[] = {'\0', '\0'};
36  int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym));
37  double *work = (double *) NULL;
38 
39  typnm[0] = La_norm_type(typstr);
40  if (*typnm == 'I') {
41  work = (double *) R_alloc(dims[0], sizeof(double));
42  }
43  return F77_CALL(dlange)(typstr, dims, dims+1,
44  REAL(GET_SLOT(obj, Matrix_xSym)),
45  dims, work);
46  }
47 }
48 
49 SEXP dgeMatrix_norm(SEXP obj, SEXP type)
50 {
51  return ScalarReal(get_norm(obj, CHAR(asChar(type))));
52 }
53 
54 SEXP dgeMatrix_rcond(SEXP obj, SEXP type)
55 {
56  SEXP LU = PROTECT(dgeMatrix_LU_(obj, FALSE));/* <- not warning about singularity */
57  char typnm[] = {'\0', '\0'};
58  int *dims = INTEGER(GET_SLOT(LU, Matrix_DimSym)), info;
59  double anorm, rcond;
60 
61  if (dims[0] != dims[1] || dims[0] < 1) {
62  UNPROTECT(1);
63  error(_("rcond requires a square, non-empty matrix"));
64  }
65  typnm[0] = La_rcond_type(CHAR(asChar(type)));
66  anorm = get_norm(obj, typnm);
67  F77_CALL(dgecon)(typnm,
68  dims, REAL(GET_SLOT(LU, Matrix_xSym)),
69  dims, &anorm, &rcond,
70  (double *) R_alloc(4*dims[0], sizeof(double)),
71  (int *) R_alloc(dims[0], sizeof(int)), &info);
72  UNPROTECT(1);
73  return ScalarReal(rcond);
74 }
75 
76 SEXP dgeMatrix_crossprod(SEXP x, SEXP trans)
77 {
78 #define DGE_CROSS_1 \
79  int tr = asLogical(trans);/* trans=TRUE: tcrossprod(x) */ \
80  SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dpoMatrix"))), \
81  nms = VECTOR_ELT(GET_SLOT(x, Matrix_DimNamesSym), tr ? 0 : 1), \
82  vDnms = ALLOC_SLOT(val, Matrix_DimNamesSym, VECSXP, 2); \
83  int *Dims = INTEGER(GET_SLOT(x, Matrix_DimSym)), \
84  *vDims = INTEGER(ALLOC_SLOT(val, Matrix_DimSym, INTSXP, 2)); \
85  int k = tr ? Dims[1] : Dims[0], \
86  n = tr ? Dims[0] : Dims[1]; \
87  double *vx = REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, n * n)), \
88  one = 1.0, zero = 0.0; \
89  \
90  Memzero(vx, n * n); \
91  SET_SLOT(val, Matrix_uploSym, mkString("U")); \
92  ALLOC_SLOT(val, Matrix_factorSym, VECSXP, 0); \
93  vDims[0] = vDims[1] = n; \
94  SET_VECTOR_ELT(vDnms, 0, duplicate(nms)); \
95  SET_VECTOR_ELT(vDnms, 1, duplicate(nms))
96 
97 #define DGE_CROSS_DO(_X_X_) \
98  if(n) \
99  F77_CALL(dsyrk)("U", tr ? "N" : "T", &n, &k, &one, \
100  _X_X_, Dims, &zero, vx, &n); \
101  UNPROTECT(1); \
102  return val
103 
104  DGE_CROSS_1;
105  DGE_CROSS_DO(REAL(GET_SLOT(x, Matrix_xSym)));
106 }
107 
108 
109 double* gematrix_real_x(SEXP x, int nn) {
110  if(class_P(x)[0] == 'd') // <<- FIXME: use R_check_class_etc(x, valid) !!!
111  return REAL(GET_SLOT(x, Matrix_xSym));
112 #ifdef _potentically_more_efficient_but_not_working
113  // else : 'l' or 'n' (for now !!)
114  int *xi = INTEGER(GET_SLOT(x, Matrix_xSym));
115  double *x_x;
116  C_or_Alloca_TO(x_x, nn, double);
117  for(int i=0; i < nn; i++)
118  x_x[i] = (double) xi[i];
119 
120  // FIXME: this is not possible either; the *caller* would have to Free(.)
121  if(nn >= SMALL_4_Alloca) Free(x_x);
122 #else
123  // ideally should be PROTECT()ed ==> make sure R does not run gc() now!
124  double *x_x = REAL(coerceVector(GET_SLOT(x, Matrix_xSym), REALSXP));
125 #endif
126  return x_x;
127 }
128 
130 SEXP _geMatrix_crossprod(SEXP x, SEXP trans)
131 {
132  DGE_CROSS_1;
133  double *x_x = gematrix_real_x(x, k * n);
134  DGE_CROSS_DO(x_x);
135 }
136 
137 SEXP geMatrix_crossprod(SEXP x, SEXP trans)
138 {
139  SEXP y = PROTECT(dup_mMatrix_as_geMatrix(x)),
140  val = _geMatrix_crossprod(y, trans);
141  UNPROTECT(1);
142  return val;
143 }
144 
145 SEXP dgeMatrix_dgeMatrix_crossprod(SEXP x, SEXP y, SEXP trans)
146 {
147 #define DGE_DGE_CROSS_1 \
148  int tr = asLogical(trans);/* trans=TRUE: tcrossprod(x,y) */ \
149  SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix"))), \
150  dn = PROTECT(allocVector(VECSXP, 2)); \
151  int *xDims = INTEGER(GET_SLOT(x, Matrix_DimSym)), \
152  *yDims = INTEGER(GET_SLOT(y, Matrix_DimSym)), \
153  *vDims; \
154  int m = xDims[!tr], n = yDims[!tr];/* -> result dim */ \
155  int xd = xDims[ tr], yd = yDims[ tr];/* the conformable dims */ \
156  double one = 1.0, zero = 0.0; \
157  \
158  if (xd != yd) \
159  error(_("Dimensions of x and y are not compatible for %s"), \
160  tr ? "tcrossprod" : "crossprod"); \
161  SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0)); \
162  /* establish dimnames */ \
163  SET_VECTOR_ELT(dn, 0, \
164  duplicate(VECTOR_ELT(GET_SLOT(x, Matrix_DimNamesSym), \
165  tr ? 0 : 1))); \
166  SET_VECTOR_ELT(dn, 1, \
167  duplicate(VECTOR_ELT(GET_SLOT(y, Matrix_DimNamesSym), \
168  tr ? 0 : 1))); \
169  SET_SLOT(val, Matrix_DimNamesSym, dn); \
170  vDims = INTEGER(ALLOC_SLOT(val, Matrix_DimSym, INTSXP, 2)); \
171  vDims[0] = m; vDims[1] = n; \
172  double *v = REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, m * n))
173 
174 #define DGE_DGE_CROSS_DO(_X_X_, _Y_Y_) \
175  if (xd > 0 && n > 0 && m > 0) \
176  F77_CALL(dgemm)(tr ? "N" : "T", tr ? "T" : "N", &m, &n, &xd, &one, \
177  _X_X_, xDims, \
178  _Y_Y_, yDims, &zero, v, &m); \
179  else \
180  Memzero(v, m * n); \
181  UNPROTECT(2); \
182  return val
183 
185  DGE_DGE_CROSS_DO(REAL(GET_SLOT(x, Matrix_xSym)),
186  REAL(GET_SLOT(y, Matrix_xSym)));
187 }
188 
190 SEXP _geMatrix__geMatrix_crossprod(SEXP x, SEXP y, SEXP trans)
191 {
193 
194  double *x_x = gematrix_real_x(x, m * xd);
195  double *y_x = gematrix_real_x(y, n * yd);
196 
197  DGE_DGE_CROSS_DO(x_x, y_x);
198 }
199 #undef DGE_DGE_CROSS_1
200 #undef DGE_DGE_CROSS_DO
201 
202 SEXP geMatrix_geMatrix_crossprod(SEXP x, SEXP y, SEXP trans)
203 {
204  SEXP gx = PROTECT(dup_mMatrix_as_geMatrix(x)),
205  gy = PROTECT(dup_mMatrix_as_geMatrix(y)),
206  val = _geMatrix__geMatrix_crossprod(gx, gy, trans);
207  UNPROTECT(2);
208  return val;
209 }
210 
211 SEXP dgeMatrix_matrix_crossprod(SEXP x, SEXP y, SEXP trans)
212 {
213 #define DGE_MAT_CROSS_1 \
214  int tr = asLogical(trans);/* trans=TRUE: tcrossprod(x,y) */ \
215  SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix"))), \
216  dn = PROTECT(allocVector(VECSXP, 2)), \
217  yDnms = R_NilValue, yD; \
218  int *xDims = INTEGER(GET_SLOT(x, Matrix_DimSym)), \
219  *yDims, *vDims, nprot = 2; \
220  int m = xDims[!tr], \
221  xd = xDims[ tr]; \
222  double one = 1.0, zero = 0.0; \
223  Rboolean y_has_dimNames; \
224  \
225  if (!isReal(y)) { \
226  if(isInteger(y) || isLogical(y)) { \
227  y = PROTECT(coerceVector(y, REALSXP)); \
228  nprot++; \
229  } \
230  else \
231  error(_("Argument y must be numeric, integer or logical")); \
232  } \
233  if(isMatrix(y)) { \
234  yDims = INTEGER(getAttrib(y, R_DimSymbol)); \
235  yDnms = getAttrib(y, R_DimNamesSymbol); \
236  y_has_dimNames = yDnms != R_NilValue; \
237  } else { /* ! matrix */ \
238  yDims = INTEGER(yD = PROTECT(allocVector(INTSXP, 2))); nprot++; \
239  if(xDims[0] == 1) { \
240  /* "new" (2014-10-10): "be tolerant" as for R 3.2.0*/ \
241  yDims[0] = 1; \
242  yDims[1] = LENGTH(y); \
243  } else { \
244  yDims[0] = LENGTH(y); \
245  yDims[1] = 1; \
246  } \
247  y_has_dimNames = FALSE; \
248  } \
249  int n = yDims[!tr],/* (m,n) -> result dim */ \
250  yd = yDims[ tr];/* (xd,yd): the conformable dims */ \
251  if (xd != yd) \
252  error(_("Dimensions of x and y are not compatible for %s"), \
253  tr ? "tcrossprod" : "crossprod"); \
254  SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0)); \
255  vDims = INTEGER(ALLOC_SLOT(val, Matrix_DimSym, INTSXP, 2)); \
256  vDims[0] = m; vDims[1] = n; \
257  /* establish dimnames */ \
258  SET_VECTOR_ELT(dn, 0, \
259  duplicate(VECTOR_ELT(GET_SLOT(x, Matrix_DimNamesSym), \
260  tr ? 0 : 1))); \
261  if(y_has_dimNames) \
262  SET_VECTOR_ELT(dn, 1, \
263  duplicate(VECTOR_ELT(yDnms, tr ? 0 : 1))); \
264  SET_SLOT(val, Matrix_DimNamesSym, dn); \
265  \
266  double *v = REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, m * n))
267 
268 #define DGE_MAT_CROSS_DO(_X_X_) \
269  if (xd > 0 && n > 0 && m > 0) \
270  F77_CALL(dgemm)(tr ? "N" : "T", tr ? "T" : "N", &m, &n, &xd, &one, \
271  _X_X_, xDims, REAL(y), yDims, \
272  &zero, v, &m); \
273  else \
274  Memzero(v, m * n); \
275  UNPROTECT(nprot); \
276  return val
277 
279  DGE_MAT_CROSS_DO(REAL(GET_SLOT(x, Matrix_xSym)));
280 }
281 
283 SEXP _geMatrix_matrix_crossprod(SEXP x, SEXP y, SEXP trans) {
285 
286  double *x_x = gematrix_real_x(x, m * xd);
287 
288  DGE_MAT_CROSS_DO(x_x);
289 }
290 
291 SEXP geMatrix_matrix_crossprod(SEXP x, SEXP y, SEXP trans) {
292  SEXP dx = PROTECT(dup_mMatrix_as_geMatrix(x)),
293  val = _geMatrix_matrix_crossprod(dx, y, trans);
294  UNPROTECT(1);
295  return val;
296 }
297 
298 // right = TRUE: %*% is called as *(y, x, right=TRUE)
299 SEXP dgeMatrix_matrix_mm(SEXP a, SEXP bP, SEXP right)
300 {
301 #define DGE_MAT_MM_1(N_PROT) \
302  SEXP val= PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix"))), \
303  dn = PROTECT(allocVector(VECSXP, 2)); \
304  int nprot = N_PROT + 2, \
305  *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)), \
306  *bdims = INTEGER(GET_SLOT(b, Matrix_DimSym)), \
307  *cdims = INTEGER(ALLOC_SLOT(val, Matrix_DimSym, INTSXP, 2)), \
308  Rt = asLogical(right), m, k, n; \
309  double one = 1., zero = 0.; \
310  \
311  if (Rt) { /* b %*% a : (m x k) (k x n) -> (m x n) */ \
312  m = bdims[0]; k = bdims[1]; n = adims[1]; \
313  if (adims[0] != k) \
314  error(_("Matrices are not conformable for multiplication")); \
315  } else { /* a %*% b : (m x k) (k x n) -> (m x n) */ \
316  m = adims[0]; k = adims[1]; n = bdims[1]; \
317  if (bdims[0] != k) \
318  error(_("Matrices are not conformable for multiplication")); \
319  } \
320  \
321  cdims[0] = m; cdims[1] = n; \
322  /* establish dimnames */ \
323  SET_VECTOR_ELT(dn, 0, duplicate( \
324  VECTOR_ELT(GET_SLOT(Rt ? b : a, \
325  Matrix_DimNamesSym), 0))); \
326  SET_VECTOR_ELT(dn, 1, \
327  duplicate( \
328  VECTOR_ELT(GET_SLOT(Rt ? a : b, \
329  Matrix_DimNamesSym), 1))); \
330  SET_SLOT(val, Matrix_DimNamesSym, dn); \
331  double *v = REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, m * n))
332 
333 #define DGE_MAT_MM_DO(_A_X_, _B_X_) \
334  if (m < 1 || n < 1 || k < 1) {/* zero extent matrices should work */ \
335  Memzero(v, m * n); \
336  } else { \
337  if (Rt) { /* b %*% a */ \
338  F77_CALL(dgemm) ("N", "N", &m, &n, &k, &one, \
339  _B_X_, &m, _A_X_, &k, &zero, v, &m); \
340  } else { /* a %*% b */ \
341  F77_CALL(dgemm) ("N", "N", &m, &n, &k, &one, \
342  _A_X_, &m, _B_X_, &k, &zero, v, &m); \
343  } \
344  } \
345  UNPROTECT(nprot); \
346  return val
347 
348  SEXP b = PROTECT(mMatrix_as_dgeMatrix(bP));
349  DGE_MAT_MM_1(1);
350  DGE_MAT_MM_DO(REAL(GET_SLOT(a, Matrix_xSym)),
351  REAL(GET_SLOT(b, Matrix_xSym)));
352 }
353 
355 SEXP _geMatrix_matrix_mm(SEXP a, SEXP b, SEXP right) {
356  DGE_MAT_MM_1(0);
357  double *a_x = gematrix_real_x(a, k * (Rt ? n : m));
358  double *b_x = gematrix_real_x(b, k * (Rt ? m : n));
359  DGE_MAT_MM_DO(a_x, b_x);
360 }
361 
363 SEXP geMatrix_matrix_mm(SEXP a, SEXP b, SEXP right) {
364  SEXP
365  da = PROTECT(dup_mMatrix_as_geMatrix(a)),
366  db = PROTECT(dup_mMatrix_as_geMatrix(b)),
367  val = _geMatrix_matrix_mm(da, db, right);
368  UNPROTECT(2);
369  return val;
370 }
371 
372 //---------------------------------------------------------------------
373 
374 SEXP dgeMatrix_getDiag(SEXP x)
375 {
376 #define geMatrix_getDiag_1 \
377  int *dims = INTEGER(GET_SLOT(x, Matrix_DimSym)); \
378  int i, m = dims[0], nret = (m < dims[1]) ? m : dims[1]; \
379  SEXP x_x = GET_SLOT(x, Matrix_xSym)
380 
382  SEXP ret = PROTECT(allocVector(REALSXP, nret));
383  double *rv = REAL(ret),
384  *xv = REAL(x_x);
385 
386 #define geMatrix_getDiag_2 \
387  for (i = 0; i < nret; i++) { \
388  rv[i] = xv[i * (m + 1)]; \
389  } \
390  UNPROTECT(1); \
391  return ret
392 
394 }
395 
396 SEXP lgeMatrix_getDiag(SEXP x)
397 {
399 
400  SEXP ret = PROTECT(allocVector(LGLSXP, nret));
401  int *rv = LOGICAL(ret),
402  *xv = LOGICAL(x_x);
403 
405 }
406 
407 #undef geMatrix_getDiag_1
408 #undef geMatrix_getDiag_2
409 
410 
411 SEXP dgeMatrix_setDiag(SEXP x, SEXP d)
412 {
413 #define geMatrix_setDiag_1 \
414  int *dims = INTEGER(GET_SLOT(x, Matrix_DimSym)); \
415  int m = dims[0], nret = (m < dims[1]) ? m : dims[1]; \
416  SEXP ret = PROTECT(duplicate(x)); \
417  SEXP r_x = GET_SLOT(ret, Matrix_xSym); \
418  int l_d = LENGTH(d); Rboolean d_full = (l_d == nret); \
419  if (!d_full && l_d != 1) \
420  error(_("replacement diagonal has wrong length"))
421 
423  double *dv = REAL(d), *rv = REAL(r_x);
424 
425 #define geMatrix_setDiag_2 \
426  if(d_full) for (int i = 0; i < nret; i++) \
427  rv[i * (m + 1)] = dv[i]; \
428  else for (int i = 0; i < nret; i++) \
429  rv[i * (m + 1)] = *dv; \
430  UNPROTECT(1); \
431  return ret
432 
434 }
435 
436 SEXP lgeMatrix_setDiag(SEXP x, SEXP d)
437 {
439  int *dv = INTEGER(d), *rv = INTEGER(r_x);
440 
442 }
443 
444 #undef geMatrix_setDiag_1
445 #undef geMatrix_setDiag_2
446 
447 SEXP dgeMatrix_addDiag(SEXP x, SEXP d)
448 {
449  int *dims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
450  m = dims[0], nret = (m < dims[1]) ? m : dims[1];
451  SEXP ret = PROTECT(duplicate(x)),
452  r_x = GET_SLOT(ret, Matrix_xSym);
453  double *dv = REAL(d), *rv = REAL(r_x);
454  int l_d = LENGTH(d); Rboolean d_full = (l_d == nret);
455  if (!d_full && l_d != 1)
456  error(_("diagonal to be added has wrong length"));
457 
458  if(d_full) for (int i = 0; i < nret; i++) rv[i * (m + 1)] += dv[i];
459  else for (int i = 0; i < nret; i++) rv[i * (m + 1)] += *dv;
460  UNPROTECT(1);
461  return ret;
462 }
463 
464 
465 
466 SEXP dgeMatrix_LU_(SEXP x, Rboolean warn_sing)
467 {
468  SEXP val = get_factors(x, "LU");
469  int *dims, npiv, info;
470 
471  if (val != R_NilValue) /* nothing to do if it's there in 'factors' slot */
472  return val;
473  dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
474  if (dims[0] < 1 || dims[1] < 1)
475  error(_("Cannot factor a matrix with zero extents"));
476  npiv = (dims[0] < dims[1]) ? dims[0] : dims[1];
477  val = PROTECT(NEW_OBJECT(MAKE_CLASS("denseLU")));
478  slot_dup(val, x, Matrix_xSym);
479  slot_dup(val, x, Matrix_DimSym);
480  slot_dup(val, x, Matrix_DimNamesSym);
481  F77_CALL(dgetrf)(dims, dims + 1, REAL(GET_SLOT(val, Matrix_xSym)),
482  dims,
483  INTEGER(ALLOC_SLOT(val, Matrix_permSym, INTSXP, npiv)),
484  &info);
485  if (info < 0)
486  error(_("Lapack routine %s returned error code %d"), "dgetrf", info);
487  else if (info > 0 && warn_sing)
488  warning(_("Exact singularity detected during LU decomposition: %s, i=%d."),
489  "U[i,i]=0", info);
490  UNPROTECT(1);
491  return set_factors(x, val, "LU");
492 }
493 // FIXME: also allow an interface to LAPACK's dgesvx() which uses LU fact.
494 // and then optionally does "equilibration" (row and column scaling)
495 // maybe also allow low-level interface to dgeEQU() ...
496 
497 SEXP dgeMatrix_LU(SEXP x, SEXP warn_singularity)
498 {
499  return dgeMatrix_LU_(x, asLogical(warn_singularity));
500 }
501 
502 SEXP dgeMatrix_determinant(SEXP x, SEXP logarithm)
503 {
504  int lg = asLogical(logarithm);
505  int *dims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
506  n = dims[0], sign = 1;
507  double modulus = lg ? 0. : 1; /* initialize; = result for n == 0 */
508 
509  if (n != dims[1])
510  error(_("Determinant requires a square matrix"));
511  if (n > 0) {
512  SEXP lu = dgeMatrix_LU_(x, /* do not warn about singular LU: */ FALSE);
513  int i, *jpvt = INTEGER(GET_SLOT(lu, Matrix_permSym));
514  double *luvals = REAL(GET_SLOT(lu, Matrix_xSym));
515 
516  for (i = 0; i < n; i++) if (jpvt[i] != (i + 1)) sign = -sign;
517  if (lg) {
518  for (i = 0; i < n; i++) {
519  double dii = luvals[i*(n + 1)]; /* ith diagonal element */
520  modulus += log(dii < 0 ? -dii : dii);
521  if (dii < 0) sign = -sign;
522  }
523  } else {
524  for (i = 0; i < n; i++)
525  modulus *= luvals[i*(n + 1)];
526  if (modulus < 0) {
527  modulus = -modulus;
528  sign = -sign;
529  }
530  }
531  }
532  return as_det_obj(modulus, lg, sign);
533 }
534 
535 SEXP dgeMatrix_solve(SEXP a)
536 {
537  /* compute the 1-norm of the matrix, which is needed
538  later for the computation of the reciprocal condition number. */
539  double aNorm = get_norm(a, "1");
540 
541  /* the LU decomposition : */
542  SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix"))),
543  lu = dgeMatrix_LU_(a, TRUE);
544  int *dims = INTEGER(GET_SLOT(lu, Matrix_DimSym)),
545  *pivot = INTEGER(GET_SLOT(lu, Matrix_permSym));
546 
547  /* prepare variables for the dgetri calls */
548  double *x, tmp;
549  int info, lwork = -1;
550 
551 
552  if (dims[0] != dims[1]) error(_("Solve requires a square matrix"));
553  slot_dup(val, lu, Matrix_xSym);
554  x = REAL(GET_SLOT(val, Matrix_xSym));
555  slot_dup(val, lu, Matrix_DimSym);
556 
557  if(dims[0]) /* the dimension is not zero */
558  {
559  /* is the matrix is *computationally* singular ? */
560  double rcond;
561  F77_CALL(dgecon)("1", dims, x, dims, &aNorm, &rcond,
562  (double *) R_alloc(4*dims[0], sizeof(double)),
563  (int *) R_alloc(dims[0], sizeof(int)), &info);
564  if (info)
565  error(_("error [%d] from Lapack 'dgecon()'"), info);
566  if(rcond < DOUBLE_EPS)
567  error(_("Lapack dgecon(): system computationally singular, reciprocal condition number = %g"),
568  rcond);
569 
570  /* only now try the inversion and check if the matrix is *exactly* singular: */
571  F77_CALL(dgetri)(dims, x, dims, pivot, &tmp, &lwork, &info);
572  lwork = (int) tmp;
573  F77_CALL(dgetri)(dims, x, dims, pivot,
574  (double *) R_alloc((size_t) lwork, sizeof(double)),
575  &lwork, &info);
576  if (info)
577  error(_("Lapack routine dgetri: system is exactly singular"));
578  }
579  UNPROTECT(1);
580  return val;
581 }
582 
583 SEXP dgeMatrix_matrix_solve(SEXP a, SEXP b)
584 {
585  SEXP val = PROTECT(dup_mMatrix_as_dgeMatrix(b)),
586  lu = PROTECT(dgeMatrix_LU_(a, TRUE));
587  int *adims = INTEGER(GET_SLOT(lu, Matrix_DimSym)),
588  *bdims = INTEGER(GET_SLOT(val, Matrix_DimSym));
589  int info, n = bdims[0], nrhs = bdims[1];
590 
591  if (adims[0] != n || adims[1] != n)
592  error(_("Dimensions of system to be solved are inconsistent"));
593  if(nrhs >= 1 && n >= 1) {
594  F77_CALL(dgetrs)("N", &n, &nrhs, REAL(GET_SLOT(lu, Matrix_xSym)), &n,
595  INTEGER(GET_SLOT(lu, Matrix_permSym)),
596  REAL(GET_SLOT(val, Matrix_xSym)), &n, &info);
597  if (info)
598  error(_("Lapack routine dgetrs: system is exactly singular"));
599  }
600  UNPROTECT(2);
601  return val;
602 }
603 
604 SEXP dgeMatrix_svd(SEXP x, SEXP nnu, SEXP nnv)
605 {
606  int /* nu = asInteger(nnu),
607  nv = asInteger(nnv), */
608  *dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
609  double *xx = REAL(GET_SLOT(x, Matrix_xSym));
610  SEXP val = PROTECT(allocVector(VECSXP, 3));
611 
612  if (dims[0] && dims[1]) {
613  int m = dims[0], n = dims[1], mm = (m < n)?m:n,
614  lwork = -1, info;
615  double tmp, *work;
616  int *iwork, n_iw = 8 * mm;
617  C_or_Alloca_TO(iwork, n_iw, int);
618 
619  SET_VECTOR_ELT(val, 0, allocVector(REALSXP, mm));
620  SET_VECTOR_ELT(val, 1, allocMatrix(REALSXP, m, mm));
621  SET_VECTOR_ELT(val, 2, allocMatrix(REALSXP, mm, n));
622  F77_CALL(dgesdd)("S", &m, &n, xx, &m,
623  REAL(VECTOR_ELT(val, 0)),
624  REAL(VECTOR_ELT(val, 1)), &m,
625  REAL(VECTOR_ELT(val, 2)), &mm,
626  &tmp, &lwork, iwork, &info);
627  lwork = (int) tmp;
628  C_or_Alloca_TO(work, lwork, double);
629 
630  F77_CALL(dgesdd)("S", &m, &n, xx, &m,
631  REAL(VECTOR_ELT(val, 0)),
632  REAL(VECTOR_ELT(val, 1)), &m,
633  REAL(VECTOR_ELT(val, 2)), &mm,
634  work, &lwork, iwork, &info);
635 
636  if(n_iw >= SMALL_4_Alloca) Free(iwork);
637  if(lwork >= SMALL_4_Alloca) Free(work);
638  }
639  UNPROTECT(1);
640  return val;
641 }
642 
643 const static double padec [] = /* for matrix exponential calculation. */
644 {
645  5.0000000000000000e-1,
646  1.1666666666666667e-1,
647  1.6666666666666667e-2,
648  1.6025641025641026e-3,
649  1.0683760683760684e-4,
650  4.8562548562548563e-6,
651  1.3875013875013875e-7,
652  1.9270852604185938e-9,
653 };
654 
662 SEXP dgeMatrix_exp(SEXP x)
663 {
664  const double one = 1.0, zero = 0.0;
665  const int i1 = 1;
666  int *Dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
667  const int n = Dims[1], nsqr = n * n, np1 = n + 1;
668 
669  SEXP val = PROTECT(duplicate(x));
670  int i, ilo, ilos, ihi, ihis, j, sqpow;
671  int *pivot = Calloc(n, int);
672  double *dpp = Calloc(nsqr, double), /* denominator power Pade' */
673  *npp = Calloc(nsqr, double), /* numerator power Pade' */
674  *perm = Calloc(n, double),
675  *scale = Calloc(n, double),
676  *v = REAL(GET_SLOT(val, Matrix_xSym)),
677  *work = Calloc(nsqr, double), inf_norm, m1_j/*= (-1)^j */, trshift;
678  R_CheckStack();
679 
680  if (n < 1 || Dims[0] != n)
681  error(_("Matrix exponential requires square, non-null matrix"));
682  if(n == 1) {
683  v[0] = exp(v[0]);
684  UNPROTECT(1);
685  return val;
686  }
687 
688  /* Preconditioning 1. Shift diagonal by average diagonal if positive. */
689  trshift = 0; /* determine average diagonal element */
690  for (i = 0; i < n; i++) trshift += v[i * np1];
691  trshift /= n;
692  if (trshift > 0.) { /* shift diagonal by -trshift */
693  for (i = 0; i < n; i++) v[i * np1] -= trshift;
694  }
695 
696  /* Preconditioning 2. Balancing with dgebal. */
697  F77_CALL(dgebal)("P", &n, v, &n, &ilo, &ihi, perm, &j);
698  if (j) error(_("dgeMatrix_exp: LAPACK routine dgebal returned %d"), j);
699  F77_CALL(dgebal)("S", &n, v, &n, &ilos, &ihis, scale, &j);
700  if (j) error(_("dgeMatrix_exp: LAPACK routine dgebal returned %d"), j);
701 
702  /* Preconditioning 3. Scaling according to infinity norm */
703  inf_norm = F77_CALL(dlange)("I", &n, &n, v, &n, work);
704  sqpow = (inf_norm > 0) ? (int) (1 + log(inf_norm)/log(2.)) : 0;
705  if (sqpow < 0) sqpow = 0;
706  if (sqpow > 0) {
707  double scale_factor = 1.0;
708  for (i = 0; i < sqpow; i++) scale_factor *= 2.;
709  for (i = 0; i < nsqr; i++) v[i] /= scale_factor;
710  }
711 
712  /* Pade' approximation. Powers v^8, v^7, ..., v^1 */
713  AZERO(npp, nsqr);
714  AZERO(dpp, nsqr);
715  m1_j = -1;
716  for (j = 7; j >=0; j--) {
717  double mult = padec[j];
718  /* npp = m * npp + padec[j] *m */
719  F77_CALL(dgemm)("N", "N", &n, &n, &n, &one, v, &n, npp, &n,
720  &zero, work, &n);
721  for (i = 0; i < nsqr; i++) npp[i] = work[i] + mult * v[i];
722  /* dpp = m * dpp + (m1_j * padec[j]) * m */
723  mult *= m1_j;
724  F77_CALL(dgemm)("N", "N", &n, &n, &n, &one, v, &n, dpp, &n,
725  &zero, work, &n);
726  for (i = 0; i < nsqr; i++) dpp[i] = work[i] + mult * v[i];
727  m1_j *= -1;
728  }
729  /* Zero power */
730  for (i = 0; i < nsqr; i++) dpp[i] *= -1.;
731  for (j = 0; j < n; j++) {
732  npp[j * np1] += 1.;
733  dpp[j * np1] += 1.;
734  }
735 
736  /* Pade' approximation is solve(dpp, npp) */
737  F77_CALL(dgetrf)(&n, &n, dpp, &n, pivot, &j);
738  if (j) error(_("dgeMatrix_exp: dgetrf returned error code %d"), j);
739  F77_CALL(dgetrs)("N", &n, &n, dpp, &n, pivot, npp, &n, &j);
740  if (j) error(_("dgeMatrix_exp: dgetrs returned error code %d"), j);
741  Memcpy(v, npp, nsqr);
742 
743  /* Now undo all of the preconditioning */
744  /* Preconditioning 3: square the result for every power of 2 */
745  while (sqpow--) {
746  F77_CALL(dgemm)("N", "N", &n, &n, &n, &one, v, &n, v, &n,
747  &zero, work, &n);
748  Memcpy(v, work, nsqr);
749  }
750  /* Preconditioning 2: apply inverse scaling */
751  for (j = 0; j < n; j++)
752  for (i = 0; i < n; i++)
753  v[i + j * n] *= scale[i]/scale[j];
754 
755 
756  /* 2 b) Inverse permutation (if not the identity permutation) */
757  if (ilo != 1 || ihi != n) {
758  /* Martin Maechler's code */
759 
760 #define SWAP_ROW(I,J) F77_CALL(dswap)(&n, &v[(I)], &n, &v[(J)], &n)
761 
762 #define SWAP_COL(I,J) F77_CALL(dswap)(&n, &v[(I)*n], &i1, &v[(J)*n], &i1)
763 
764 #define RE_PERMUTE(I) \
765  int p_I = (int) (perm[I]) - 1; \
766  SWAP_COL(I, p_I); \
767  SWAP_ROW(I, p_I)
768 
769  /* reversion of "leading permutations" : in reverse order */
770  for (i = (ilo - 1) - 1; i >= 0; i--) {
771  RE_PERMUTE(i);
772  }
773 
774  /* reversion of "trailing permutations" : applied in forward order */
775  for (i = (ihi + 1) - 1; i < n; i++) {
776  RE_PERMUTE(i);
777  }
778  }
779 
780  /* Preconditioning 1: Trace normalization */
781  if (trshift > 0.) {
782  double mult = exp(trshift);
783  for (i = 0; i < nsqr; i++) v[i] *= mult;
784  }
785 
786  /* Clean up */
787  Free(work); Free(scale); Free(perm); Free(npp); Free(dpp); Free(pivot);
788  UNPROTECT(1);
789  return val;
790 }
791 
792 SEXP dgeMatrix_Schur(SEXP x, SEXP vectors, SEXP isDGE)
793 {
794 // 'x' is either a traditional matrix or a dgeMatrix, as indicated by isDGE.
795  int *dims, n, vecs = asLogical(vectors), is_dge = asLogical(isDGE),
796  info, izero = 0, lwork = -1, nprot = 1;
797 
798  if(is_dge) {
799  dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
800  } else { // traditional matrix
801  dims = INTEGER(getAttrib(x, R_DimSymbol));
802  if(!isReal(x)) { // may not be "numeric" ..
803  x = PROTECT(coerceVector(x, REALSXP)); // -> maybe error
804  nprot++;
805  }
806  }
807  double *work, tmp;
808  const char *nms[] = {"WR", "WI", "T", "Z", ""};
809  SEXP val = PROTECT(Rf_mkNamed(VECSXP, nms));
810 
811  n = dims[0];
812  if (n != dims[1] || n < 1)
813  error(_("dgeMatrix_Schur: argument x must be a non-null square matrix"));
814  SET_VECTOR_ELT(val, 0, allocVector(REALSXP, n));
815  SET_VECTOR_ELT(val, 1, allocVector(REALSXP, n));
816  SET_VECTOR_ELT(val, 2, allocMatrix(REALSXP, n, n));
817  Memcpy(REAL(VECTOR_ELT(val, 2)),
818  REAL(is_dge ? GET_SLOT(x, Matrix_xSym) : x),
819  n * n);
820  SET_VECTOR_ELT(val, 3, allocMatrix(REALSXP, vecs ? n : 0, vecs ? n : 0));
821  F77_CALL(dgees)(vecs ? "V" : "N", "N", NULL, dims, (double *) NULL, dims, &izero,
822  (double *) NULL, (double *) NULL, (double *) NULL, dims,
823  &tmp, &lwork, (int *) NULL, &info);
824  if (info) error(_("dgeMatrix_Schur: first call to dgees failed"));
825  lwork = (int) tmp;
826  C_or_Alloca_TO(work, lwork, double);
827 
828  F77_CALL(dgees)(vecs ? "V" : "N", "N", NULL, dims, REAL(VECTOR_ELT(val, 2)), dims,
829  &izero, REAL(VECTOR_ELT(val, 0)), REAL(VECTOR_ELT(val, 1)),
830  REAL(VECTOR_ELT(val, 3)), dims, work, &lwork,
831  (int *) NULL, &info);
832  if(lwork >= SMALL_4_Alloca) Free(work);
833  if (info) error(_("dgeMatrix_Schur: dgees returned code %d"), info);
834  UNPROTECT(nprot);
835  return val;
836 } // dgeMatrix_Schur
837 
838 // colSums(), colMeans(), rowSums() and rowMeans() -- called from ../R/colSums.R
839 SEXP dgeMatrix_colsums(SEXP x, SEXP naRmP, SEXP cols, SEXP mean)
840 {
841  int keepNA = !asLogical(naRmP); // <==> na.rm = FALSE, the default
842  int doMean = asLogical(mean);
843  int useCols = asLogical(cols);
844  int *dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
845  int i, j, m = dims[0], n = dims[1];
846  SEXP ans = PROTECT(allocVector(REALSXP, (useCols) ? n : m));
847  double *aa = REAL(ans), *xx = REAL(GET_SLOT(x, Matrix_xSym));
848 
849  if (useCols) { /* col(Sums|Means) : */
850  int cnt = m; // := number of 'valid' entries in current column
851  for (j = 0; j < n; j++) { // column j :
852  double *x_j = xx + m * j, s = 0.;
853 
854  if (keepNA)
855  for (i = 0; i < m; i++) s += x_j[i];
856  else {
857  cnt = 0;
858  for (i = 0; i < m; i++)
859  if (!ISNAN(x_j[i])) {cnt++; s += x_j[i];}
860  }
861  if (doMean) {
862  if (cnt > 0) s /= cnt; else s = NA_REAL;
863  }
864  aa[j] = s;
865  }
866  } else { /* row(Sums|Means) : */
867  Rboolean do_count = (!keepNA) && doMean;
868  int *cnt = (int*) NULL;
869  if(do_count) { C_or_Alloca_TO(cnt, m, int); }
870 
871  // (taking care to access x contiguously: vary i inside j)
872  for (i = 0; i < m; i++) {
873  aa[i] = 0.;
874  if(do_count) cnt[i] = 0;
875  }
876  for (j = 0; j < n; j++) {
877  if (keepNA)
878  for (i = 0; i < m; i++) aa[i] += xx[i + j * m];
879  else
880  for (i = 0; i < m; i++) {
881  double el = xx[i + j * m];
882  if (!ISNAN(el)) {
883  aa[i] += el;
884  if (doMean) cnt[i]++;
885  }
886  }
887  }
888  if (doMean) {
889  if (keepNA)
890  for (i = 0; i < m; i++) aa[i] /= n;
891  else
892  for (i = 0; i < m; i++)
893  aa[i] = (cnt[i] > 0) ? aa[i]/cnt[i] : NA_REAL;
894  }
895  if(do_count && m >= SMALL_4_Alloca) Free(cnt);
896  }
897 
898  SEXP nms = VECTOR_ELT(GET_SLOT(x, Matrix_DimNamesSym), useCols ? 1 : 0);
899  if(!isNull(nms))
900  setAttrib(ans, R_NamesSymbol, duplicate(nms));
901  UNPROTECT(1);
902  return ans;
903 }
SEXP Matrix_DimSym
Definition: Syms.h:2
#define class_P(_x_)
Definition: Mutils.h:178
SEXP dgeMatrix_LU_(SEXP x, Rboolean warn_sing)
Definition: dgeMatrix.c:466
#define C_or_Alloca_TO(_VAR_, _N_, _TYPE_)
Definition: Mutils.h:50
#define DGE_MAT_CROSS_1
SEXP as_det_obj(double val, int log, int sign)
Definition: Mutils.c:89
SEXP lgeMatrix_setDiag(SEXP x, SEXP d)
Definition: dgeMatrix.c:436
SEXP dgeMatrix_colsums(SEXP x, SEXP naRmP, SEXP cols, SEXP mean)
Definition: dgeMatrix.c:839
SEXP Matrix_xSym
Definition: Syms.h:2
#define slot_dup(dest, src, sym)
Definition: Mutils.h:149
SEXP Matrix_factorSym
Definition: Syms.h:2
SEXP geMatrix_geMatrix_crossprod(SEXP x, SEXP y, SEXP trans)
Definition: dgeMatrix.c:202
#define RE_PERMUTE(I)
SEXP Matrix_DimNamesSym
Definition: Syms.h:2
#define geMatrix_setDiag_1
#define DGE_DGE_CROSS_DO(_X_X_, _Y_Y_)
SEXP _geMatrix_matrix_crossprod(SEXP x, SEXP y, SEXP trans)
as dgeMatrix_matrix_crossprod() but x can be [dln]geMatrix
Definition: dgeMatrix.c:283
static R_INLINE SEXP mMatrix_as_dgeMatrix(SEXP A)
Definition: Mutils.h:334
SEXP dgeMatrix_matrix_mm(SEXP a, SEXP bP, SEXP right)
Definition: dgeMatrix.c:299
SEXP dup_mMatrix_as_dgeMatrix(SEXP A)
Definition: Mutils.c:852
SEXP geMatrix_crossprod(SEXP x, SEXP trans)
Definition: dgeMatrix.c:137
SEXP dim_validate(SEXP Dim, const char *name)
Definition: Mutils.c:318
SEXP dgeMatrix_svd(SEXP x, SEXP nnu, SEXP nnv)
Definition: dgeMatrix.c:604
#define geMatrix_getDiag_1
void F77_NAME() dgesdd(const char *jobz, const int *m, const int *n, double *a, const int *lda, double *s, double *u, const int *ldu, double *vt, const int *ldvt, double *work, const int *lwork, int *iwork, int *info)
SEXP dgeMatrix_norm(SEXP obj, SEXP type)
Definition: dgeMatrix.c:49
SEXP dgeMatrix_solve(SEXP a)
Definition: dgeMatrix.c:535
char La_rcond_type(const char *typstr)
Definition: Mutils.c:27
#define SMALL_4_Alloca
Definition: Mutils.h:47
static double get_norm(SEXP obj, const char *typstr)
Definition: dgeMatrix.c:30
SEXP lgeMatrix_getDiag(SEXP x)
Definition: dgeMatrix.c:396
SEXP dgeMatrix_validate(SEXP obj)
Definition: dgeMatrix.c:16
SEXP set_factors(SEXP obj, SEXP val, char *nm)
Caches 'val' in the 'factors' slot of obj, i.e.
Definition: Mutils.c:130
#define DGE_CROSS_1
SEXP get_factors(SEXP obj, char *nm)
Definition: Mutils.c:106
#define DGE_MAT_CROSS_DO(_X_X_)
#define geMatrix_getDiag_2
SEXP _geMatrix__geMatrix_crossprod(SEXP x, SEXP y, SEXP trans)
As dgeMatrix_dgeMatrix_crossprod(), but x and y can be [dln]geMatrix.
Definition: dgeMatrix.c:190
#define DGE_DGE_CROSS_1
#define AZERO(x, n)
Definition: Mutils.h:140
SEXP dgeMatrix_matrix_crossprod(SEXP x, SEXP y, SEXP trans)
Definition: dgeMatrix.c:211
SEXP Matrix_permSym
Definition: Syms.h:2
#define geMatrix_setDiag_2
#define _(String)
Definition: Mutils.h:32
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
SEXP dgeMatrix_crossprod(SEXP x, SEXP trans)
Definition: dgeMatrix.c:76
static const double padec[]
Definition: dgeMatrix.c:643
SEXP dgeMatrix_exp(SEXP x)
Matrix exponential - based on the corrected code for Octave's expm function.
Definition: dgeMatrix.c:662
char La_norm_type(const char *typstr)
Definition: Mutils.c:8
#define DGE_CROSS_DO(_X_X_)
double * gematrix_real_x(SEXP x, int nn)
Definition: dgeMatrix.c:109
#define DGE_MAT_MM_1(N_PROT)
SEXP _geMatrix_crossprod(SEXP x, SEXP trans)
As dgeMatrix_crossprod(), but x can be [dln]geMatrix.
Definition: dgeMatrix.c:130
SEXP dgeMatrix_matrix_solve(SEXP a, SEXP b)
Definition: dgeMatrix.c:583
SEXP geMatrix_matrix_crossprod(SEXP x, SEXP y, SEXP trans)
Definition: dgeMatrix.c:291
SEXP dgeMatrix_determinant(SEXP x, SEXP logarithm)
Definition: dgeMatrix.c:502
SEXP dgeMatrix_dgeMatrix_crossprod(SEXP x, SEXP y, SEXP trans)
Definition: dgeMatrix.c:145
SEXP dgeMatrix_Schur(SEXP x, SEXP vectors, SEXP isDGE)
Definition: dgeMatrix.c:792
#define DGE_MAT_MM_DO(_A_X_, _B_X_)
SEXP dgeMatrix_setDiag(SEXP x, SEXP d)
Definition: dgeMatrix.c:411
SEXP geMatrix_matrix_mm(SEXP a, SEXP b, SEXP right)
%*% – generalized from dge to *ge():
Definition: dgeMatrix.c:363
SEXP _geMatrix_matrix_mm(SEXP a, SEXP b, SEXP right)
as dgeMatrix_matrix_mm() but a can be [dln]geMatrix
Definition: dgeMatrix.c:355
static R_INLINE Rboolean any_NA_in_x(SEXP obj)
Check if slot(obj, "x") contains any NA (or NaN).
Definition: Mutils.h:277
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 dMatrix_validate(SEXP obj)
Definition: dgeMatrix.c:4
SEXP dgeMatrix_LU(SEXP x, SEXP warn_singularity)
Definition: dgeMatrix.c:497
SEXP dgeMatrix_addDiag(SEXP x, SEXP d)
Definition: dgeMatrix.c:447
SEXP dense_nonpacked_validate(SEXP obj)
Definition: Mutils.c:310
SEXP dgeMatrix_getDiag(SEXP x)
Definition: dgeMatrix.c:374
SEXP dgeMatrix_rcond(SEXP obj, SEXP type)
Definition: dgeMatrix.c:54