Matrix  \$Rev: 3071 \$ at \$LastChangedDate: 2015-03-26 15:35:47 +0100 (Thu, 26 Mar 2015) \$
dense.c
Go to the documentation of this file.
1 #include "dense.h"
2 #include "Mutils.h"
3 #include "chm_common.h"
4
21 static
22 int left_cyclic(double x[], int ldx, int j, int k,
23  double cosines[], double sines[])
24 {
25  double *lastcol;
26  int i, jj;
27
28  if (j >= k)
29  error(_("incorrect left cyclic shift, j (%d) >= k (%d)"), j, k);
30  if (j < 0)
31  error(_("incorrect left cyclic shift, j (%d) < 0"), j, k);
32  if (ldx < k)
33  error(_("incorrect left cyclic shift, k (%d) > ldx (%d)"), k, ldx);
34  lastcol = (double*) R_alloc(k+1, sizeof(double));
35  /* keep a copy of column j */
36  for(i = 0; i <= j; i++) lastcol[i] = x[i + j*ldx];
37  /* For safety, zero the rest */
38  for(i = j+1; i <= k; i++) lastcol[i] = 0.;
39  for(jj = j+1; jj <= k; jj++) { /* columns to be shifted */
40  int diagind = jj*(ldx+1), ind = (jj-j) - 1;
41  double tmp = x[diagind], cc, ss;
42  /* Calculate the Givens rotation. */
43  /* This modified the super-diagonal element */
44  F77_CALL(drotg)(x + diagind-1, &tmp, cosines + ind, sines + ind);
45  cc = cosines[ind]; ss = sines[ind];
46  /* Copy column jj+1 to column jj. */
47  for(i = 0; i < jj; i++) x[i + (jj-1)*ldx] = x[i+jj*ldx];
48  /* Apply rotation to columns up to k */
49  for(i = jj; i < k; i++) {
50  tmp = cc*x[(jj-1)+i*ldx] + ss*x[jj+i*ldx];
51  x[jj+i*ldx] = cc*x[jj+i*ldx] - ss*x[(jj-1)+i*ldx];
52  x[(jj-1)+i*ldx] = tmp;
53  }
54  /* Apply rotation to lastcol */
55  lastcol[jj] = -ss*lastcol[jj-1]; lastcol[jj-1] *= cc;
56  }
57  /* Copy lastcol to column k */
58  for(i = 0; i <= k; i++) x[i+k*ldx] = lastcol[i];
59  return 0;
60 }
61
62 static
63 SEXP getGivens(double x[], int ldx, int jmin, int rank)
64 {
65  int shiftlen = (rank - jmin) - 1;
66  SEXP ans = PROTECT(allocVector(VECSXP, 4)), nms, cosines, sines;
67
68  SET_VECTOR_ELT(ans, 0, ScalarInteger(jmin));
69  SET_VECTOR_ELT(ans, 1, ScalarInteger(rank));
70  SET_VECTOR_ELT(ans, 2, cosines = allocVector(REALSXP, shiftlen));
71  SET_VECTOR_ELT(ans, 3, sines = allocVector(REALSXP, shiftlen));
72  setAttrib(ans, R_NamesSymbol, nms = allocVector(STRSXP, 4));
73  SET_STRING_ELT(nms, 0, mkChar("jmin"));
74  SET_STRING_ELT(nms, 1, mkChar("rank"));
75  SET_STRING_ELT(nms, 2, mkChar("cosines"));
76  SET_STRING_ELT(nms, 3, mkChar("sines"));
77  if (left_cyclic(x, ldx, jmin, rank - 1, REAL(cosines), REAL(sines)))
78  error(_("Unknown error in getGivens"));
79  UNPROTECT(1);
80  return ans;
81 }
82
83 SEXP checkGivens(SEXP X, SEXP jmin, SEXP rank)
84 {
85  SEXP ans = PROTECT(allocVector(VECSXP, 2)),
86  Xcp = PROTECT(duplicate(X));
87  int *Xdims;
88
89  if (!(isReal(X) & isMatrix(X)))
90  error(_("X must be a numeric (double precision) matrix"));
91  Xdims = INTEGER(coerceVector(getAttrib(X, R_DimSymbol), INTSXP));
92  SET_VECTOR_ELT(ans, 1, getGivens(REAL(Xcp), Xdims[0],
93  asInteger(jmin), asInteger(rank)));
94  SET_VECTOR_ELT(ans, 0, Xcp);
95  UNPROTECT(2);
96  return ans;
97 }
98
99 SEXP lsq_dense_Chol(SEXP X, SEXP y)
100 {
101  SEXP ans;
102  int info, n, p, k, *Xdims, *ydims;
103  double *xpx, d_one = 1., d_zero = 0.;
104
105  if (!(isReal(X) & isMatrix(X)))
106  error(_("X must be a numeric (double precision) matrix"));
107  Xdims = INTEGER(coerceVector(getAttrib(X, R_DimSymbol), INTSXP));
108  n = Xdims[0];
109  p = Xdims[1];
110  if (!(isReal(y) & isMatrix(y)))
111  error(_("y must be a numeric (double precision) matrix"));
112  ydims = INTEGER(coerceVector(getAttrib(y, R_DimSymbol), INTSXP));
113  if (ydims[0] != n)
114  error(_(
115  "number of rows in y (%d) does not match number of rows in X (%d)"),
116  ydims[0], n);
117  k = ydims[1];
118  if (k < 1 || p < 1) return allocMatrix(REALSXP, p, k);
119  ans = PROTECT(allocMatrix(REALSXP, p, k));
120  F77_CALL(dgemm)("T", "N", &p, &k, &n, &d_one, REAL(X), &n, REAL(y), &n,
121  &d_zero, REAL(ans), &p);
122  xpx = (double *) R_alloc(p * p, sizeof(double));
123  F77_CALL(dsyrk)("U", "T", &p, &n, &d_one, REAL(X), &n, &d_zero,
124  xpx, &p);
125  F77_CALL(dposv)("U", &p, &k, xpx, &p, REAL(ans), &p, &info);
126  if (info) error(_("Lapack routine dposv returned error code %d"), info);
127  UNPROTECT(1);
128  return ans;
129 }
130
131
132 SEXP lsq_dense_QR(SEXP X, SEXP y)
133 {
134  SEXP ans;
135  int info, n, p, k, *Xdims, *ydims, lwork;
136  double *work, tmp, *xvals;
137
138  if (!(isReal(X) & isMatrix(X)))
139  error(_("X must be a numeric (double precision) matrix"));
140  Xdims = INTEGER(coerceVector(getAttrib(X, R_DimSymbol), INTSXP));
141  n = Xdims[0];
142  p = Xdims[1];
143  if (!(isReal(y) & isMatrix(y)))
144  error(_("y must be a numeric (double precision) matrix"));
145  ydims = INTEGER(coerceVector(getAttrib(y, R_DimSymbol), INTSXP));
146  if (ydims[0] != n)
147  error(_(
148  "number of rows in y (%d) does not match number of rows in X (%d)"),
149  ydims[0], n);
150  k = ydims[1];
151  if (k < 1 || p < 1) return allocMatrix(REALSXP, p, k);
152  xvals = (double *) Memcpy(R_alloc(n * p, sizeof(double)), REAL(X), n * p);
153  ans = PROTECT(duplicate(y));
154  lwork = -1;
155  F77_CALL(dgels)("N", &n, &p, &k, xvals, &n, REAL(ans), &n,
156  &tmp, &lwork, &info);
157  if (info)
158  error(_("First call to Lapack routine dgels returned error code %d"),
159  info);
160  lwork = (int) tmp;
161  work = (double *) R_alloc(lwork, sizeof(double));
162  F77_CALL(dgels)("N", &n, &p, &k, xvals, &n, REAL(ans), &n,
163  work, &lwork, &info);
164  if (info)
165  error(_("Second call to Lapack routine dgels returned error code %d"),
166  info);
167  UNPROTECT(1);
168  return ans;
169 }
170
171 SEXP lapack_qr(SEXP Xin, SEXP tl)
172 {
173  SEXP ans, Givens, Gcpy, nms, pivot, qraux, X;
174  int i, n, nGivens = 0, p, trsz, *Xdims, rank;
175  double rcond = 0., tol = asReal(tl), *work;
176
177  if (!(isReal(Xin) & isMatrix(Xin)))
178  error(_("X must be a real (numeric) matrix"));
179  if (tol < 0.) error(_("tol, given as %g, must be non-negative"), tol);
180  if (tol > 1.) error(_("tol, given as %g, must be <= 1"), tol);
181  ans = PROTECT(allocVector(VECSXP,5));
182  SET_VECTOR_ELT(ans, 0, X = duplicate(Xin));
183  Xdims = INTEGER(coerceVector(getAttrib(X, R_DimSymbol), INTSXP));
184  n = Xdims[0]; p = Xdims[1];
185  SET_VECTOR_ELT(ans, 2, qraux = allocVector(REALSXP, (n < p) ? n : p));
186  SET_VECTOR_ELT(ans, 3, pivot = allocVector(INTSXP, p));
187  for (i = 0; i < p; i++) INTEGER(pivot)[i] = i + 1;
188  trsz = (n < p) ? n : p; /* size of triangular part of decomposition */
189  rank = trsz;
190  Givens = PROTECT(allocVector(VECSXP, rank - 1));
191  setAttrib(ans, R_NamesSymbol, nms = allocVector(STRSXP, 5));
192  SET_STRING_ELT(nms, 0, mkChar("qr"));
193  SET_STRING_ELT(nms, 1, mkChar("rank"));
194  SET_STRING_ELT(nms, 2, mkChar("qraux"));
195  SET_STRING_ELT(nms, 3, mkChar("pivot"));
196  SET_STRING_ELT(nms, 4, mkChar("Givens"));
197  if (n > 0 && p > 0) {
198  int info, *iwork, lwork;
199  double *xpt = REAL(X), tmp;
200
201  lwork = -1;
202  F77_CALL(dgeqrf)(&n, &p, xpt, &n, REAL(qraux), &tmp, &lwork, &info);
203  if (info)
204  error(_("First call to dgeqrf returned error code %d"), info);
205  lwork = (int) tmp;
206  work = (double *) R_alloc((lwork < 3*trsz) ? 3*trsz : lwork,
207  sizeof(double));
208  F77_CALL(dgeqrf)(&n, &p, xpt, &n, REAL(qraux), work, &lwork, &info);
209  if (info)
210  error(_("Second call to dgeqrf returned error code %d"), info);
211  iwork = (int *) R_alloc(trsz, sizeof(int));
212  F77_CALL(dtrcon)("1", "U", "N", &rank, xpt, &n, &rcond,
213  work, iwork, &info);
214  if (info)
215  error(_("Lapack routine dtrcon returned error code %d"), info);
216  while (rcond < tol) { /* check diagonal elements */
217  double minabs = (xpt[0] < 0.) ? -xpt[0]: xpt[0];
218  int jmin = 0;
219  for (i = 1; i < rank; i++) {
220  double el = xpt[i*(n+1)];
221  el = (el < 0.) ? -el: el;
222  if (el < minabs) {
223  jmin = i;
224  minabs = el;
225  }
226  }
227  if (jmin < (rank - 1)) {
228  SET_VECTOR_ELT(Givens, nGivens, getGivens(xpt, n, jmin, rank));
229  nGivens++;
230  }
231  rank--;
232  F77_CALL(dtrcon)("1", "U", "N", &rank, xpt, &n, &rcond,
233  work, iwork, &info);
234  if (info)
235  error(_("Lapack routine dtrcon returned error code %d"), info);
236  }
237  }
238  SET_VECTOR_ELT(ans, 4, Gcpy = allocVector(VECSXP, nGivens));
239  for (i = 0; i < nGivens; i++)
240  SET_VECTOR_ELT(Gcpy, i, VECTOR_ELT(Givens, i));
241  SET_VECTOR_ELT(ans, 1, ScalarInteger(rank));
242  setAttrib(ans, install("useLAPACK"), ScalarLogical(1));
243  setAttrib(ans, install("rcond"), ScalarReal(rcond));
244  UNPROTECT(2);
245  return ans;
246 }
247
248 SEXP dense_to_Csparse(SEXP x)
249 {
250  CHM_DN chxd = AS_CHM_xDN(PROTECT(mMatrix_as_geMatrix(x)));
251  /* cholmod_dense_to_sparse() in CHOLMOD/Core/ below does only work for
252  "REAL" 'xtypes', i.e. *not* for "nMatrix".
253  ===> need "_x" in above AS_CHM_xDN() call.
254
255  Also it cannot keep symmetric / triangular, hence the
256  as_geMatrix() above. Note that this is already a *waste* for
257  symmetric matrices; However, we could conceivably use an
258  enhanced cholmod_dense_to_sparse(), with an extra boolean
259  argument for symmetry.
260  */
261  // TODO: When prod(dim(.)) > INT_MAX, we *must* use but cannot --> need itype = CHOLMOD_LONG
262  // ---- CHM_SP chxs = cholmod_l_dense_to_sparse(chxd, 1, &c);
263  CHM_SP chxs = cholmod_dense_to_sparse(chxd, 1, &c);
264
265  int Rkind = (chxd->xtype == CHOLMOD_REAL) ? Real_KIND2(x) : 0;
266  /* Note: when 'x' was integer Matrix, Real_KIND(x) = -1, but *_KIND2(.) = 0 */
267  R_CheckStack();
268
269  UNPROTECT(1);
270  /* chm_sparse_to_SEXP() *could* deal with symmetric
271  * if chxs had such an stype; and we should be able to use uplo below */
272  return chm_sparse_to_SEXP(chxs, 1, 0/*TODO: uplo_P(x) if x has an uplo slot*/,
273  Rkind, "",
274  isMatrix(x) ? getAttrib(x, R_DimNamesSymbol)
275  : GET_SLOT(x, Matrix_DimNamesSym));
276 }
277
278
279 SEXP dense_band(SEXP x, SEXP k1P, SEXP k2P)
280 /* Always returns a full matrix with entries outside the band zeroed
281  * Class of the value can be [dln]trMatrix or [dln]geMatrix
282  */
283 {
284  int k1 = asInteger(k1P), k2 = asInteger(k2P);
285
286  if (k1 > k2) {
287  error(_("Lower band %d > upper band %d"), k1, k2);
288  return R_NilValue; /* -Wall */
289  }
290  else {
291  SEXP ans = PROTECT(dup_mMatrix_as_geMatrix(x));
292  int *adims = INTEGER(GET_SLOT(ans, Matrix_DimSym)),
295  tru = (k1 >= 0), trl = (k2 <= 0);
296  const char *cl = class_P(ans);
297  enum dense_enum M_type = ( (cl[0] == 'd') ? ddense :
298  ((cl[0] == 'l') ? ldense : ndense));
299
300
301 #define SET_ZERO_OUTSIDE \
302  for (j = 0; j < n; j++) { \
303  int i, i1 = j - k2, i2 = j + 1 - k1; \
304  if(i1 > m) i1 = m; \
305  if(i2 < 0) i2 = 0; \
306  for (i = 0; i < i1; i++) xx[i + j * m] = 0; \
307  for (i = i2; i < m; i++) xx[i + j * m] = 0; \
308  }
309
310  if(M_type == ddense) {
311  double *xx = REAL(GET_SLOT(ans, Matrix_xSym));
313  }
314  else { /* (M_type == ldense || M_type == ndense) */
315  int *xx = LOGICAL(GET_SLOT(ans, Matrix_xSym));
317  }
318
319  if (!sqr || (!tru && !trl)) { /* return the *geMatrix */
320  UNPROTECT(1);
321  return ans;
322  }
323  else {
324  /* Copy ans to a *trMatrix object (must be square) */
325  SEXP aa= PROTECT(NEW_OBJECT(MAKE_CLASS(M_type == ddense? "dtrMatrix":
326  (M_type== ldense? "ltrMatrix"
327  : "ntrMatrix"))));
328  /* Because slots of ans are freshly allocated and ans will not be
329  * used, we use the slots themselves and don't duplicate */
330  SET_SLOT(aa, Matrix_xSym, GET_SLOT(ans, Matrix_xSym));
331  SET_SLOT(aa, Matrix_DimSym, GET_SLOT(ans, Matrix_DimSym));
332  SET_SLOT(aa, Matrix_DimNamesSym,GET_SLOT(ans, Matrix_DimNamesSym));
333  SET_SLOT(aa, Matrix_diagSym, mkString("N"));
334  SET_SLOT(aa, Matrix_uploSym, mkString(tru ? "U" : "L"));
335  UNPROTECT(2);
336  return aa;
337  }
338  }
339 }
340
341 SEXP dense_to_symmetric(SEXP x, SEXP uplo, SEXP symm_test)
342 /* Class of result will be [dln]syMatrix */
343 {
344 /*== FIXME: allow uplo = NA and then behave a bit like symmpart():
345  *== ----- would use the *dimnames* to determine U or L (??)
346  */
347
348  int symm_tst = asLogical(symm_test);
349  SEXP dx = PROTECT(dup_mMatrix_as_geMatrix(x));
350  SEXP ans, dns, nms_dns;
351  const char *cl = class_P(dx);
352  /* same as in ..._geMatrix() above:*/
353  enum dense_enum M_type = ( (cl[0] == 'd') ? ddense :
354  ((cl[0] == 'l') ? ldense : ndense));
357  UNPROTECT(1);
358  error(_("ddense_to_symmetric(): matrix is not square!"));
359  return R_NilValue; /* -Wall */
360  }
361
362  if(symm_tst) {
363  int i,j;
364 # define CHECK_SYMMETRIC \
365  for (j = 0; j < n; j++) \
366  for (i = 0; i < j; i++) \
367  if(xx[j * n + i] != xx[i * n + j]) { \
368  UNPROTECT(1); \
369  error(_("matrix is not symmetric [%d,%d]"), i+1, j+1); \
370  return R_NilValue; /* -Wall */ \
371  }
372  if(M_type == ddense) {
373
374  double *xx = REAL(GET_SLOT(dx, Matrix_xSym));
376  }
377  else { /* (M_type == ldense || M_type == ndense) */
378
379  int *xx = LOGICAL(GET_SLOT(dx, Matrix_xSym));
381  }
382  }
383 # undef CHECK_SYMMETRIC
384
385  ans = PROTECT(NEW_OBJECT(MAKE_CLASS( M_type == ddense ? "dsyMatrix" :
386  (M_type == ldense ? "lsyMatrix" :
387  "nsyMatrix"))));
388
389
390 // --- FIXME: Use MK_SYMMETRIC_DIMNAMES_AND_RETURN from below -- with "uplo" argument
391
392  /* need _symmetric_ dimnames */
393  dns = GET_SLOT(dx, Matrix_DimNamesSym);
394  if(!equal_string_vectors(VECTOR_ELT(dns,0),
395  VECTOR_ELT(dns,1))) {
396  if(*CHAR(asChar(uplo)) == 'U')
397  SET_VECTOR_ELT(dns,0, VECTOR_ELT(dns,1));
398  else
399  SET_VECTOR_ELT(dns,1, VECTOR_ELT(dns,0));
400  }
401  if(!isNull(nms_dns = getAttrib(dns, R_NamesSymbol)) &&
402  !R_compute_identical(STRING_ELT(nms_dns, 0),
403  STRING_ELT(nms_dns, 1), 16)) { // names(dimnames(.)) :
404  if(*CHAR(asChar(uplo)) == 'U')
405  SET_STRING_ELT(nms_dns, 0, STRING_ELT(nms_dns,1));
406  else
407  SET_STRING_ELT(nms_dns, 1, STRING_ELT(nms_dns,0));
408  setAttrib(dns, R_NamesSymbol, nms_dns);
409  }
410
411  /* Copy dx to ans;
412  * Because slots of dx are freshly allocated and dx will not be
413  * used, we use the slots themselves and don't duplicate */
414  SET_SLOT(ans, Matrix_xSym, GET_SLOT(dx, Matrix_xSym));
415  SET_SLOT(ans, Matrix_DimSym, GET_SLOT(dx, Matrix_DimSym));
416  SET_SLOT(ans, Matrix_DimNamesSym, dns);
417  SET_SLOT(ans, Matrix_uploSym, ScalarString(asChar(uplo)));
418
419  UNPROTECT(2);
420  return ans;
421 }
422
423 SEXP ddense_symmpart(SEXP x)
424 /* Class of the value will be dsyMatrix */
425 {
426  SEXP dx = dup_mMatrix_as_dgeMatrix(x);
428
430  error(_("matrix is not square! (symmetric part)"));
431  return R_NilValue; /* -Wall */
432  } else {
433  PROTECT(dx);
434  SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dsyMatrix"))), dns, nms_dns;
435  double *xx = REAL(GET_SLOT(dx, Matrix_xSym));
436
437  /* only need to assign the *upper* triangle (uplo = "U");
438  * noting that diagonal remains unchanged */
439  for (int j = 0; j < n; j++) {
440  for (int i = 0; i < j; i++) {
441  xx[j * n + i] = (xx[j * n + i] + xx[i * n + j]) / 2.;
442  }
443  }
444
445 // FIXME?: Compare and synchronize with symmetric_DimNames() in ./Mutils.c
446 # define MK_SYMMETRIC_DIMNAMES_AND_RETURN \
447  \
448  dns = GET_SLOT(dx, Matrix_DimNamesSym); \
449  int J = 1; \
450  if(!equal_string_vectors(VECTOR_ELT(dns,0), \
451  VECTOR_ELT(dns,1))) { \
452  /* _symmetric_ dimnames: behave as symmDN(*, col=TRUE) */ \
453  if(isNull(VECTOR_ELT(dns, J))) \
454  J = !J; \
455  SET_VECTOR_ELT(dns, !J, VECTOR_ELT(dns, J)); \
456  } \
457  /* names(dimnames(.)): */ \
458  if(!isNull(nms_dns = getAttrib(dns, R_NamesSymbol)) && \
459  !R_compute_identical(STRING_ELT(nms_dns, 0), \
460  STRING_ELT(nms_dns, 1), 16)) { \
461  SET_STRING_ELT(nms_dns, !J, STRING_ELT(nms_dns, J)); \
462  setAttrib(dns, R_NamesSymbol, nms_dns); \
463  } \
464  \
465  /* Copy dx to ans; \
466  * Because slots of dx are freshly allocated and dx will not be \
467  * used, we use the slots themselves and don't duplicate */ \
468  \
469  SET_SLOT(ans, Matrix_xSym, GET_SLOT(dx, Matrix_xSym)); \
470  SET_SLOT(ans, Matrix_DimSym, GET_SLOT(dx, Matrix_DimSym)); \
471  SET_SLOT(ans, Matrix_DimNamesSym, dns); \
472  SET_SLOT(ans, Matrix_uploSym, mkString("U")); \
473  \
474  UNPROTECT(2); \
475  return ans
476
478  }
479 }
480
481 SEXP ddense_skewpart(SEXP x)
482 /* Class of the value will be dgeMatrix */
483 {
484  SEXP dx = dup_mMatrix_as_dgeMatrix(x);
486
488  error(_("matrix is not square! (skew-symmetric part)"));
489  return R_NilValue; /* -Wall */
490  } else {
491  PROTECT(dx);
492  SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix"))), dns, nms_dns;
493  double *xx = REAL(GET_SLOT(dx, Matrix_xSym));
494
495  for (int j = 0; j < n; j++) {
496  xx[j * n + j] = 0.;
497  for (int i = 0; i < j; i++) {
498  double s = (xx[j * n + i] - xx[i * n + j]) / 2.;
499  xx[j * n + i] = s;
500  xx[i * n + j] = -s;
501  }
502  }
503
505  }
506 }
SEXP Matrix_DimSym
Definition: Syms.h:2
#define class_P(_x_)
Definition: Mutils.h:178
cholmod_sparse * CHM_SP
Definition: chm_common.h:25
#define Real_KIND2(_x_)
Definition: Mutils.h:191
SEXP Matrix_xSym
Definition: Syms.h:2
SEXP dense_to_symmetric(SEXP x, SEXP uplo, SEXP symm_test)
Definition: dense.c:341
SEXP Matrix_DimNamesSym
Definition: Syms.h:2
SEXP chm_sparse_to_SEXP(CHM_SP a, int dofree, int uploT, int Rkind, const char *diag, SEXP dn)
Copy the contents of a to an appropriate CsparseMatrix object and, optionally, free a or free both a ...
Definition: chm_common.c:335
cholmod_common cl
Definition: chm_common.c:16
SEXP Matrix_uploSym
Definition: Syms.h:2
SEXP dup_mMatrix_as_dgeMatrix(SEXP A)
Definition: Mutils.c:852
SEXP dense_band(SEXP x, SEXP k1P, SEXP k2P)
Definition: dense.c:279
dense_enum
Definition: Mutils.h:181
SEXP lapack_qr(SEXP Xin, SEXP tl)
Definition: dense.c:171
static SEXP getGivens(double x[], int ldx, int jmin, int rank)
Definition: dense.c:63
Rboolean equal_string_vectors(SEXP s1, SEXP s2)
Definition: Mutils.c:286
#define MK_SYMMETRIC_DIMNAMES_AND_RETURN
#define CHECK_SYMMETRIC
Definition: Mutils.h:181
static R_INLINE SEXP mMatrix_as_geMatrix(SEXP A)
Definition: Mutils.h:343
#define AS_CHM_xDN(x)
Definition: chm_common.h:44
SEXP checkGivens(SEXP X, SEXP jmin, SEXP rank)
Definition: dense.c:83
SEXP lsq_dense_QR(SEXP X, SEXP y)
Definition: dense.c:132
SEXP lsq_dense_Chol(SEXP X, SEXP y)
Definition: dense.c:99
#define _(String)
Definition: Mutils.h:32
SEXP ddense_skewpart(SEXP x)
Definition: dense.c:481
cholmod_dense * CHM_DN
Definition: chm_common.h:21
static void * xpt(int ctype, SEXP x)
Definition: chm_common.c:154
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 dense_to_Csparse(SEXP x)
Definition: dense.c:248
Definition: Mutils.h:181
static int left_cyclic(double x[], int ldx, int j, int k, double cosines[], double sines[])
Perform a left cyclic shift of columns j to k in the upper triangular matrix x, then restore it to up...
Definition: dense.c:22
SEXP Matrix_diagSym
Definition: Syms.h:2
Definition: Mutils.h:181
cholmod_common c
Definition: chm_common.c:15
#define SET_ZERO_OUTSIDE
SEXP ddense_symmpart(SEXP x)
Definition: dense.c:423