Matrix  $Rev: 3071 $ at $LastChangedDate: 2015-03-26 15:35:47 +0100 (Thu, 26 Mar 2015) $
dsyMatrix.c
Go to the documentation of this file.
1 #include "dsyMatrix.h"
2 
4 {
5  SEXP val = GET_SLOT(obj, Matrix_DimSym);
6  if (LENGTH(val) < 2)
7  return mkString(_("'Dim' slot has length less than two"));
8  if (INTEGER(val)[0] != INTEGER(val)[1])
9  return mkString(_("Matrix is not square"));
10  if (isString(val = check_scalar_string(GET_SLOT(obj, Matrix_uploSym),
11  "LU", "uplo"))) return val;
12  return ScalarLogical(1);
13 }
14 
15 double get_norm_sy(SEXP obj, const char *typstr)
16 {
17  char typnm[] = {'\0', '\0'};
18  int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym));
19  double *work = (double *) NULL;
20 
21  typnm[0] = La_norm_type(typstr);
22  if (*typnm == 'I' || *typnm == 'O') {
23  work = (double *) R_alloc(dims[0], sizeof(double));
24  }
25  return F77_CALL(dlansy)(typnm, uplo_P(obj),
26  dims, REAL(GET_SLOT(obj, Matrix_xSym)),
27  dims, work);
28 }
29 
30 SEXP dsyMatrix_norm(SEXP obj, SEXP type)
31 {
32  return ScalarReal(get_norm_sy(obj, CHAR(asChar(type))));
33 }
34 
35 
36 SEXP dsyMatrix_rcond(SEXP obj, SEXP type)
37 {
38  SEXP trf = dsyMatrix_trf(obj);
39  int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym)), info;
40  double anorm = get_norm_sy(obj, "O");
41  double rcond;
42 
43  F77_CALL(dsycon)(uplo_P(trf), dims,
44  REAL (GET_SLOT(trf, Matrix_xSym)), dims,
45  INTEGER(GET_SLOT(trf, Matrix_permSym)),
46  &anorm, &rcond,
47  (double *) R_alloc(2*dims[0], sizeof(double)),
48  (int *) R_alloc(dims[0], sizeof(int)), &info);
49  return ScalarReal(rcond);
50 }
51 
52 SEXP dsyMatrix_solve(SEXP a)
53 {
54  SEXP trf = dsyMatrix_trf(a);
55  SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dsyMatrix")));
56  int *dims = INTEGER(GET_SLOT(trf, Matrix_DimSym)), info;
57 
58  slot_dup(val, trf, Matrix_uploSym);
59  slot_dup(val, trf, Matrix_xSym);
60  slot_dup(val, trf, Matrix_DimSym);
61  F77_CALL(dsytri)(uplo_P(val), dims,
62  REAL(GET_SLOT(val, Matrix_xSym)), dims,
63  INTEGER(GET_SLOT(trf, Matrix_permSym)),
64  (double *) R_alloc((long) dims[0], sizeof(double)),
65  &info);
66  UNPROTECT(1);
67  return val;
68 }
69 
70 SEXP dsyMatrix_matrix_solve(SEXP a, SEXP b)
71 {
72  SEXP trf = dsyMatrix_trf(a),
73  val = PROTECT(dup_mMatrix_as_dgeMatrix(b));
74  int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
75  *bdims = INTEGER(GET_SLOT(val, Matrix_DimSym)),
76  info;
77 
78  if (*adims != *bdims || bdims[1] < 1 || *adims < 1)
79  error(_("Dimensions of system to be solved are inconsistent"));
80  F77_CALL(dsytrs)(uplo_P(trf), adims, bdims + 1,
81  REAL(GET_SLOT(trf, Matrix_xSym)), adims,
82  INTEGER(GET_SLOT(trf, Matrix_permSym)),
83  REAL(GET_SLOT(val, Matrix_xSym)),
84  bdims, &info);
85  UNPROTECT(1);
86  return val;
87 }
88 
89 SEXP dsyMatrix_as_matrix(SEXP from, SEXP keep_dimnames)
90 {
91  int n = INTEGER(GET_SLOT(from, Matrix_DimSym))[0];
92  SEXP val = PROTECT(allocMatrix(REALSXP, n, n));
93 
94  make_d_matrix_symmetric(Memcpy(REAL(val),
95  REAL(GET_SLOT(from, Matrix_xSym)), n * n),
96  from);
97  if(asLogical(keep_dimnames))
98  setAttrib(val, R_DimNamesSymbol, R_symmetric_Dimnames(from));
99  UNPROTECT(1);
100  return val;
101 }
102 
103 SEXP dsyMatrix_matrix_mm(SEXP a, SEXP b, SEXP rtP)
104 {
105  SEXP val = PROTECT(dup_mMatrix_as_dgeMatrix(b));// incl. its dimnames
106  int rt = asLogical(rtP); /* if(rt), compute b %*% a, else a %*% b */
107  int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
108  *bdims = INTEGER(GET_SLOT(val, Matrix_DimSym)),
109  m = bdims[0], n = bdims[1];
110  double one = 1., zero = 0., mn = ((double) m) * ((double) n);
111  if (mn > INT_MAX)
112  error(_("Matrix dimension %d x %d (= %g) is too large"), m, n, mn);
113  // else: m * n will not overflow below
114  double *bcp, *vx = REAL(GET_SLOT(val, Matrix_xSym));
115  C_or_Alloca_TO(bcp, m * n, double);
116  Memcpy(bcp, vx, m * n);
117 
118  if ((rt && n != adims[0]) || (!rt && m != adims[0]))
119  error(_("Matrices are not conformable for multiplication"));
120  if (m >=1 && n >= 1)
121  F77_CALL(dsymm)(rt ? "R" :"L", uplo_P(a), &m, &n, &one,
122  REAL(GET_SLOT(a, Matrix_xSym)), adims, bcp,
123  &m, &zero, vx, &m);
124  // add dimnames:
125  if(rt) { // v <- b %*% a : rownames(v) == rownames(b) are already there
126  SET_VECTOR_ELT(GET_SLOT(val, Matrix_DimNamesSym), 1,
127  duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 1)));
128  } else { // v <- a %*% b : colnames(v) == colnames(b) are already there
129  SET_VECTOR_ELT(GET_SLOT(val, Matrix_DimNamesSym), 0,
130  duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 0)));
131  }
132  if(mn >= SMALL_4_Alloca) Free(bcp);
133  UNPROTECT(1);
134  return val;
135 }
136 
137 SEXP dsyMatrix_trf(SEXP x)
138 {
139  SEXP val = get_factors(x, "BunchKaufman"),
140  dimP = GET_SLOT(x, Matrix_DimSym),
141  uploP = GET_SLOT(x, Matrix_uploSym);
142  int *dims = INTEGER(dimP), *perm, info;
143  int lwork = -1, n = dims[0];
144  const char *uplo = CHAR(STRING_ELT(uploP, 0));
145  double tmp, *vx, *work;
146 
147  if (val != R_NilValue) return val;
148  dims = INTEGER(dimP);
149  val = PROTECT(NEW_OBJECT(MAKE_CLASS("BunchKaufman")));
150  SET_SLOT(val, Matrix_uploSym, duplicate(uploP));
151  SET_SLOT(val, Matrix_diagSym, mkString("N"));
152  SET_SLOT(val, Matrix_DimSym, duplicate(dimP));
153  vx = REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, n * n));
154  AZERO(vx, n * n);
155  F77_CALL(dlacpy)(uplo, &n, &n, REAL(GET_SLOT(x, Matrix_xSym)), &n, vx, &n);
156  perm = INTEGER(ALLOC_SLOT(val, Matrix_permSym, INTSXP, n));
157  F77_CALL(dsytrf)(uplo, &n, vx, &n, perm, &tmp, &lwork, &info);
158  lwork = (int) tmp;
159  C_or_Alloca_TO(work, lwork, double);
160 
161  F77_CALL(dsytrf)(uplo, &n, vx, &n, perm, work, &lwork, &info);
162 
163  if(lwork >= SMALL_4_Alloca) Free(work);
164  if (info) error(_("Lapack routine dsytrf returned error code %d"), info);
165  UNPROTECT(1);
166  return set_factors(x, val, "BunchKaufman");
167 }
168 
169 // this is very close to lsyMatrix_as_lsp*() in ./ldense.c -- keep synced !
170 SEXP dsyMatrix_as_dspMatrix(SEXP from)
171 {
172  SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dspMatrix"))),
173  uplo = GET_SLOT(from, Matrix_uploSym),
174  dimP = GET_SLOT(from, Matrix_DimSym);
175  int n = *INTEGER(dimP);
176 
177  SET_SLOT(val, Matrix_DimSym, duplicate(dimP));
178  SET_SLOT(val, Matrix_uploSym, duplicate(uplo));
179  full_to_packed_double(
180  REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, (n*(n+1))/2)),
181  REAL( GET_SLOT(from, Matrix_xSym)), n,
182  *CHAR(STRING_ELT(uplo, 0)) == 'U' ? UPP : LOW, NUN);
183  SET_SLOT(val, Matrix_DimNamesSym,
184  duplicate(GET_SLOT(from, Matrix_DimNamesSym)));
185  SET_SLOT(val, Matrix_factorSym,
186  duplicate(GET_SLOT(from, Matrix_factorSym)));
187  UNPROTECT(1);
188  return val;
189 }
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 slot_dup(dest, src, sym)
Definition: Mutils.h:149
SEXP Matrix_factorSym
Definition: Syms.h:2
SEXP Matrix_DimNamesSym
Definition: Syms.h:2
SEXP dsyMatrix_as_dspMatrix(SEXP from)
Definition: dsyMatrix.c:170
SEXP Matrix_uploSym
Definition: Syms.h:2
SEXP dup_mMatrix_as_dgeMatrix(SEXP A)
Definition: Mutils.c:852
SEXP dsyMatrix_matrix_solve(SEXP a, SEXP b)
Definition: dsyMatrix.c:70
SEXP check_scalar_string(SEXP sP, char *vals, char *nm)
Check validity of 1-letter string from a set of possible values (typically used in S4 validity method...
Definition: Mutils.c:254
#define SMALL_4_Alloca
Definition: Mutils.h:47
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 UPP
Definition: Mutils.h:81
double get_norm_sy(SEXP obj, const char *typstr)
Definition: dsyMatrix.c:15
#define AZERO(x, n)
Definition: Mutils.h:140
SEXP Matrix_permSym
Definition: Syms.h:2
#define uplo_P(_x_)
Definition: Mutils.h:174
#define _(String)
Definition: Mutils.h:32
char La_norm_type(const char *typstr)
Definition: Mutils.c:8
#define LOW
Definition: Mutils.h:82
SEXP dsyMatrix_matrix_mm(SEXP a, SEXP b, SEXP rtP)
Definition: dsyMatrix.c:103
void make_d_matrix_symmetric(double *to, SEXP from) void make_i_matrix_symmetric(int *to
SEXP Matrix_diagSym
Definition: Syms.h:2
SEXP symmetricMatrix_validate(SEXP obj)
Definition: dsyMatrix.c:3
#define NUN
Definition: Mutils.h:83
SEXP dsyMatrix_norm(SEXP obj, SEXP type)
Definition: dsyMatrix.c:30
SEXP dsyMatrix_rcond(SEXP obj, SEXP type)
Definition: dsyMatrix.c:36
SEXP dsyMatrix_solve(SEXP a)
Definition: dsyMatrix.c:52
SEXP dsyMatrix_as_matrix(SEXP from, SEXP keep_dimnames)
Definition: dsyMatrix.c:89
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 R_symmetric_Dimnames(SEXP x)
Even if the Dimnames slot is list(NULL, ) etc, return symmetric dimnames: Get ...
Definition: Mutils.c:1315
SEXP dsyMatrix_trf(SEXP x)
Definition: dsyMatrix.c:137