Matrix  $Rev: 3071 $ at $LastChangedDate: 2015-03-26 15:35:47 +0100 (Thu, 26 Mar 2015) $
dspMatrix.c
Go to the documentation of this file.
1 #include "dspMatrix.h"
2 
3 /* Note: also used for lspMatrix */
4 SEXP dspMatrix_validate(SEXP obj)
5 {
6  SEXP val = symmetricMatrix_validate(obj);
7  if(isString(val))
8  return(val);
9  else { /* identical to the test in dtpMatrix_validate() : */
10  int d = INTEGER(GET_SLOT(obj, Matrix_DimSym))[0],
11  lx = length(GET_SLOT(obj, Matrix_xSym));
12  if(lx * 2 != d*(d+1))
13  return(mkString(_("Incorrect length of 'x' slot")));
14  return ScalarLogical(1);
15  }
16 }
17 
18 double get_norm_sp(SEXP obj, const char *typstr)
19 {
20  char typnm[] = {'\0', '\0'};
21  int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym));
22  double *work = (double *) NULL;
23 
24  typnm[0] = La_norm_type(typstr);
25  if (*typnm == 'I' || *typnm == 'O') {
26  work = (double *) R_alloc(dims[0], sizeof(double));
27  }
28  return F77_CALL(dlansp)(typnm, uplo_P(obj), dims,
29  REAL(GET_SLOT(obj, Matrix_xSym)), work);
30 }
31 
32 SEXP dspMatrix_norm(SEXP obj, SEXP type)
33 {
34  return ScalarReal(get_norm_sp(obj, CHAR(asChar(type))));
35 }
36 
37 SEXP dspMatrix_rcond(SEXP obj, SEXP type)
38 {
39  SEXP trf = dspMatrix_trf(obj);
40  int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym)), info;
41  double anorm = get_norm_sp(obj, "O"), rcond;
42 
43  F77_CALL(dspcon)(uplo_P(trf), dims,
44  REAL (GET_SLOT(trf, Matrix_xSym)),
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 dspMatrix_solve(SEXP a)
53 {
54  SEXP trf = dspMatrix_trf(a);
55  SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dspMatrix")));
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(dsptri)(uplo_P(val), dims, REAL(GET_SLOT(val, Matrix_xSym)),
62  INTEGER(GET_SLOT(trf, Matrix_permSym)),
63  (double *) R_alloc((long) dims[0], sizeof(double)),
64  &info);
65  UNPROTECT(1);
66  return val;
67 }
68 
69 SEXP dspMatrix_matrix_solve(SEXP a, SEXP b)
70 {
71  SEXP trf = dspMatrix_trf(a),
72  val = PROTECT(dup_mMatrix_as_dgeMatrix(b));
73  int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
74  *bdims = INTEGER(GET_SLOT(val, Matrix_DimSym));
75  int n = bdims[0], nrhs = bdims[1], info;
76 
77  if (adims[0] != n || nrhs < 1 || n < 1)
78  error(_("Dimensions of system to be solved are inconsistent"));
79  F77_CALL(dsptrs)(uplo_P(trf),
80  &n, &nrhs, REAL(GET_SLOT(trf, Matrix_xSym)),
81  INTEGER(GET_SLOT(trf, Matrix_permSym)),
82  REAL(GET_SLOT(val, Matrix_xSym)), &n, &info);
83  UNPROTECT(1);
84  return val;
85 }
86 
87 SEXP dspMatrix_getDiag(SEXP x)
88 
89 {
90  int n = *INTEGER(GET_SLOT(x, Matrix_DimSym));
91  SEXP val = PROTECT(allocVector(REALSXP, n));
92 
93  d_packed_getDiag(REAL(val), x, n);
94  UNPROTECT(1);
95  return val;
96 }
97 
98 SEXP lspMatrix_getDiag(SEXP x)
99 {
100  int n = *INTEGER(GET_SLOT(x, Matrix_DimSym));
101  SEXP val = PROTECT(allocVector(LGLSXP, n));
102 
103  l_packed_getDiag(LOGICAL(val), x, n);
104  UNPROTECT(1);
105  return val;
106 }
107 
108 SEXP dspMatrix_setDiag(SEXP x, SEXP d)
109 {
110  int n = INTEGER(GET_SLOT(x, Matrix_DimSym))[0];
111  return d_packed_setDiag(REAL(d), LENGTH(d), x, n);
112 }
113 
114 SEXP lspMatrix_setDiag(SEXP x, SEXP d)
115 {
116  int n = INTEGER(GET_SLOT(x, Matrix_DimSym))[0];
117  return l_packed_setDiag(INTEGER(d), LENGTH(d), x, n);
118 }
119 
120 
121 SEXP dspMatrix_as_dsyMatrix(SEXP from)
122 {
123  SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dsyMatrix"))),
124  uplo = GET_SLOT(from, Matrix_uploSym),
125  dimP = GET_SLOT(from, Matrix_DimSym),
126  dmnP = GET_SLOT(from, Matrix_DimNamesSym);
127  int n = *INTEGER(dimP);
128 
129  SET_SLOT(val, Matrix_DimSym, duplicate(dimP));
130  SET_SLOT(val, Matrix_DimNamesSym, duplicate(dmnP));
131  SET_SLOT(val, Matrix_uploSym, duplicate(uplo));
132  packed_to_full_double(REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, n*n)),
133  REAL(GET_SLOT(from, Matrix_xSym)), n,
134  *CHAR(STRING_ELT(uplo, 0)) == 'U' ? UPP : LOW);
135  UNPROTECT(1);
136  return val;
137 }
138 
139 SEXP dspMatrix_matrix_mm(SEXP a, SEXP b)
140 {
141  SEXP val = PROTECT(dup_mMatrix_as_dgeMatrix(b));
142  int *bdims = INTEGER(GET_SLOT(val, Matrix_DimSym));
143  int i, ione = 1, n = bdims[0], nrhs = bdims[1];
144  double nn = ((double) n) * ((double) nrhs);
145  if (nn > INT_MAX)
146  error(_("Matrix dimension %d x %d (= %g) is too large"), n, nrhs, nn);
147  const char *uplo = uplo_P(a);
148  double *ax = REAL(GET_SLOT(a, Matrix_xSym)), one = 1., zero = 0.,
149  *vx = REAL(GET_SLOT(val, Matrix_xSym)), *bx;
150  C_or_Alloca_TO(bx, n * nrhs, double);
151 
152  Memcpy(bx, vx, n * nrhs);
153  if (bdims[0] != n)
154  error(_("Matrices are not conformable for multiplication"));
155  if (nrhs >= 1 && n >= 1) {
156  for (i = 0; i < nrhs; i++)
157  F77_CALL(dspmv)(uplo, &n, &one, ax, bx + i * n, &ione,
158  &zero, vx + i * n, &ione);
159  if(n * nrhs >= SMALL_4_Alloca) Free(bx);
160  }
161  UNPROTECT(1);
162  return val;
163 }
164 
165 SEXP dspMatrix_trf(SEXP x)
166 {
167  SEXP val = get_factors(x, "pBunchKaufman"),
168  dimP = GET_SLOT(x, Matrix_DimSym),
169  uploP = GET_SLOT(x, Matrix_uploSym);
170  int *dims = INTEGER(dimP), *perm, info;
171  int n = dims[0];
172  const char *uplo = CHAR(STRING_ELT(uploP, 0));
173 
174  if (val != R_NilValue) return val;
175  dims = INTEGER(dimP);
176  val = PROTECT(NEW_OBJECT(MAKE_CLASS("pBunchKaufman")));
177  SET_SLOT(val, Matrix_uploSym, duplicate(uploP));
178  SET_SLOT(val, Matrix_diagSym, mkString("N"));
179  SET_SLOT(val, Matrix_DimSym, duplicate(dimP));
180  slot_dup(val, x, Matrix_xSym);
181  perm = INTEGER(ALLOC_SLOT(val, Matrix_permSym, INTSXP, n));
182  F77_CALL(dsptrf)(uplo, dims, REAL(GET_SLOT(val, Matrix_xSym)), perm, &info);
183  if (info) error(_("Lapack routine %s returned error code %d"), "dsptrf", info);
184  UNPROTECT(1);
185  return set_factors(x, val, "pBunchKaufman");
186 }
187 
SEXP Matrix_DimSym
Definition: Syms.h:2
#define C_or_Alloca_TO(_VAR_, _N_, _TYPE_)
Definition: Mutils.h:50
SEXP dspMatrix_trf(SEXP x)
Definition: dspMatrix.c:165
SEXP Matrix_xSym
Definition: Syms.h:2
#define slot_dup(dest, src, sym)
Definition: Mutils.h:149
SEXP Matrix_DimNamesSym
Definition: Syms.h:2
SEXP dspMatrix_validate(SEXP obj)
Definition: dspMatrix.c:4
void l_packed_getDiag(int *dest, SEXP x, int n)
Definition: Mutils.c:450
SEXP Matrix_uploSym
Definition: Syms.h:2
SEXP dup_mMatrix_as_dgeMatrix(SEXP A)
Definition: Mutils.c:852
SEXP dspMatrix_setDiag(SEXP x, SEXP d)
Definition: dspMatrix.c:108
SEXP dspMatrix_getDiag(SEXP x)
Definition: dspMatrix.c:87
#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 lspMatrix_getDiag(SEXP x)
Definition: dspMatrix.c:98
SEXP get_factors(SEXP obj, char *nm)
Definition: Mutils.c:106
#define UPP
Definition: Mutils.h:81
SEXP dspMatrix_rcond(SEXP obj, SEXP type)
Definition: dspMatrix.c:37
SEXP Matrix_permSym
Definition: Syms.h:2
#define uplo_P(_x_)
Definition: Mutils.h:174
#define _(String)
Definition: Mutils.h:32
SEXP d_packed_setDiag(double *diag, int l_d, SEXP x, int n)
Definition: Mutils.c:461
char La_norm_type(const char *typstr)
Definition: Mutils.c:8
#define LOW
Definition: Mutils.h:82
SEXP dspMatrix_matrix_mm(SEXP a, SEXP b)
Definition: dspMatrix.c:139
SEXP dspMatrix_norm(SEXP obj, SEXP type)
Definition: dspMatrix.c:32
SEXP Matrix_diagSym
Definition: Syms.h:2
SEXP dspMatrix_solve(SEXP a)
Definition: dspMatrix.c:52
SEXP lspMatrix_setDiag(SEXP x, SEXP d)
Definition: dspMatrix.c:114
SEXP symmetricMatrix_validate(SEXP obj)
Definition: dsyMatrix.c:3
double get_norm_sp(SEXP obj, const char *typstr)
Definition: dspMatrix.c:18
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
void d_packed_getDiag(double *dest, SEXP x, int n)
Copy the diagonal elements of the packed denseMatrix x to dest.
Definition: Mutils.c:433
SEXP dspMatrix_matrix_solve(SEXP a, SEXP b)
Definition: dspMatrix.c:69
SEXP l_packed_setDiag(int *diag, int l_d, SEXP x, int n)
Definition: Mutils.c:491
SEXP dspMatrix_as_dsyMatrix(SEXP from)
Definition: dspMatrix.c:121