Matrix  $Rev: 3071 $ at $LastChangedDate: 2015-03-26 15:35:47 +0100 (Thu, 26 Mar 2015) $
dpoMatrix.c
Go to the documentation of this file.
1 #include "dpoMatrix.h"
2 
3 SEXP dpoMatrix_validate(SEXP obj)
4 {
5  SEXP val;
6  if (isString(val = dense_nonpacked_validate(obj)))
7  return(val);
8 
9  int n = INTEGER(GET_SLOT(obj, Matrix_DimSym))[0],
10  np1 = n + 1;
11  double *x = REAL(GET_SLOT(obj, Matrix_xSym));
12 
13  /* quick but nondefinitive check on positive definiteness */
14  for (int i = 0; i < n; i++)
15  if (x[i * np1] < 0)
16  return mkString(_("dpoMatrix is not positive definite"));
17  return ScalarLogical(1);
18 }
19 
20 SEXP dpoMatrix_chol(SEXP x)
21 {
22  SEXP val = get_factors(x, "Cholesky"),
23  dimP = GET_SLOT(x, Matrix_DimSym),
24  uploP = GET_SLOT(x, Matrix_uploSym);
25  const char *uplo = CHAR(STRING_ELT(uploP, 0));
26  int *dims = INTEGER(dimP), info;
27  int n = dims[0];
28  double *vx;
29 
30  if (val != R_NilValue) return val;// use x@factors$Cholesky if available
31  dims = INTEGER(dimP);
32  val = PROTECT(NEW_OBJECT(MAKE_CLASS("Cholesky")));
33  SET_SLOT(val, Matrix_uploSym, duplicate(uploP));
34  SET_SLOT(val, Matrix_diagSym, mkString("N"));
35  SET_SLOT(val, Matrix_DimSym, duplicate(dimP));
36  vx = REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, n * n));
37  AZERO(vx, n * n);
38  F77_CALL(dlacpy)(uplo, &n, &n, REAL(GET_SLOT(x, Matrix_xSym)), &n, vx, &n);
39  if (n > 0) {
40  F77_CALL(dpotrf)(uplo, &n, vx, &n, &info);
41  if (info) {
42  if(info > 0)
43  error(_("the leading minor of order %d is not positive definite"),
44  info);
45  else /* should never happen! */
46  error(_("Lapack routine %s returned error code %d"), "dpotrf", info);
47  }
48  }
49  UNPROTECT(1);
50  return set_factors(x, val, "Cholesky");
51 }
52 
53 SEXP dpoMatrix_rcond(SEXP obj, SEXP type)
54 {
55  SEXP Chol = dpoMatrix_chol(obj);
56  const char typnm[] = {'O', '\0'}; /* always use the one norm */
57  int *dims = INTEGER(GET_SLOT(Chol, Matrix_DimSym)), info;
58  double anorm = get_norm_sy(obj, typnm), rcond;
59 
60  F77_CALL(dpocon)(uplo_P(Chol),
61  dims, REAL(GET_SLOT(Chol, Matrix_xSym)),
62  dims, &anorm, &rcond,
63  (double *) R_alloc(3*dims[0], sizeof(double)),
64  (int *) R_alloc(dims[0], sizeof(int)), &info);
65  return ScalarReal(rcond);
66 }
67 
68 SEXP dpoMatrix_solve(SEXP x)
69 {
70  SEXP Chol = dpoMatrix_chol(x);
71  SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dpoMatrix")));
72  int *dims = INTEGER(GET_SLOT(x, Matrix_DimSym)), info;
73 
74  SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));
75  slot_dup(val, Chol, Matrix_uploSym);
76  slot_dup(val, Chol, Matrix_xSym);
77  slot_dup(val, Chol, Matrix_DimSym);
78  SET_SLOT(val, Matrix_DimNamesSym,
79  duplicate(GET_SLOT(x, Matrix_DimNamesSym)));
80  F77_CALL(dpotri)(uplo_P(val), dims,
81  REAL(GET_SLOT(val, Matrix_xSym)), dims, &info);
82  UNPROTECT(1);
83  return val;
84 }
85 
86 SEXP dpoMatrix_dgeMatrix_solve(SEXP a, SEXP b)
87 {
88  SEXP Chol = dpoMatrix_chol(a),
89  val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
90  int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
91  *bdims = INTEGER(GET_SLOT(b, Matrix_DimSym)),
92  info;
93 
94  if (adims[1] != bdims[0])
95  error(_("Dimensions of system to be solved are inconsistent"));
96  if (adims[0] < 1 || bdims[1] < 1)
97  error(_("Cannot solve() for matrices with zero extents"));
98  SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));
99  slot_dup(val, b, Matrix_DimSym);
100  slot_dup(val, b, Matrix_xSym);
101  F77_CALL(dpotrs)(uplo_P(Chol), adims, bdims + 1,
102  REAL(GET_SLOT(Chol, Matrix_xSym)), adims,
103  REAL(GET_SLOT(val, Matrix_xSym)),
104  bdims, &info);
105  UNPROTECT(1);
106  return val;
107 }
108 
109 SEXP dpoMatrix_matrix_solve(SEXP a, SEXP b)
110 {
111  SEXP Chol = dpoMatrix_chol(a),
112  val = PROTECT(duplicate(b));
113  int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
114  *bdims = INTEGER(getAttrib(b, R_DimSymbol)),
115  info;
116 
117  if (!(isReal(b) && isMatrix(b)))
118  error(_("Argument b must be a numeric matrix"));
119  if (*adims != *bdims || bdims[1] < 1 || *adims < 1)
120  error(_("Dimensions of system to be solved are inconsistent"));
121  F77_CALL(dpotrs)(uplo_P(Chol), adims, bdims + 1,
122  REAL(GET_SLOT(Chol, Matrix_xSym)), adims,
123  REAL(val), bdims, &info);
124  UNPROTECT(1);
125  return val;
126 }
SEXP Matrix_DimSym
Definition: Syms.h:2
SEXP Matrix_xSym
Definition: Syms.h:2
double get_norm_sy(SEXP obj, const char *typstr)
Definition: dsyMatrix.c:15
#define slot_dup(dest, src, sym)
Definition: Mutils.h:149
SEXP Matrix_factorSym
Definition: Syms.h:2
SEXP Matrix_DimNamesSym
Definition: Syms.h:2
SEXP dpoMatrix_matrix_solve(SEXP a, SEXP b)
Definition: dpoMatrix.c:109
SEXP Matrix_uploSym
Definition: Syms.h:2
SEXP set_factors(SEXP obj, SEXP val, char *nm)
Caches 'val' in the 'factors' slot of obj, i.e.
Definition: Mutils.c:130
SEXP get_factors(SEXP obj, char *nm)
Definition: Mutils.c:106
SEXP dpoMatrix_solve(SEXP x)
Definition: dpoMatrix.c:68
#define AZERO(x, n)
Definition: Mutils.h:140
#define uplo_P(_x_)
Definition: Mutils.h:174
#define _(String)
Definition: Mutils.h:32
SEXP dpoMatrix_chol(SEXP x)
Definition: dpoMatrix.c:20
SEXP dpoMatrix_validate(SEXP obj)
Definition: dpoMatrix.c:3
SEXP Matrix_diagSym
Definition: Syms.h:2
SEXP dpoMatrix_dgeMatrix_solve(SEXP a, SEXP b)
Definition: dpoMatrix.c:86
SEXP dpoMatrix_rcond(SEXP obj, SEXP type)
Definition: dpoMatrix.c:53
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 dense_nonpacked_validate(SEXP obj)
Definition: Mutils.c:310