Matrix r4655
Loading...
Searching...
No Matches
kappa.c
Go to the documentation of this file.
1#include "Lapack-etc.h"
2#include "Mdefines.h"
3#include "kappa.h"
4
5static
6char La_norm_type(SEXP s)
7{
8#define ARGNAME "type"
9 if (TYPEOF(s) != STRSXP)
10 error(_("argument '%s' is not of type \"%s\""),
11 ARGNAME, "character");
12 if (LENGTH(s) == 0)
13 error(_("argument '%s' has length %d"),
14 ARGNAME, 0);
15 const char *type = CHAR(STRING_ELT(s, 0));
16 if (type[0] == '\0' || type[1] != '\0')
17 error(_("argument '%s' (\"%s\") does not have string length %d"),
18 ARGNAME, type, 1);
19 char t = '\0';
20 switch (type[0]) {
21 case 'M':
22 case 'm':
23 t = 'M';
24 break;
25 case 'O':
26 case 'o':
27 case '1':
28 t = 'O';
29 break;
30 case 'I':
31 case 'i':
32 t = 'I';
33 break;
34 case 'F':
35 case 'f':
36 case 'E':
37 case 'e':
38 t = 'F';
39 break;
40 default:
41 error(_("argument '%s' (\"%s\") is not \"%s\", \"%s\", \"%s\", \"%s\", \"%s\", or \"%s\""),
42 ARGNAME, type, "M", "O", "1", "I", "F", "E");
43 break;
44 }
45 return t;
46#undef ARGNAME
47}
48
49static
50char La_rcond_type(SEXP s)
51{
52#define ARGNAME "norm"
53 if (TYPEOF(s) != STRSXP)
54 error(_("argument '%s' is not of type \"%s\""),
55 ARGNAME, "character");
56 if (LENGTH(s) == 0)
57 error(_("argument '%s' has length %d"),
58 ARGNAME, 0);
59 const char *type = CHAR(STRING_ELT(s, 0));
60 if (type[0] == '\0' || type[1] != '\0')
61 error(_("argument '%s' (\"%s\") does not have string length %d"),
62 ARGNAME, type, 1);
63 char t = '\0';
64 switch (type[0]) {
65 case 'O':
66 case 'o':
67 case '1':
68 t = 'O';
69 break;
70 case 'I':
71 case 'i':
72 t = 'I';
73 break;
74 default:
75 error(_("argument '%s' (\"%s\") is not \"%s\", \"%s\", or \"%s\""),
76 ARGNAME, type, "O", "1", "I");
77 break;
78 }
79 return t;
80#undef ARGNAME
81}
82
83SEXP dgeMatrix_norm(SEXP obj, SEXP type)
84{
85 char t = La_norm_type(type);
86
87 SEXP dim = GET_SLOT(obj, Matrix_DimSym);
88 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1];
89 if (m == 0 || n == 0)
90 return ScalarReal(0.0);
91
92 SEXP x = PROTECT(GET_SLOT(obj, Matrix_xSym));
93 double norm, *work = NULL;
94 if (t == 'I')
95 work = (double *) R_alloc((size_t) m, sizeof(double));
96#ifdef MATRIX_ENABLE_ZMATRIX
97 if (TYPEOF(x) == CPLXSXP)
98 norm =
99 F77_CALL(zlange)(&t, &m, &n, COMPLEX(x), &m, work FCONE);
100 else
101#endif
102 norm =
103 F77_CALL(dlange)(&t, &m, &n, REAL(x), &m, work FCONE);
104 UNPROTECT(1); /* x */
105
106 return ScalarReal(norm);
107}
108
109SEXP dsyMatrix_norm(SEXP obj, SEXP type)
110{
111 char t = La_norm_type(type);
112
113 SEXP dim = GET_SLOT(obj, Matrix_DimSym);
114 int n = INTEGER(dim)[1];
115 if (n == 0)
116 return ScalarReal(0.0);
117
118 SEXP uplo = GET_SLOT(obj, Matrix_uploSym);
119 char ul = *CHAR(STRING_ELT(uplo, 0));
120
121 SEXP x = PROTECT(GET_SLOT(obj, Matrix_xSym));
122 double norm, *work = NULL;
123 if (t == 'O' || t == 'I')
124 work = (double *) R_alloc((size_t) n, sizeof(double));
125#ifdef MATRIX_ENABLE_ZMATRIX
126 if (TYPEOF(x) == CPLXSXP)
127 norm =
128 F77_CALL(zlansy)(&t, &ul, &n, COMPLEX(x), &n, work FCONE FCONE);
129 else
130#endif
131 norm =
132 F77_CALL(dlansy)(&t, &ul, &n, REAL(x), &n, work FCONE FCONE);
133 UNPROTECT(1); /* x */
134
135 return ScalarReal(norm);
136}
137
138SEXP dspMatrix_norm(SEXP obj, SEXP type)
139{
140 char t = La_norm_type(type);
141
142 SEXP dim = GET_SLOT(obj, Matrix_DimSym);
143 int n = INTEGER(dim)[1];
144 if (n == 0)
145 return ScalarReal(0.0);
146
147 SEXP uplo = GET_SLOT(obj, Matrix_uploSym);
148 char ul = CHAR(STRING_ELT(uplo, 0))[0];
149
150 SEXP x = PROTECT(GET_SLOT(obj, Matrix_xSym));
151 double norm, *work = NULL;
152 if (t == 'O' || t == 'I')
153 work = (double *) R_alloc((size_t) n, sizeof(double));
154#ifdef MATRIX_ENABLE_ZMATRIX
155 if (TYPEOF(x) == CPLXSXP)
156 norm =
157 F77_CALL(zlansp)(&t, &ul, &n, COMPLEX(x), work FCONE FCONE);
158 else
159#endif
160 norm =
161 F77_CALL(dlansp)(&t, &ul, &n, REAL(x), work FCONE FCONE);
162 UNPROTECT(1); /* x */
163
164 return ScalarReal(norm);
165}
166
167SEXP dtrMatrix_norm(SEXP obj, SEXP type)
168{
169 char t = La_norm_type(type);
170
171 SEXP dim = GET_SLOT(obj, Matrix_DimSym);
172 int n = INTEGER(dim)[0];
173 if (n == 0)
174 return ScalarReal(0.0);
175
176 SEXP uplo = GET_SLOT(obj, Matrix_uploSym);
177 char ul = CHAR(STRING_ELT(uplo, 0))[0];
178
179 SEXP diag = GET_SLOT(obj, Matrix_diagSym);
180 char di = CHAR(STRING_ELT(diag, 0))[0];
181
182 SEXP x = PROTECT(GET_SLOT(obj, Matrix_xSym));
183 double norm, *work = NULL;
184 if (t == 'I')
185 work = (double *) R_alloc((size_t) n, sizeof(double));
186#ifdef MATRIX_ENABLE_ZMATRIX
187 if (TYPEOF(x) == CPLXSXP)
188 norm =
189 F77_CALL(zlantr)(&t, &ul, &di, &n, &n, COMPLEX(x), &n, work FCONE FCONE FCONE);
190 else
191#endif
192 norm =
193 F77_CALL(dlantr)(&t, &ul, &di, &n, &n, REAL(x), &n, work FCONE FCONE FCONE);
194 UNPROTECT(1); /* x */
195
196 return ScalarReal(norm);
197}
198
199SEXP dtpMatrix_norm(SEXP obj, SEXP type)
200{
201 char t = La_norm_type(type);
202
203 SEXP dim = GET_SLOT(obj, Matrix_DimSym);
204 int n = INTEGER(dim)[0];
205 if (n == 0)
206 return ScalarReal(0.0);
207
208 SEXP uplo = GET_SLOT(obj, Matrix_uploSym);
209 char ul = CHAR(STRING_ELT(uplo, 0))[0];
210
211 SEXP diag = GET_SLOT(obj, Matrix_diagSym);
212 char di = CHAR(STRING_ELT(diag, 0))[0];
213
214 SEXP x = PROTECT(GET_SLOT(obj, Matrix_xSym));
215 double norm, *work = NULL;
216 if (t == 'I')
217 work = (double *) R_alloc((size_t) n, sizeof(double));
218#ifdef MATRIX_ENABLE_ZMATRIX
219 if (TYPEOF(x) == CPLXSXP)
220 norm =
221 F77_CALL(zlantp)(&t, &ul, &di, &n, COMPLEX(x), work FCONE FCONE FCONE);
222 else
223#endif
224 norm =
225 F77_CALL(dlantp)(&t, &ul, &di, &n, REAL(x), work FCONE FCONE FCONE);
226 UNPROTECT(1); /* x */
227
228 return ScalarReal(norm);
229}
230
231SEXP dgeMatrix_rcond(SEXP obj, SEXP trf, SEXP type)
232{
233 char t = La_rcond_type(type);
234
235 SEXP dim = GET_SLOT(obj, Matrix_DimSym);
236 int *pdim = INTEGER(dim), m = pdim[0], n = pdim[1];
237 if (m != n)
238 error(_("%s(%s) is undefined: '%s' is not square"), "rcond", "x", "x");
239 if (n == 0)
240 return(ScalarReal(R_PosInf));
241
242 SEXP x = PROTECT(GET_SLOT(obj, Matrix_xSym)),
243 y = PROTECT(GET_SLOT(trf, Matrix_xSym));
244 double norm, rcond;
245 int info;
246 double * work = (double *) R_alloc((size_t) 4 * n, sizeof(double));
247#ifdef MATRIX_ENABLE_ZMATRIX
248 if (TYPEOF(x) == CPLXSXP) {
249 double *rwork = (double *) R_alloc((size_t) 2 * n, sizeof(double));
250 norm =
251 F77_CALL(zlange)(&t, &n, &n, COMPLEX(x), &n, work FCONE);
252 F77_CALL(zgecon)(&t, &n, COMPLEX(y), &n, &norm, &rcond,
253 (Rcomplex *) work, rwork, &info FCONE);
254 } else {
255#endif
256 int *iwork = (int *) R_alloc((size_t) n, sizeof(int ));
257 norm =
258 F77_CALL(dlange)(&t, &n, &n, REAL(x), &n, work FCONE);
259 F77_CALL(dgecon)(&t, &n, REAL(y), &n, &norm, &rcond,
260 (double *) work, iwork, &info FCONE);
261#ifdef MATRIX_ENABLE_ZMATRIX
262 }
263#endif
264 UNPROTECT(2); /* x, y */
265
266 return ScalarReal(rcond);
267}
268
269SEXP dsyMatrix_rcond(SEXP obj, SEXP trf, SEXP type)
270{
271 char t = La_rcond_type(type);
272
273 SEXP dim = GET_SLOT(obj, Matrix_DimSym);
274 int n = INTEGER(dim)[0];
275 if (n == 0)
276 return(ScalarReal(R_PosInf));
277
278 SEXP uplo = GET_SLOT(obj, Matrix_uploSym);
279 char ul = CHAR(STRING_ELT(uplo, 0))[0];
280
281 SEXP x = PROTECT(GET_SLOT(obj, Matrix_xSym)),
282 y = PROTECT(GET_SLOT(trf, Matrix_xSym)),
283 pivot = PROTECT(GET_SLOT(trf, Matrix_permSym));
284 double norm, rcond;
285 int info;
286#ifdef MATRIX_ENABLE_ZMATRIX
287 if (TYPEOF(x) == CPLXSXP) {
288 double * work = (double *) R_alloc((size_t) 4 * n, sizeof(double));
289 norm =
290 F77_CALL(zlansy)(&t, &ul, &n, COMPLEX(x), &n, work FCONE FCONE);
291 F77_CALL(zsycon)( &ul, &n, COMPLEX(y), &n, INTEGER(pivot), &norm, &rcond,
292 (Rcomplex *) work, &info FCONE);
293 } else {
294#endif
295 double * work = (double *) R_alloc((size_t) 2 * n, sizeof(double));
296 int *iwork = (int *) R_alloc((size_t) n, sizeof(int ));
297 norm =
298 F77_CALL(dlansy)(&t, &ul, &n, REAL(x), &n, work FCONE FCONE);
299 F77_CALL(dsycon)( &ul, &n, REAL(y), &n, INTEGER(pivot), &norm, &rcond,
300 (double *) work, iwork, &info FCONE);
301#ifdef MATRIX_ENABLE_ZMATRIX
302 }
303#endif
304 UNPROTECT(3); /* x, y, pivot */
305
306 return ScalarReal(rcond);
307}
308
309SEXP dspMatrix_rcond(SEXP obj, SEXP trf, SEXP type)
310{
311 char t = La_rcond_type(type);
312
313 SEXP dim = GET_SLOT(obj, Matrix_DimSym);
314 int n = INTEGER(dim)[0];
315 if (n == 0)
316 return(ScalarReal(R_PosInf));
317
318 SEXP uplo = GET_SLOT(obj, Matrix_uploSym);
319 char ul = CHAR(STRING_ELT(uplo, 0))[0];
320
321 SEXP x = PROTECT(GET_SLOT(obj, Matrix_xSym)),
322 y = PROTECT(GET_SLOT(trf, Matrix_xSym)),
323 pivot = PROTECT(GET_SLOT(trf, Matrix_permSym));
324 double norm, rcond;
325 int info;
326#ifdef MATRIX_ENABLE_ZMATRIX
327 if (TYPEOF(x) == CPLXSXP) {
328 double * work = (double *) R_alloc((size_t) 4 * n, sizeof(double));
329 norm =
330 F77_CALL(zlansp)(&t, &ul, &n, COMPLEX(x), work FCONE FCONE);
331 F77_CALL(zspcon)( &ul, &n, COMPLEX(y), INTEGER(pivot), &norm, &rcond,
332 (Rcomplex *) work, &info FCONE);
333 } else {
334#endif
335 double * work = (double *) R_alloc((size_t) 2 * n, sizeof(double));
336 int *iwork = (int *) R_alloc((size_t) n, sizeof(int ));
337 norm =
338 F77_CALL(dlansp)(&t, &ul, &n, REAL(x), work FCONE FCONE);
339 F77_CALL(dspcon)( &ul, &n, REAL(y), INTEGER(pivot), &norm, &rcond,
340 (double *) work, iwork, &info FCONE);
341#ifdef MATRIX_ENABLE_ZMATRIX
342 }
343#endif
344 UNPROTECT(3); /* x, y, pivot */
345
346 return ScalarReal(rcond);
347}
348
349SEXP dpoMatrix_rcond(SEXP obj, SEXP trf, SEXP type)
350{
351 char t = La_rcond_type(type);
352
353 SEXP dim = GET_SLOT(obj, Matrix_DimSym);
354 int n = INTEGER(dim)[0];
355 if (n == 0)
356 return(ScalarReal(R_PosInf));
357
358 SEXP uplo = GET_SLOT(obj, Matrix_uploSym);
359 char ul = CHAR(STRING_ELT(uplo, 0))[0];
360
361 SEXP x = PROTECT(GET_SLOT(obj, Matrix_xSym)),
362 y = PROTECT(GET_SLOT(trf, Matrix_xSym));
363 double norm, rcond;
364 int info;
365#ifdef MATRIX_ENABLE_ZMATRIX
366 if (TYPEOF(x) == CPLXSXP) {
367 double * work = (double *) R_alloc((size_t) 4 * n, sizeof(double));
368 double *rwork = (double *) R_alloc((size_t) n, sizeof(double));
369 norm =
370 F77_CALL(zlansy)(&t, &ul, &n, COMPLEX(x), &n, work FCONE FCONE);
371 F77_CALL(zpocon)( &ul, &n, COMPLEX(y), &n, &norm, &rcond,
372 (Rcomplex *) work, rwork, &info FCONE);
373 } else {
374#endif
375 double * work = (double *) R_alloc((size_t) 3 * n, sizeof(double));
376 int *iwork = (int *) R_alloc((size_t) n, sizeof(int ));
377 norm =
378 F77_CALL(dlansy)(&t, &ul, &n, REAL(x), &n, work FCONE FCONE);
379 F77_CALL(dpocon)( &ul, &n, REAL(y), &n, &norm, &rcond,
380 (double *) work, iwork, &info FCONE);
381#ifdef MATRIX_ENABLE_ZMATRIX
382 }
383#endif
384 UNPROTECT(2); /* x, y */
385
386 return ScalarReal(rcond);
387}
388
389SEXP dppMatrix_rcond(SEXP obj, SEXP trf, SEXP type)
390{
391 char t = La_rcond_type(type);
392
393 SEXP dim = GET_SLOT(obj, Matrix_DimSym);
394 int n = INTEGER(dim)[0];
395 if (n == 0)
396 return(ScalarReal(R_PosInf));
397
398 SEXP uplo = GET_SLOT(obj, Matrix_uploSym);
399 char ul = CHAR(STRING_ELT(uplo, 0))[0];
400
401 SEXP x = PROTECT(GET_SLOT(obj, Matrix_xSym)),
402 y = PROTECT(GET_SLOT(trf, Matrix_xSym));
403 double norm, rcond;
404 int info;
405#ifdef MATRIX_ENABLE_ZMATRIX
406 if (TYPEOF(x) == CPLXSXP) {
407 double * work = (double *) R_alloc((size_t) 4 * n, sizeof(double));
408 double *rwork = (double *) R_alloc((size_t) n, sizeof(double));
409 norm =
410 F77_CALL(zlansp)(&t, &ul, &n, COMPLEX(x), work FCONE FCONE);
411 F77_CALL(zppcon)( &ul, &n, COMPLEX(y), &norm, &rcond,
412 (Rcomplex *) work, rwork, &info FCONE);
413 } else {
414#endif
415 double * work = (double *) R_alloc((size_t) 3 * n, sizeof(double));
416 int *iwork = (int *) R_alloc((size_t) n, sizeof(int ));
417 norm =
418 F77_CALL(dlansp)(&t, &ul, &n, REAL(x), work FCONE FCONE);
419 F77_CALL(dppcon)( &ul, &n, REAL(y), &norm, &rcond,
420 (double *) work, iwork, &info FCONE);
421#ifdef MATRIX_ENABLE_ZMATRIX
422 }
423#endif
424 UNPROTECT(2); /* x, y */
425
426 return ScalarReal(rcond);
427}
428
429SEXP dtrMatrix_rcond(SEXP obj, SEXP type)
430{
431 char t = La_rcond_type(type);
432
433 SEXP dim = GET_SLOT(obj, Matrix_DimSym);
434 int n = INTEGER(dim)[0];
435 if (n == 0)
436 return(ScalarReal(R_PosInf));
437
438 SEXP uplo = GET_SLOT(obj, Matrix_uploSym);
439 char ul = CHAR(STRING_ELT(uplo, 0))[0];
440
441 SEXP diag = GET_SLOT(obj, Matrix_diagSym);
442 char di = CHAR(STRING_ELT(diag, 0))[0];
443
444 SEXP x = PROTECT(GET_SLOT(obj, Matrix_xSym));
445 double rcond;
446 int info;
447#ifdef MATRIX_ENABLE_ZMATRIX
448 if (TYPEOF(x) == CPLXSXP) {
449 double * work = (double *) R_alloc((size_t) 4 * n, sizeof(double));
450 double *rwork = (double *) R_alloc((size_t) n, sizeof(double));
451 F77_CALL(ztrcon)(&t, &ul, &di, &n, COMPLEX(x), &n, &rcond,
452 (Rcomplex *) work, rwork, &info FCONE FCONE FCONE);
453 } else {
454#endif
455 double * work = (double *) R_alloc((size_t) 3 * n, sizeof(double));
456 int *iwork = (int *) R_alloc((size_t) n, sizeof(int ));
457 F77_CALL(dtrcon)(&t, &ul, &di, &n, REAL(x), &n, &rcond,
458 (double *) work, iwork, &info FCONE FCONE FCONE);
459#ifdef MATRIX_ENABLE_ZMATRIX
460 }
461#endif
462 UNPROTECT(1); /* x */
463
464 return ScalarReal(rcond);
465}
466
467SEXP dtpMatrix_rcond(SEXP obj, SEXP type)
468{
469 char t = La_rcond_type(type);
470
471 SEXP dim = GET_SLOT(obj, Matrix_DimSym);
472 int n = INTEGER(dim)[0];
473 if (n == 0)
474 return(ScalarReal(R_PosInf));
475
476 SEXP uplo = GET_SLOT(obj, Matrix_uploSym);
477 char ul = CHAR(STRING_ELT(uplo, 0))[0];
478
479 SEXP diag = GET_SLOT(obj, Matrix_diagSym);
480 char di = CHAR(STRING_ELT(diag, 0))[0];
481
482 SEXP x = PROTECT(GET_SLOT(obj, Matrix_xSym));
483 double rcond;
484 int info;
485#ifdef MATRIX_ENABLE_ZMATRIX
486 if (TYPEOF(x) == CPLXSXP) {
487 double * work = (double *) R_alloc((size_t) 4 * n, sizeof(double));
488 double *rwork = (double *) R_alloc((size_t) n, sizeof(double));
489 F77_CALL(ztpcon)(&t, &ul, &di, &n, COMPLEX(x), &rcond,
490 (Rcomplex *) work, rwork, &info FCONE FCONE FCONE);
491 } else {
492#endif
493 double * work = (double *) R_alloc((size_t) 3 * n, sizeof(double));
494 int *iwork = (int *) R_alloc((size_t) n, sizeof(int ));
495 F77_CALL(dtpcon)(&t, &ul, &di, &n, REAL(x), &rcond,
496 (double *) work, iwork, &info FCONE FCONE FCONE);
497#ifdef MATRIX_ENABLE_ZMATRIX
498 }
499#endif
500 UNPROTECT(1); /* x */
501
502 return ScalarReal(rcond);
503}
FCONE
#define FCONE
Definition Lapack-etc.h:18
_
#define _(String)
Definition Mdefines.h:44
GET_SLOT
#define GET_SLOT(x, what)
Definition Mdefines.h:85
Matrix_permSym
SEXP Matrix_permSym
Definition Msymbols.h:18
Matrix_DimSym
SEXP Matrix_DimSym
Definition Msymbols.h:3
Matrix_xSym
SEXP Matrix_xSym
Definition Msymbols.h:22
Matrix_diagSym
SEXP Matrix_diagSym
Definition Msymbols.h:11
Matrix_uploSym
SEXP Matrix_uploSym
Definition Msymbols.h:21
dsyMatrix_rcond
SEXP dsyMatrix_rcond(SEXP obj, SEXP trf, SEXP type)
Definition kappa.c:269
dspMatrix_norm
SEXP dspMatrix_norm(SEXP obj, SEXP type)
Definition kappa.c:138
dtrMatrix_rcond
SEXP dtrMatrix_rcond(SEXP obj, SEXP type)
Definition kappa.c:429
La_norm_type
static char La_norm_type(SEXP s)
Definition kappa.c:6
dtpMatrix_rcond
SEXP dtpMatrix_rcond(SEXP obj, SEXP type)
Definition kappa.c:467
dpoMatrix_rcond
SEXP dpoMatrix_rcond(SEXP obj, SEXP trf, SEXP type)
Definition kappa.c:349
dgeMatrix_norm
SEXP dgeMatrix_norm(SEXP obj, SEXP type)
Definition kappa.c:83
dppMatrix_rcond
SEXP dppMatrix_rcond(SEXP obj, SEXP trf, SEXP type)
Definition kappa.c:389
dsyMatrix_norm
SEXP dsyMatrix_norm(SEXP obj, SEXP type)
Definition kappa.c:109
dtpMatrix_norm
SEXP dtpMatrix_norm(SEXP obj, SEXP type)
Definition kappa.c:199
dgeMatrix_rcond
SEXP dgeMatrix_rcond(SEXP obj, SEXP trf, SEXP type)
Definition kappa.c:231
ARGNAME
#define ARGNAME
dtrMatrix_norm
SEXP dtrMatrix_norm(SEXP obj, SEXP type)
Definition kappa.c:167
La_rcond_type
static char La_rcond_type(SEXP s)
Definition kappa.c:50
dspMatrix_rcond
SEXP dspMatrix_rcond(SEXP obj, SEXP trf, SEXP type)
Definition kappa.c:309
Generated by doxygen 1.9.8