Matrix  $Rev: 3071 $ at $LastChangedDate: 2015-03-26 15:35:47 +0100 (Thu, 26 Mar 2015) $
dppMatrix.c
Go to the documentation of this file.
1 #include "dppMatrix.h"
2 
3 SEXP dppMatrix_validate(SEXP obj)
4 {
5 /* int i, n = INTEGER(GET_SLOT(obj, Matrix_DimSym))[0]; */
6 /* double *x = REAL(GET_SLOT(obj, Matrix_xSym)); */
7 
8  /* quick but nondefinitive check on positive definiteness */
9 /* for (i = 0; i < n; i++) */
10 /* if (x[i * np1] < 0) */
11 /* return mkString(_("dppMatrix is not positive definite")); */
12  return dspMatrix_validate(obj);
13 }
14 
15 SEXP dppMatrix_chol(SEXP x)
16 {
17  SEXP val = get_factors(x, "pCholesky"),
18  dimP = GET_SLOT(x, Matrix_DimSym),
19  uploP = GET_SLOT(x, Matrix_uploSym);
20  const char *uplo = CHAR(STRING_ELT(uploP, 0));
21  int *dims = INTEGER(dimP), info;
22 
23  if (val != R_NilValue) return val;
24  dims = INTEGER(dimP);
25  val = PROTECT(NEW_OBJECT(MAKE_CLASS("pCholesky")));
26  SET_SLOT(val, Matrix_uploSym, duplicate(uploP));
27  SET_SLOT(val, Matrix_diagSym, mkString("N"));
28  SET_SLOT(val, Matrix_DimSym, duplicate(dimP));
29  slot_dup(val, x, Matrix_xSym);
30  F77_CALL(dpptrf)(uplo, dims, REAL(GET_SLOT(val, Matrix_xSym)), &info);
31  if (info) {
32  if(info > 0) /* e.g. x singular */
33  error(_("the leading minor of order %d is not positive definite"),
34  info);
35  else /* should never happen! */
36  error(_("Lapack routine %s returned error code %d"), "dpptrf", info);
37  }
38  UNPROTECT(1);
39  return set_factors(x, val, "pCholesky");
40 }
41 
42 SEXP dppMatrix_rcond(SEXP obj, SEXP type)
43 {
44  SEXP Chol = dppMatrix_chol(obj);
45  char typnm[] = {'O', '\0'}; /* always use the one norm */
46  int *dims = INTEGER(GET_SLOT(Chol, Matrix_DimSym)), info;
47  double anorm = get_norm_sp(obj, typnm), rcond;
48 
49  F77_CALL(dppcon)(uplo_P(Chol), dims,
50  REAL(GET_SLOT(Chol, Matrix_xSym)), &anorm, &rcond,
51  (double *) R_alloc(3*dims[0], sizeof(double)),
52  (int *) R_alloc(dims[0], sizeof(int)), &info);
53  return ScalarReal(rcond);
54 }
55 
56 SEXP dppMatrix_solve(SEXP x)
57 {
58  SEXP Chol = dppMatrix_chol(x);
59  SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dppMatrix")));
60  int *dims = INTEGER(GET_SLOT(x, Matrix_DimSym)), info;
61 
62  slot_dup(val, Chol, Matrix_uploSym);
63  slot_dup(val, Chol, Matrix_xSym);
64  slot_dup(val, Chol, Matrix_DimSym);
65  F77_CALL(dpptri)(uplo_P(val), dims,
66  REAL(GET_SLOT(val, Matrix_xSym)), &info);
67  UNPROTECT(1);
68  return val;
69 }
70 
71 SEXP dppMatrix_matrix_solve(SEXP a, SEXP b)
72 {
73  SEXP val = PROTECT(dup_mMatrix_as_dgeMatrix(b));
74  SEXP Chol = dppMatrix_chol(a);
75  int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
76  *bdims = INTEGER(GET_SLOT(val, Matrix_DimSym));
77  int n = bdims[0], nrhs = bdims[1], info;
78 
79  if (*adims != *bdims || bdims[1] < 1 || *adims < 1)
80  error(_("Dimensions of system to be solved are inconsistent"));
81  F77_CALL(dpptrs)(uplo_P(Chol), &n, &nrhs,
82  REAL(GET_SLOT(Chol, Matrix_xSym)),
83  REAL(GET_SLOT(val, Matrix_xSym)), &n, &info);
84  UNPROTECT(1);
85  return val;
86 }
SEXP Matrix_DimSym
Definition: Syms.h:2
SEXP dppMatrix_solve(SEXP x)
Definition: dppMatrix.c:56
SEXP Matrix_xSym
Definition: Syms.h:2
#define slot_dup(dest, src, sym)
Definition: Mutils.h:149
SEXP dspMatrix_validate(SEXP obj)
Definition: dspMatrix.c:4
SEXP Matrix_uploSym
Definition: Syms.h:2
SEXP dup_mMatrix_as_dgeMatrix(SEXP A)
Definition: Mutils.c:852
SEXP dppMatrix_matrix_solve(SEXP a, SEXP b)
Definition: dppMatrix.c:71
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
#define uplo_P(_x_)
Definition: Mutils.h:174
#define _(String)
Definition: Mutils.h:32
SEXP dppMatrix_chol(SEXP x)
Definition: dppMatrix.c:15
SEXP Matrix_diagSym
Definition: Syms.h:2
double get_norm_sp(SEXP obj, const char *typstr)
Definition: dspMatrix.c:18
SEXP dppMatrix_validate(SEXP obj)
Definition: dppMatrix.c:3
SEXP dppMatrix_rcond(SEXP obj, SEXP type)
Definition: dppMatrix.c:42