Matrix r5059
Loading...
Searching...
No Matches
rcond.c
Go to the documentation of this file.
1/* C implementation of methods for rcond */
2
3#include "Lapack-etc.h"
4#include "Mdefines.h"
5
6static
7char La_rcond_type(SEXP s)
8{
9#define ARGNAME "norm"
10 if (TYPEOF(s) != STRSXP)
11 Rf_error(_("argument '%s' is not of type \"%s\""),
12 ARGNAME, "character");
13 if (LENGTH(s) == 0)
14 Rf_error(_("argument '%s' has length %d"),
15 ARGNAME, 0);
16 s = STRING_ELT(s, 0);
17 if (CHAR(s)[0] == '\0' || CHAR(s)[1] != '\0')
18 Rf_error(_("argument '%s' (\"%s\") does not have string length %d"),
19 ARGNAME, CHAR(s), 1);
20 char type = '\0';
21 switch (CHAR(s)[0]) {
22 case 'O':
23 case 'o':
24 case '1':
25 type = 'O';
26 break;
27 case 'I':
28 case 'i':
29 type = 'I';
30 break;
31 default:
32 Rf_error(_("argument '%s' (\"%s\") is not \"%s\", \"%s\", or \"%s\""),
33 ARGNAME, CHAR(s), "O", "1", "I");
34 break;
35 }
36 return type;
37#undef ARGNAME
38}
39
40SEXP geMatrix_rcond(SEXP s_obj, SEXP trf, SEXP s_type)
41{
42 char type = La_rcond_type(s_type);
43
44 int *pdim = DIM(s_obj), n = pdim[1];
45 if (pdim[0] != n)
46 Rf_error(_("non-square matrix arguments not supported by '%s'"),
47 __func__);
48 if (n == 0)
49 return(Rf_ScalarReal(R_PosInf));
50
51 SEXP x = PROTECT(GET_SLOT(s_obj, Matrix_xSym)),
52 y = PROTECT(GET_SLOT(trf, Matrix_xSym));
53 double norm, rcond;
54 int info;
55 double * work = (double *) R_alloc((size_t) n * 4, sizeof(double));
56 if (TYPEOF(x) == CPLXSXP) {
57 double *rwork = (double *) R_alloc((size_t) n * 2, sizeof(double));
58 norm =
59 F77_CALL(zlange)(&type, &n, &n, COMPLEX(x), &n, work FCONE);
60 F77_CALL(zgecon)(&type, &n, COMPLEX(y), &n, &norm, &rcond,
61 (Rcomplex *) work, rwork, &info FCONE);
62 } else {
63 int *iwork = (int *) R_alloc((size_t) n , sizeof(int ));
64 norm =
65 F77_CALL(dlange)(&type, &n, &n, REAL(x), &n, work FCONE);
66 F77_CALL(dgecon)(&type, &n, REAL(y), &n, &norm, &rcond,
67 (double *) work, iwork, &info FCONE);
68 }
69 UNPROTECT(2); /* x, y */
70
71 return Rf_ScalarReal(rcond);
72}
73
74SEXP syMatrix_rcond(SEXP s_obj, SEXP trf, SEXP s_type)
75{
76 char type = La_rcond_type(s_type);
77
78 int n = DIM(s_obj)[0];
79 if (n == 0)
80 return(Rf_ScalarReal(R_PosInf));
81 char ul = UPLO(s_obj);
82
83 SEXP x = PROTECT(GET_SLOT(s_obj, Matrix_xSym)),
84 y = PROTECT(GET_SLOT(trf, Matrix_xSym)),
85 pivot = PROTECT(GET_SLOT(trf, Matrix_permSym));
86 double norm, rcond;
87 int info;
88 if (TYPEOF(x) == CPLXSXP) {
89 double * work = (double *) R_alloc((size_t) n * 4, sizeof(double));
90 if (TRANS(s_obj) == 'C') {
91 norm =
92 F77_CALL(zlanhe)(&type, &ul, &n, COMPLEX(x), &n, work FCONE FCONE);
93 F77_CALL(zhecon)( &ul, &n, COMPLEX(y), &n, INTEGER(pivot), &norm, &rcond,
94 (Rcomplex *) work, &info FCONE);
95 } else {
96 norm =
97 F77_CALL(zlansy)(&type, &ul, &n, COMPLEX(x), &n, work FCONE FCONE);
98 F77_CALL(zsycon)( &ul, &n, COMPLEX(y), &n, INTEGER(pivot), &norm, &rcond,
99 (Rcomplex *) work, &info FCONE);
100 }
101 } else {
102 double * work = (double *) R_alloc((size_t) n * 2, sizeof(double));
103 int *iwork = (int *) R_alloc((size_t) n , sizeof(int ));
104 norm =
105 F77_CALL(dlansy)(&type, &ul, &n, REAL(x), &n, work FCONE FCONE);
106 F77_CALL(dsycon)( &ul, &n, REAL(y), &n, INTEGER(pivot), &norm, &rcond,
107 (double *) work, iwork, &info FCONE);
108 }
109 UNPROTECT(3); /* x, y, pivot */
110
111 return Rf_ScalarReal(rcond);
112}
113
114SEXP spMatrix_rcond(SEXP s_obj, SEXP trf, SEXP s_type)
115{
116 char type = La_rcond_type(s_type);
117
118 int n = DIM(s_obj)[1];
119 if (n == 0)
120 return(Rf_ScalarReal(R_PosInf));
121 char ul = UPLO(s_obj);
122
123 SEXP x = PROTECT(GET_SLOT(s_obj, Matrix_xSym)),
124 y = PROTECT(GET_SLOT(trf, Matrix_xSym)),
125 pivot = PROTECT(GET_SLOT(trf, Matrix_permSym));
126 double norm, rcond;
127 int info;
128 if (TYPEOF(x) == CPLXSXP) {
129 double * work = (double *) R_alloc((size_t) n * 4, sizeof(double));
130 if (TRANS(s_obj) == 'C') {
131 norm =
132 F77_CALL(zlanhp)(&type, &ul, &n, COMPLEX(x), work FCONE FCONE);
133 F77_CALL(zhpcon)( &ul, &n, COMPLEX(y), INTEGER(pivot), &norm, &rcond,
134 (Rcomplex *) work, &info FCONE);
135 } else {
136 norm =
137 F77_CALL(zlansp)(&type, &ul, &n, COMPLEX(x), work FCONE FCONE);
138 F77_CALL(zspcon)( &ul, &n, COMPLEX(y), INTEGER(pivot), &norm, &rcond,
139 (Rcomplex *) work, &info FCONE);
140 }
141 } else {
142 double * work = (double *) R_alloc((size_t) n * 2, sizeof(double));
143 int *iwork = (int *) R_alloc((size_t) n , sizeof(int ));
144 norm =
145 F77_CALL(dlansp)(&type, &ul, &n, REAL(x), work FCONE FCONE);
146 F77_CALL(dspcon)( &ul, &n, REAL(y), INTEGER(pivot), &norm, &rcond,
147 (double *) work, iwork, &info FCONE);
148 }
149 UNPROTECT(3); /* x, y, pivot */
150
151 return Rf_ScalarReal(rcond);
152}
153
154SEXP poMatrix_rcond(SEXP s_obj, SEXP trf, SEXP s_type)
155{
156 char type = La_rcond_type(s_type);
157
158 int n = DIM(s_obj)[1];
159 if (n == 0)
160 return(Rf_ScalarReal(R_PosInf));
161 char ul = UPLO(s_obj);
162
163 SEXP x = PROTECT(GET_SLOT(s_obj, Matrix_xSym)),
164 y = PROTECT(GET_SLOT(trf, Matrix_xSym));
165 double norm, rcond;
166 int info;
167 if (TYPEOF(x) == CPLXSXP) {
168 double * work = (double *) R_alloc((size_t) n * 4, sizeof(double));
169 double *rwork = (double *) R_alloc((size_t) n , sizeof(double));
170 norm =
171 F77_CALL(zlansy)(&type, &ul, &n, COMPLEX(x), &n, work FCONE FCONE);
172 F77_CALL(zpocon)( &ul, &n, COMPLEX(y), &n, &norm, &rcond,
173 (Rcomplex *) work, rwork, &info FCONE);
174 } else {
175 double * work = (double *) R_alloc((size_t) n * 3, sizeof(double));
176 int *iwork = (int *) R_alloc((size_t) n , sizeof(int ));
177 norm =
178 F77_CALL(dlansy)(&type, &ul, &n, REAL(x), &n, work FCONE FCONE);
179 F77_CALL(dpocon)( &ul, &n, REAL(y), &n, &norm, &rcond,
180 (double *) work, iwork, &info FCONE);
181 }
182 UNPROTECT(2); /* x, y */
183
184 return Rf_ScalarReal(rcond);
185}
186
187SEXP ppMatrix_rcond(SEXP s_obj, SEXP trf, SEXP s_type)
188{
189 char type = La_rcond_type(s_type);
190
191 int n = DIM(s_obj)[1];
192 if (n == 0)
193 return(Rf_ScalarReal(R_PosInf));
194 char ul = UPLO(s_obj);
195
196 SEXP x = PROTECT(GET_SLOT(s_obj, Matrix_xSym)),
197 y = PROTECT(GET_SLOT(trf, Matrix_xSym));
198 double norm, rcond;
199 int info;
200 if (TYPEOF(x) == CPLXSXP) {
201 double * work = (double *) R_alloc((size_t) n * 4, sizeof(double));
202 double *rwork = (double *) R_alloc((size_t) n , sizeof(double));
203 norm =
204 F77_CALL(zlansp)(&type, &ul, &n, COMPLEX(x), work FCONE FCONE);
205 F77_CALL(zppcon)( &ul, &n, COMPLEX(y), &norm, &rcond,
206 (Rcomplex *) work, rwork, &info FCONE);
207 } else {
208 double * work = (double *) R_alloc((size_t) n * 4, sizeof(double));
209 int *iwork = (int *) R_alloc((size_t) n , sizeof(int ));
210 norm =
211 F77_CALL(dlansp)(&type, &ul, &n, REAL(x), work FCONE FCONE);
212 F77_CALL(dppcon)( &ul, &n, REAL(y), &norm, &rcond,
213 (double *) work, iwork, &info FCONE);
214 }
215 UNPROTECT(2); /* x, y */
216
217 return Rf_ScalarReal(rcond);
218}
219
220SEXP trMatrix_rcond(SEXP s_obj, SEXP s_type)
221{
222 char type = La_rcond_type(s_type);
223
224 int n = DIM(s_obj)[1];
225 if (n == 0)
226 return(Rf_ScalarReal(R_PosInf));
227 char ul = UPLO(s_obj), nu = DIAG(s_obj);
228
229 SEXP x = PROTECT(GET_SLOT(s_obj, Matrix_xSym));
230 double rcond;
231 int info;
232 if (TYPEOF(x) == CPLXSXP) {
233 double * work = (double *) R_alloc((size_t) n * 4, sizeof(double));
234 double *rwork = (double *) R_alloc((size_t) n , sizeof(double));
235 F77_CALL(ztrcon)(&type, &ul, &nu, &n, COMPLEX(x), &n, &rcond,
236 (Rcomplex *) work, rwork, &info FCONE FCONE FCONE);
237 } else {
238 double * work = (double *) R_alloc((size_t) n * 3, sizeof(double));
239 int *iwork = (int *) R_alloc((size_t) n , sizeof(int ));
240 F77_CALL(dtrcon)(&type, &ul, &nu, &n, REAL(x), &n, &rcond,
241 (double *) work, iwork, &info FCONE FCONE FCONE);
242 }
243 UNPROTECT(1); /* x */
244
245 return Rf_ScalarReal(rcond);
246}
247
248SEXP tpMatrix_rcond(SEXP s_obj, SEXP s_type)
249{
250 char type = La_rcond_type(s_type);
251
252 int n = DIM(s_obj)[1];
253 if (n == 0)
254 return(Rf_ScalarReal(R_PosInf));
255 char ul = UPLO(s_obj), nu = DIAG(s_obj);
256
257 SEXP x = PROTECT(GET_SLOT(s_obj, Matrix_xSym));
258 double rcond;
259 int info;
260 if (TYPEOF(x) == CPLXSXP) {
261 double * work = (double *) R_alloc((size_t) n * 4, sizeof(double));
262 double *rwork = (double *) R_alloc((size_t) n , sizeof(double));
263 F77_CALL(ztpcon)(&type, &ul, &nu, &n, COMPLEX(x), &rcond,
264 (Rcomplex *) work, rwork, &info FCONE FCONE FCONE);
265 } else {
266 double * work = (double *) R_alloc((size_t) n * 3, sizeof(double));
267 int *iwork = (int *) R_alloc((size_t) n , sizeof(int ));
268 F77_CALL(dtpcon)(&type, &ul, &nu, &n, REAL(x), &rcond,
269 (double *) work, iwork, &info FCONE FCONE FCONE);
270 }
271 UNPROTECT(1); /* x */
272
273 return Rf_ScalarReal(rcond);
274}
FCONE
#define FCONE
Definition Lapack-etc.h:13
_
#define _(String)
Definition Mdefines.h:66
DIAG
#define DIAG(x)
Definition Mdefines.h:111
UPLO
#define UPLO(x)
Definition Mdefines.h:101
TRANS
#define TRANS(x)
Definition Mdefines.h:106
DIM
#define DIM(x)
Definition Mdefines.h:85
GET_SLOT
#define GET_SLOT(x, name)
Definition Mdefines.h:72
TYPEOF
#define TYPEOF(s)
Definition Mdefines.h:123
Matrix_permSym
SEXP Matrix_permSym
Definition init.c:623
Matrix_xSym
SEXP Matrix_xSym
Definition init.c:635
trMatrix_rcond
SEXP trMatrix_rcond(SEXP s_obj, SEXP s_type)
Definition rcond.c:220
ppMatrix_rcond
SEXP ppMatrix_rcond(SEXP s_obj, SEXP trf, SEXP s_type)
Definition rcond.c:187
poMatrix_rcond
SEXP poMatrix_rcond(SEXP s_obj, SEXP trf, SEXP s_type)
Definition rcond.c:154
tpMatrix_rcond
SEXP tpMatrix_rcond(SEXP s_obj, SEXP s_type)
Definition rcond.c:248
geMatrix_rcond
SEXP geMatrix_rcond(SEXP s_obj, SEXP trf, SEXP s_type)
Definition rcond.c:40
ARGNAME
#define ARGNAME
syMatrix_rcond
SEXP syMatrix_rcond(SEXP s_obj, SEXP trf, SEXP s_type)
Definition rcond.c:74
La_rcond_type
static char La_rcond_type(SEXP s)
Definition rcond.c:7
spMatrix_rcond
SEXP spMatrix_rcond(SEXP s_obj, SEXP trf, SEXP s_type)
Definition rcond.c:114
Generated by doxygen 1.12.0