1 /* R : A Computer Language for Statistical Data Analysis
2 *
3 * Copyright (C) 2012-2013 The R Core Team
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, a copy is available at
17 * https://www.R-project.org/Licenses/.
18 */
19
20 #include <R.h>
21 #include <Rinternals.h>
22 #include <R_ext/Applic.h>
23
24 #include "statsR.h"
25
26 #ifdef ENABLE_NLS
27 #include <libintl.h>
28 #define _(String) dgettext ("stats", String)
29 #else
30 #define _(String) (String)
31 #endif
32
33 /* A wrapper to replace
34
35 z <- .Fortran("dqrls",
36 qr = x, n = n, p = p,
37 y = y, ny = ny,
38 tol = as.double(tol),
39 coefficients = mat.or.vec(p, ny),
40 residuals = y, effects = y, rank = integer(1L),
41 pivot = 1L:p, qraux = double(p), work = double(2*p),
42 PACKAGE="base")
43
44 with less allocation.
45 */
46
47
Cdqrls(SEXP x,SEXP y,SEXP tol,SEXP chk)48 SEXP Cdqrls(SEXP x, SEXP y, SEXP tol, SEXP chk)
49 {
50 SEXP ans;
51 SEXP qr, coefficients, residuals, effects, pivot, qraux;
52 int n, ny = 0, p, rank, nprotect = 4, pivoted = 0;
53 double rtol = asReal(tol), *work;
54 Rboolean check = asLogical(chk);
55
56 ans = getAttrib(x, R_DimSymbol);
57 if(check && length(ans) != 2) error(_("'x' is not a matrix"));
58 int *dims = INTEGER(ans);
59 n = dims[0]; p = dims[1];
60 if(n) ny = (int)(XLENGTH(y)/n); /* y : n x ny, or an n - vector */
61 if(check && n * ny != XLENGTH(y))
62 error(_("dimensions of 'x' (%d,%d) and 'y' (%d) do not match"),
63 n,p, XLENGTH(y));
64
65 /* These lose attributes, so do after we have extracted dims */
66 if (TYPEOF(x) != REALSXP) {
67 PROTECT(x = coerceVector(x, REALSXP));
68 nprotect++;
69 }
70 if (TYPEOF(y) != REALSXP) {
71 PROTECT(y = coerceVector(y, REALSXP));
72 nprotect++;
73 }
74
75 double *rptr = REAL(x);
76 for (R_xlen_t i = 0 ; i < XLENGTH(x) ; i++)
77 if(!R_FINITE(rptr[i])) error(_("NA/NaN/Inf in '%s'"), "x");
78
79 rptr = REAL(y);
80 for (R_xlen_t i = 0 ; i < XLENGTH(y) ; i++)
81 if(!R_FINITE(rptr[i])) error(_("NA/NaN/Inf in '%s'"), "y");
82
83 const char *ansNms[] = {"qr", "coefficients", "residuals", "effects",
84 "rank", "pivot", "qraux", "tol", "pivoted", ""};
85 PROTECT(ans = mkNamed(VECSXP, ansNms));
86 SET_VECTOR_ELT(ans, 0, qr = shallow_duplicate(x));
87 coefficients = (ny > 1) ? allocMatrix(REALSXP, p, ny) : allocVector(REALSXP, p);
88 PROTECT(coefficients);
89 SET_VECTOR_ELT(ans, 1, coefficients);
90 SET_VECTOR_ELT(ans, 2, residuals = shallow_duplicate(y));
91 SET_VECTOR_ELT(ans, 3, effects = shallow_duplicate(y));
92 PROTECT(pivot = allocVector(INTSXP, p));
93 int *ip = INTEGER(pivot);
94 for(int i = 0; i < p; i++) ip[i] = i+1;
95 SET_VECTOR_ELT(ans, 5, pivot);
96 PROTECT(qraux = allocVector(REALSXP, p));
97 SET_VECTOR_ELT(ans, 6, qraux);
98 SET_VECTOR_ELT(ans, 7, tol);
99
100 work = (double *) R_alloc(2 * p, sizeof(double));
101 F77_CALL(dqrls)(REAL(qr), &n, &p, REAL(y), &ny, &rtol,
102 REAL(coefficients), REAL(residuals), REAL(effects),
103 &rank, INTEGER(pivot), REAL(qraux), work);
104 SET_VECTOR_ELT(ans, 4, ScalarInteger(rank));
105 for(int i = 0; i < p; i++)
106 if(ip[i] != i+1) { pivoted = 1; break; }
107 SET_VECTOR_ELT(ans, 8, ScalarLogical(pivoted));
108 UNPROTECT(nprotect);
109
110 return ans;
111 }
112