Matrix  $Rev: 3071 $ at $LastChangedDate: 2015-03-26 15:35:47 +0100 (Thu, 26 Mar 2015) $
factorizations.c
Go to the documentation of this file.
1 #include "factorizations.h"
2 
4 {
5  SEXP val;
6  if (isString(val = dim_validate(GET_SLOT(obj, Matrix_DimSym),
7  "MatrixFactorization")))
8  return(val);
9  return ScalarLogical(1);
10 }
11 
12 SEXP LU_validate(SEXP obj)
13 {
14  SEXP x = GET_SLOT(obj, Matrix_xSym),
15  Dim = GET_SLOT(obj, Matrix_DimSym);
16  int m = INTEGER(Dim)[0], n = INTEGER(Dim)[1]; // checked already in MatrixF.._validate()
17  if(TYPEOF(x) != REALSXP)
18  return mkString(_("x slot is not \"double\""));
19  if(XLENGTH(x) != ((double) m) * n)
20  return mkString(_("x slot is not of correct length"));
21  return dimNames_validate(obj);
22 }
23 
24 SEXP BunchKaufman_validate(SEXP obj)
25 {
26  // TODO
27  return ScalarLogical(1);
28 }
29 
30 SEXP pBunchKaufman_validate(SEXP obj)
31 {
32  // TODO
33  return ScalarLogical(1);
34 }
35 
36 SEXP Cholesky_validate(SEXP obj)
37 {
38  // TODO
39  return ScalarLogical(1);
40 }
41 
42 SEXP pCholesky_validate(SEXP obj)
43 {
44  // TODO
45  return ScalarLogical(1);
46 }
47 
48 #ifdef _Matrix_has_SVD_
49 SEXP SVD_validate(SEXP obj)
50 {
51  return ScalarLogical(1);
52 }
53 #endif
54 
55 SEXP LU_expand(SEXP x)
56 {
57  const char *nms[] = {"L", "U", "P", ""};
58  // x[,] is m x n (using LAPACK dgetrf notation)
59  SEXP L, U, P, val = PROTECT(Rf_mkNamed(VECSXP, nms)),
60  lux = GET_SLOT(x, Matrix_xSym),
61  dd = GET_SLOT(x, Matrix_DimSym);
62  int *iperm, *perm, *pivot = INTEGER(GET_SLOT(x, Matrix_permSym)),
63  *dim = INTEGER(dd), m = dim[0], n = dim[1], nn = m, i;
64  size_t m_ = (size_t) m; // to prevent integer (multiplication) overflow
65  Rboolean is_sq = (n == m), L_is_tri = TRUE, U_is_tri = TRUE;
66 
67  // nn := min(n,m) == length(pivot[])
68  if(!is_sq) {
69  if(n < m) { // "long"
70  nn = n;
71  L_is_tri = FALSE;
72  } else { // m < n : "wide"
73  U_is_tri = FALSE;
74  }
75  }
76 
77  SET_VECTOR_ELT(val, 0, NEW_OBJECT(MAKE_CLASS(L_is_tri ? "dtrMatrix":"dgeMatrix")));
78  SET_VECTOR_ELT(val, 1, NEW_OBJECT(MAKE_CLASS(U_is_tri ? "dtrMatrix":"dgeMatrix")));
79  SET_VECTOR_ELT(val, 2, NEW_OBJECT(MAKE_CLASS("pMatrix")));
80  L = VECTOR_ELT(val, 0);
81  U = VECTOR_ELT(val, 1);
82  P = VECTOR_ELT(val, 2);
83  if(is_sq || !L_is_tri) {
84  SET_SLOT(L, Matrix_xSym, duplicate(lux));
85  SET_SLOT(L, Matrix_DimSym, duplicate(dd));
86  } else { // !is_sq && L_is_tri -- m < n -- "wide" -- L is m x m
87  size_t m2 = m_ * m;
88  double *Lx = REAL(ALLOC_SLOT(L, Matrix_xSym, REALSXP, m2));
89  int *dL = INTEGER(ALLOC_SLOT(L, Matrix_DimSym, INTSXP, 2));
90  dL[0] = dL[1] = m;
91  // fill lower-diagonal (non-{0,1}) part -- remainder by make_d_matrix*() below:
92  Memcpy(Lx, REAL(lux), m2);
93  }
94  if(is_sq || !U_is_tri) {
95  SET_SLOT(U, Matrix_xSym, duplicate(lux));
96  SET_SLOT(U, Matrix_DimSym, duplicate(dd));
97  } else { // !is_sq && U_is_tri -- m > n -- "long" -- U is n x n
98  double *Ux = REAL(ALLOC_SLOT(U, Matrix_xSym, REALSXP, ((size_t) n) * n)),
99  *xx = REAL(lux);
100  int *dU = INTEGER(ALLOC_SLOT(U, Matrix_DimSym, INTSXP, 2));
101  dU[0] = dU[1] = n;
102  /* fill upper-diagonal (non-0) part -- remainder by make_d_matrix*() below:
103  * this is more complicated than in the L case, as the x / lux part we need
104  * is *not* continguous: Memcpy(Ux, REAL(lux), n * n); -- is WRONG */
105  for (size_t j = 0; j < n; j++) {
106  Memcpy(Ux+j*n, xx+j*m, j+1);
107  // for (i = 0; i <= j; i++)
108  // Ux[i + j*n] = xx[i + j*m];
109  }
110  }
111  if(L_is_tri) {
112  SET_SLOT(L, Matrix_uploSym, mkString("L"));
113  SET_SLOT(L, Matrix_diagSym, mkString("U"));
114  make_d_matrix_triangular(REAL(GET_SLOT(L, Matrix_xSym)), L);
115  } else { // L is "unit-diagonal" trapezoidal -- m > n -- "long"
116  // fill the upper right part with 0 *and* the diagonal with 1
117  double *Lx = REAL(GET_SLOT(L, Matrix_xSym));
118  size_t ii;
119  for (i = 0, ii = 0; i < n; i++, ii+=(m+1)) { // ii = i*(m+1)
120  Lx[ii] = 1.;
121  for (size_t j = i*m_; j < ii; j++)
122  Lx[j] = 0.;
123  }
124  }
125 
126  if(U_is_tri) {
127  SET_SLOT(U, Matrix_uploSym, mkString("U"));
128  SET_SLOT(U, Matrix_diagSym, mkString("N"));
129  make_d_matrix_triangular(REAL(GET_SLOT(U, Matrix_xSym)), U);
130  } else { // U is trapezoidal -- m < n
131  // fill the lower left part with 0
132  double *Ux = REAL(GET_SLOT(U, Matrix_xSym));
133  for (i = 0; i < m; i++)
134  for (size_t j = i*(m_+1) +1; j < (i+1)*m_; j++)
135  Ux[j] = 0.;
136  }
137 
138  SET_SLOT(P, Matrix_DimSym, duplicate(dd));
139  if(!is_sq) // m != n -- P is m x m
140  INTEGER(GET_SLOT(P, Matrix_DimSym))[1] = m;
141  perm = INTEGER(ALLOC_SLOT(P, Matrix_permSym, INTSXP, m));
142  C_or_Alloca_TO(iperm, m, int);
143 
144  for (i = 0; i < m; i++) iperm[i] = i + 1; /* initialize permutation*/
145  for (i = 0; i < nn; i++) { /* generate inverse permutation */
146  int newp = pivot[i] - 1;
147  if (newp != i) { // swap
148  int tmp = iperm[i]; iperm[i] = iperm[newp]; iperm[newp] = tmp;
149  }
150  }
151  // invert the inverse
152  for (i = 0; i < m; i++) perm[iperm[i] - 1] = i + 1;
153 
154  if(m >= SMALL_4_Alloca) Free(iperm);
155  UNPROTECT(1);
156  return val;
157 }
SEXP Matrix_DimSym
Definition: Syms.h:2
#define C_or_Alloca_TO(_VAR_, _N_, _TYPE_)
Definition: Mutils.h:50
SEXP Matrix_xSym
Definition: Syms.h:2
#define XLENGTH(x)
Definition: Mutils.h:39
SEXP Matrix_uploSym
Definition: Syms.h:2
SEXP Cholesky_validate(SEXP obj)
SEXP dim_validate(SEXP Dim, const char *name)
Definition: Mutils.c:318
SEXP pBunchKaufman_validate(SEXP obj)
void make_d_matrix_triangular(double *to, SEXP from) void make_i_matrix_triangular(int *to
#define SMALL_4_Alloca
Definition: Mutils.h:47
SEXP MatrixFactorization_validate(SEXP obj)
Definition: factorizations.c:3
SEXP Matrix_permSym
Definition: Syms.h:2
#define _(String)
Definition: Mutils.h:32
SEXP LU_expand(SEXP x)
SEXP BunchKaufman_validate(SEXP obj)
SEXP dimNames_validate(SEXP obj)
Definition: Mutils.c:343
SEXP Matrix_diagSym
Definition: Syms.h:2
SEXP pCholesky_validate(SEXP obj)
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 LU_validate(SEXP obj)