| /* |
| * R : A Computer Language for Statistical Data Analysis |
| * Copyright (C) 2001--2019 The R Core Team. |
| * Copyright (C) 2003--2010 The R Foundation |
| * |
| * This program is free software; you can redistribute it and/or modify |
| * it under the terms of the GNU General Public License as published by |
| * the Free Software Foundation; either version 2 of the License, or |
| * (at your option) any later version. |
| * |
| * This program is distributed in the hope that it will be useful, |
| * but WITHOUT ANY WARRANTY; without even the implied warranty of |
| * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| * GNU General Public License for more details. |
| * |
| * You should have received a copy of the GNU General Public License |
| * along with this program; if not, a copy is available at |
| * https://www.R-project.org/Licenses/ |
| */ |
| |
| /* Interface routines, callable from R using .Internal, for Lapack code */ |
| |
| #ifdef HAVE_CONFIG_H |
| # include <config.h> |
| #endif |
| |
| #include <Defn.h> |
| #include <ctype.h> /* for toupper */ |
| #include <limits.h> /* for PATH_MAX */ |
| #include <stdlib.h> /* for realpath */ |
| |
| #ifdef HAVE_UNISTD_H |
| # include <unistd.h> /* for realpath on some systems */ |
| #endif |
| |
| #ifdef HAVE_DLFCN_H |
| # include <dlfcn.h> /* for dladdr */ |
| #endif |
| |
| #if defined(HAVE_REALPATH) && defined(HAVE_DECL_REALPATH) && !HAVE_DECL_REALPATH |
| extern char *realpath(const char *path, char *resolved_path); |
| #endif |
| |
| #if defined(HAVE_DLADDR) && defined(HAVE_DECL_DLADDR) && !HAVE_DECL_DLADDR |
| extern int dladdr(void *addr, Dl_info *info); |
| #endif |
| |
| #include "Lapack.h" |
| |
| /* NB: the handling of dims is odd here. Most are coerced to be |
| * integers (which dimgets currently guarantees), but a couple were |
| * used unchecked. |
| */ |
| |
| /* FIXME: MM would want to make these available to packages; |
| * BUT: 1) cannot be in libRlapack {since that may be external} |
| * 2) Pkgs cannot get it from the C-Lapack interface code {lapack.so} |
| * since that is R-internal |
| */ |
| static char La_norm_type(const char *typstr) |
| { |
| char typup; |
| |
| if (strlen(typstr) != 1) |
| error( |
| _("argument type[1]='%s' must be a character string of string length 1"), |
| typstr); |
| typup = (char) toupper(*typstr); |
| if (typup == '1') |
| typup = 'O'; /* aliases */ |
| else if (typup == 'E') |
| typup = 'F'; |
| else if (typup != 'M' && typup != 'O' && typup != 'I' && typup != 'F') |
| error(_("argument type[1]='%s' must be one of 'M','1','O','I','F' or 'E'"), |
| typstr); |
| return typup; |
| } |
| |
| /* Lapack condition number approximation: currently only supports _1 or _Inf norm : */ |
| static char La_rcond_type(const char *typstr) |
| { |
| char typup; |
| |
| if (strlen(typstr) != 1) |
| error(_("argument type[1]='%s' must be a character string of string length 1"), |
| typstr); |
| typup = (char) toupper(*typstr); |
| if (typup == '1') |
| typup = 'O'; /* alias */ |
| else if (typup != 'O' && typup != 'I') |
| error(_("argument type[1]='%s' must be one of '1','O', or 'I'"), |
| typstr); |
| return typup; /* 'O' or 'I' */ |
| } |
| |
| /* La.svd, called from svd */ |
| static SEXP La_svd(SEXP jobu, SEXP x, SEXP s, SEXP u, SEXP vt) |
| { |
| if (!isString(jobu)) |
| error("'jobu' must be a character string"); |
| int *xdims = INTEGER(coerceVector(getAttrib(x, R_DimSymbol), INTSXP)), |
| n = xdims[0], p = xdims[1], nprot = 2; |
| |
| /* work on a copy of x */ |
| double *xvals; |
| if (!isReal(x)) { |
| x = PROTECT(coerceVector(x, REALSXP)); nprot++; |
| xvals = REAL(x); |
| } else { |
| xvals = (double *) R_alloc(n * (size_t) p, sizeof(double)); |
| Memcpy(xvals, REAL(x), n * (size_t) p); |
| } |
| |
| SEXP dims = getAttrib(u, R_DimSymbol); |
| if (TYPEOF(dims) != INTSXP) error("non-integer dim(u)"); |
| int ldu = INTEGER(dims)[0]; |
| dims = getAttrib(vt, R_DimSymbol); |
| if (TYPEOF(dims) != INTSXP) error("non-integer dim(vt)"); |
| int ldvt = INTEGER(dims)[0]; |
| double tmp; |
| /* min(n,p) large is implausible, but cast to be sure */ |
| int *iwork= (int *) R_alloc(8*(size_t)(n < p ? n : p), sizeof(int)); |
| |
| /* ask for optimal size of work array */ |
| const char *ju = CHAR(STRING_ELT(jobu, 0)); |
| int info = 0, lwork = -1; |
| F77_CALL(dgesdd)(ju, &n, &p, xvals, &n, REAL(s), |
| REAL(u), &ldu, REAL(vt), &ldvt, |
| &tmp, &lwork, iwork, &info FCONE); |
| if (info != 0) |
| error(_("error code %d from Lapack routine '%s'"), info, "dgesdd"); |
| lwork = (int) tmp; |
| double *work = (double *) R_alloc(lwork, sizeof(double)); |
| F77_CALL(dgesdd)(ju, &n, &p, xvals, &n, REAL(s), |
| REAL(u), &ldu, REAL(vt), &ldvt, |
| work, &lwork, iwork, &info FCONE); |
| if (info != 0) |
| error(_("error code %d from Lapack routine '%s'"), info, "dgesdd"); |
| |
| SEXP val = PROTECT(allocVector(VECSXP, 3)); |
| SEXP nm = PROTECT(allocVector(STRSXP, 3)); |
| SET_STRING_ELT(nm, 0, mkChar("d")); |
| SET_STRING_ELT(nm, 1, mkChar("u")); |
| SET_STRING_ELT(nm, 2, mkChar("vt")); |
| setAttrib(val, R_NamesSymbol, nm); |
| SET_VECTOR_ELT(val, 0, s); |
| SET_VECTOR_ELT(val, 1, u); |
| SET_VECTOR_ELT(val, 2, vt); |
| UNPROTECT(nprot); |
| return val; |
| } |
| |
| /* Real, symmetric case of eigen */ |
| static SEXP La_rs(SEXP x, SEXP only_values) |
| { |
| int *xdims, n, lwork, info = 0, ov; |
| char jobv[2] = "U", uplo[2] = "L", range[2] = "A"; |
| SEXP z = R_NilValue; |
| double *work, *rx, *rvalues, tmp, *rz = NULL; |
| int liwork, *iwork, itmp, m; |
| double vl = 0.0, vu = 0.0, abstol = 0.0; |
| /* valgrind seems to think vu should be set, but it is documented |
| not to be used if range='a' */ |
| int il, iu, *isuppz; |
| |
| xdims = INTEGER(coerceVector(getAttrib(x, R_DimSymbol), INTSXP)); |
| n = xdims[0]; |
| if (n != xdims[1]) error(_("'x' must be a square numeric matrix")); |
| ov = asLogical(only_values); |
| if (ov == NA_LOGICAL) error(_("invalid '%s' argument"), "only.values"); |
| if (ov) jobv[0] = 'N'; else jobv[0] = 'V'; |
| |
| /* work on a copy of x, since LAPACK trashes it */ |
| if (!isReal(x)) { |
| x = coerceVector(x, REALSXP); |
| rx = REAL(x); |
| } else { |
| rx = (double *) R_alloc(n * (size_t) n, sizeof(double)); |
| Memcpy(rx, REAL(x), (size_t) n * n); |
| } |
| PROTECT(x); |
| SEXP values = PROTECT(allocVector(REALSXP, n)); |
| rvalues = REAL(values); |
| |
| if (!ov) { |
| z = PROTECT(allocMatrix(REALSXP, n, n)); |
| rz = REAL(z); |
| } |
| isuppz = (int *) R_alloc(2*(size_t)n, sizeof(int)); |
| /* ask for optimal size of work arrays */ |
| lwork = -1; liwork = -1; |
| F77_CALL(dsyevr)(jobv, range, uplo, &n, rx, &n, |
| &vl, &vu, &il, &iu, &abstol, &m, rvalues, |
| rz, &n, isuppz, |
| &tmp, &lwork, &itmp, &liwork, &info FCONE FCONE FCONE); |
| if (info != 0) |
| error(_("error code %d from Lapack routine '%s'"), info, "dsyevr"); |
| lwork = (int) tmp; |
| liwork = itmp; |
| |
| work = (double *) R_alloc(lwork, sizeof(double)); |
| iwork = (int *) R_alloc(liwork, sizeof(int)); |
| F77_CALL(dsyevr)(jobv, range, uplo, &n, rx, &n, |
| &vl, &vu, &il, &iu, &abstol, &m, rvalues, |
| rz, &n, isuppz, |
| work, &lwork, iwork, &liwork, &info FCONE FCONE FCONE); |
| if (info != 0) |
| error(_("error code %d from Lapack routine '%s'"), info, "dsyevr"); |
| |
| SEXP ret, nm; |
| if (!ov) { |
| ret = PROTECT(allocVector(VECSXP, 2)); |
| nm = PROTECT(allocVector(STRSXP, 2)); |
| SET_STRING_ELT(nm, 1, mkChar("vectors")); |
| SET_VECTOR_ELT(ret, 1, z); |
| } else { |
| ret = PROTECT(allocVector(VECSXP, 1)); |
| nm = PROTECT(allocVector(STRSXP, 1)); |
| } |
| SET_STRING_ELT(nm, 0, mkChar("values")); |
| setAttrib(ret, R_NamesSymbol, nm); |
| SET_VECTOR_ELT(ret, 0, values); |
| UNPROTECT(ov ? 4 : 5); |
| return ret; |
| } |
| |
| static SEXP unscramble(const double* imaginary, int n, |
| const double* vecs) |
| { |
| int i, j; |
| SEXP s = allocMatrix(CPLXSXP, n, n); |
| size_t N = n; |
| for (j = 0; j < n; j++) { |
| if (imaginary[j] != 0) { |
| int j1 = j + 1; |
| for (i = 0; i < n; i++) { |
| COMPLEX(s)[i+N*j].r = COMPLEX(s)[i+N*j1].r = vecs[i + j * N]; |
| COMPLEX(s)[i+N*j1].i = -(COMPLEX(s)[i+N*j].i = vecs[i + j1 * N]); |
| } |
| j = j1; |
| } else { |
| for (i = 0; i < n; i++) { |
| COMPLEX(s)[i+N*j].r = vecs[i + j * N]; |
| COMPLEX(s)[i+N*j].i = 0.0; |
| } |
| } |
| } |
| return s; |
| } |
| |
| /* Real, general case of eigen */ |
| static SEXP La_rg(SEXP x, SEXP only_values) |
| { |
| Rboolean vectors, complexValues; |
| int i, n, lwork, info, *xdims, ov; |
| double *work, *wR, *wI, *left, *right, *xvals, tmp; |
| char jobVL[2] = "N", jobVR[2] = "N"; |
| |
| xdims = INTEGER(coerceVector(getAttrib(x, R_DimSymbol), INTSXP)); |
| n = xdims[0]; |
| if (n != xdims[1]) |
| error(_("'x' must be a square numeric matrix")); |
| |
| /* work on a copy of x */ |
| if (!isReal(x)) { |
| x = coerceVector(x, REALSXP); |
| xvals = REAL(x); |
| } else { |
| xvals = (double *) R_alloc(n * (size_t)n, sizeof(double)); |
| Memcpy(xvals, REAL(x), (size_t) n * n); |
| } |
| PROTECT(x); |
| |
| ov = asLogical(only_values); |
| if (ov == NA_LOGICAL) error(_("invalid '%s' argument"), "only.values"); |
| vectors = !ov; |
| left = right = (double *) 0; |
| if (vectors) { |
| jobVR[0] = 'V'; |
| right = (double *) R_alloc(n * (size_t)n, sizeof(double)); |
| } |
| wR = (double *) R_alloc(n, sizeof(double)); |
| wI = (double *) R_alloc(n, sizeof(double)); |
| /* ask for optimal size of work array */ |
| lwork = -1; |
| F77_CALL(dgeev)(jobVL, jobVR, &n, xvals, &n, wR, wI, |
| left, &n, right, &n, &tmp, &lwork, &info FCONE FCONE); |
| if (info != 0) |
| error(_("error code %d from Lapack routine '%s'"), info, "dgeev"); |
| lwork = (int) tmp; |
| work = (double *) R_alloc(lwork, sizeof(double)); |
| F77_CALL(dgeev)(jobVL, jobVR, &n, xvals, &n, wR, wI, |
| left, &n, right, &n, work, &lwork, &info FCONE FCONE); |
| if (info != 0) |
| error(_("error code %d from Lapack routine '%s'"), info, "dgeev"); |
| |
| complexValues = FALSE; |
| for (i = 0; i < n; i++) |
| /* This test used to be !=0 for R < 2.3.0. This is OK for 0+0i */ |
| if (fabs(wI[i]) > 10 * R_AccuracyInfo.eps * fabs(wR[i])) { |
| complexValues = TRUE; |
| break; |
| } |
| SEXP ret = PROTECT(allocVector(VECSXP, 2)); |
| SEXP nm = PROTECT(allocVector(STRSXP, 2)); |
| SET_STRING_ELT(nm, 0, mkChar("values")); |
| SET_STRING_ELT(nm, 1, mkChar("vectors")); |
| setAttrib(ret, R_NamesSymbol, nm); |
| SET_VECTOR_ELT(ret, 1, R_NilValue); |
| if (complexValues) { |
| SEXP val = allocVector(CPLXSXP, n); |
| for (i = 0; i < n; i++) { |
| COMPLEX(val)[i].r = wR[i]; |
| COMPLEX(val)[i].i = wI[i]; |
| } |
| SET_VECTOR_ELT(ret, 0, val); |
| if (vectors) SET_VECTOR_ELT(ret, 1, unscramble(wI, n, right)); |
| } else { |
| SEXP val = allocVector(REALSXP, n); |
| for (i = 0; i < n; i++) REAL(val)[i] = wR[i]; |
| SET_VECTOR_ELT(ret, 0, val); |
| if(vectors) { |
| val = allocMatrix(REALSXP, n, n); |
| for (i = 0; i < (n * n); i++) REAL(val)[i] = right[i]; |
| SET_VECTOR_ELT(ret, 1, val); |
| } |
| } |
| UNPROTECT(3); |
| return ret; |
| } |
| |
| /* norm() except for 2-norm */ |
| static SEXP La_dlange(SEXP A, SEXP type) |
| { |
| int *xdims, m, n, nprot = 1; |
| double *work = NULL; |
| char typNorm[] = {'\0', '\0'}; |
| |
| if (!isMatrix(A)) error(_("'A' must be a numeric matrix")); |
| if (!isString(type)) |
| error(_("'type' must be a character string")); |
| if (!isReal(A)) { |
| A = PROTECT(coerceVector(A, REALSXP)); nprot++; |
| } |
| |
| xdims = INTEGER(coerceVector(getAttrib(A, R_DimSymbol), INTSXP)); |
| m = xdims[0]; |
| n = xdims[1]; /* m x n matrix {using Lapack naming convention} */ |
| |
| typNorm[0] = La_norm_type(CHAR(asChar(type))); |
| |
| SEXP val = PROTECT(allocVector(REALSXP, 1)); |
| if(*typNorm == 'I') work = (double *) R_alloc(m, sizeof(double)); |
| REAL(val)[0] = F77_CALL(dlange)(typNorm, &m, &n, REAL(A), &m, work FCONE); |
| |
| UNPROTECT(nprot); |
| return val; |
| } |
| |
| |
| /* ------------------------------------------------------------ */ |
| /* Real case of rcond, for general matrices */ |
| static SEXP La_dgecon(SEXP A, SEXP norm) |
| { |
| int *xdims, m, n, info, *iwork; |
| double anorm, *work; |
| char typNorm[] = {'\0', '\0'}; |
| |
| if (!isMatrix(A)) error(_("'A' must be a numeric matrix")); |
| if (!isString(norm)) |
| error(_("'norm' must be a character string")); |
| A = PROTECT(isReal(A) ? duplicate(A) : coerceVector(A, REALSXP)); |
| |
| xdims = INTEGER(coerceVector(getAttrib(A, R_DimSymbol), INTSXP)); |
| m = xdims[0]; n = xdims[1]; |
| |
| typNorm[0] = La_rcond_type(CHAR(asChar(norm))); |
| |
| SEXP val = PROTECT(allocVector(REALSXP, 1)); |
| work = (double *) R_alloc((*typNorm == 'I' && m > 4*(size_t)n) ? m : 4*(size_t)n, |
| sizeof(double)); |
| iwork = (int *) R_alloc(m, sizeof(int)); |
| |
| anorm = F77_CALL(dlange)(typNorm, &m, &n, REAL(A), &m, work FCONE); |
| |
| /* Compute the LU-decomposition and overwrite 'A' with result :*/ |
| F77_CALL(dgetrf)(&m, &n, REAL(A), &m, iwork, &info); |
| if (info) { |
| if (info < 0) { |
| UNPROTECT(2); |
| error(_("error [%d] from Lapack 'dgetrf()'"), info); |
| } |
| else { /* i := info > 0: LU decomp. is completed, but U[i,i] = 0 |
| * <==> singularity */ |
| #if 0 |
| warning(_("exact singularity: U[%d,%d] = 0 in LU-decomposition {Lapack 'dgetrf()'}"), |
| info,info); |
| #endif |
| REAL(val)[0] = 0.; /* rcond = 0 <==> singularity */ |
| UNPROTECT(2); |
| return val; |
| } |
| } |
| F77_CALL(dgecon)(typNorm, &n, REAL(A), &n, &anorm, |
| /* rcond = */ REAL(val), |
| work, iwork, &info FCONE); |
| UNPROTECT(2); |
| if (info) error(_("error [%d] from Lapack 'dgecon()'"), info); |
| return val; |
| } |
| |
| /* Real case of kappa.tri (and also rcond for a triangular matrix) */ |
| static SEXP La_dtrcon(SEXP A, SEXP norm) |
| { |
| int *xdims, n, nprot = 0, info; |
| char typNorm[] = {'\0', '\0'}; |
| |
| if (!isMatrix(A)) error(_("'A' must be a numeric matrix")); |
| if (!isString(norm)) error(_("'norm' must be a character string")); |
| if (!isReal(A)) { |
| nprot++; |
| A = PROTECT(coerceVector(A, REALSXP)); |
| } |
| xdims = INTEGER(coerceVector(getAttrib(A, R_DimSymbol), INTSXP)); |
| n = xdims[0]; |
| if(n != xdims[1]) { |
| UNPROTECT(nprot); |
| error(_("'A' must be a *square* matrix")); |
| } |
| |
| typNorm[0] = La_rcond_type(CHAR(asChar(norm))); |
| |
| nprot++; |
| SEXP val = PROTECT(allocVector(REALSXP, 1)); |
| |
| F77_CALL(dtrcon)(typNorm, "U", "N", &n, REAL(A), &n, |
| REAL(val), |
| /* work : */ (double *) R_alloc(3*(size_t)n, sizeof(double)), |
| /* iwork: */ (int *) R_alloc(n, sizeof(int)), |
| &info FCONE FCONE FCONE); |
| UNPROTECT(nprot); |
| if (info) error(_("error [%d] from Lapack 'dtrcon()'"), info); |
| return val; |
| } |
| |
| /* Complex case of rcond, for general matrices */ |
| static SEXP La_zgecon(SEXP A, SEXP norm) |
| { |
| #ifdef HAVE_FORTRAN_DOUBLE_COMPLEX |
| Rcomplex *avals; /* copy of A, to be modified */ |
| int *dims, n, info; |
| double anorm, *rwork; |
| char typNorm[] = {'\0', '\0'}; |
| |
| if (!isString(norm)) error(_("'norm' must be a character string")); |
| if (!(isMatrix(A) && isComplex(A))) |
| error(_("'A' must be a complex matrix")); |
| dims = INTEGER(coerceVector(getAttrib(A, R_DimSymbol), INTSXP)); |
| n = dims[0]; |
| if(n != dims[1]) error(_("'A' must be a *square* matrix")); |
| |
| typNorm[0] = La_rcond_type(CHAR(asChar(norm))); |
| |
| SEXP val = PROTECT(allocVector(REALSXP, 1)); |
| |
| rwork = (double *) R_alloc(2*(size_t)n, sizeof(Rcomplex)); |
| anorm = F77_CALL(zlange)(typNorm, &n, &n, COMPLEX(A), &n, rwork FCONE); |
| |
| /* Compute the LU-decomposition and overwrite 'x' with result; |
| * working on a copy of A : */ |
| avals = (Rcomplex *) R_alloc((size_t)n * n, sizeof(Rcomplex)); |
| Memcpy(avals, COMPLEX(A), (size_t) n * n); |
| F77_CALL(zgetrf)(&n, &n, avals, &n, |
| /* iwork: */(int *) R_alloc(n, sizeof(int)), |
| &info); |
| if (info) { |
| if (info < 0) { |
| UNPROTECT(1); |
| error(_("error [%d] from Lapack 'zgetrf()'"), info); |
| } else { |
| REAL(val)[0] = 0.; /* rcond = 0 <==> singularity */ |
| UNPROTECT(1); |
| return val; |
| } |
| UNPROTECT(1); |
| error(_("error [%d] from Lapack 'zgetrf()'"), info); |
| } |
| |
| F77_CALL(zgecon)(typNorm, &n, avals, &n, &anorm, |
| /* rcond = */ REAL(val), |
| /* work : */ (Rcomplex *) R_alloc(4*(size_t)n, sizeof(Rcomplex)), |
| rwork, &info FCONE); |
| UNPROTECT(1); |
| if (info) error(_("error [%d] from Lapack 'zgecon()'"), info); |
| return val; |
| |
| #else |
| error(_("Fortran complex functions are not available on this platform")); |
| return R_NilValue; /* -Wall */ |
| #endif |
| } |
| |
| /* Complex case of kappa.tri (and also rcond for a triangular matrix) */ |
| static SEXP La_ztrcon(SEXP A, SEXP norm) |
| { |
| #ifdef HAVE_FORTRAN_DOUBLE_COMPLEX |
| SEXP val; |
| int *dims, n, info; |
| char typNorm[] = {'\0', '\0'}; |
| |
| if (!isString(norm)) |
| error(_("'norm' must be a character string")); |
| if (!(isMatrix(A) && isComplex(A))) |
| error(_("'A' must be a complex matrix")); |
| dims = INTEGER(coerceVector(getAttrib(A, R_DimSymbol), INTSXP)); |
| n = dims[0]; |
| if(n != dims[1]) |
| error(_("'A' must be a *square* matrix")); |
| |
| typNorm[0] = La_rcond_type(CHAR(asChar(norm))); |
| |
| val = PROTECT(allocVector(REALSXP, 1)); |
| |
| F77_CALL(ztrcon)(typNorm, "U", "N", &n, COMPLEX(A), &n, |
| REAL(val), |
| /* work : */ (Rcomplex *) R_alloc(2*(size_t)n, sizeof(Rcomplex)), |
| /* rwork: */ (double *) R_alloc(n, sizeof(double)), |
| &info FCONE FCONE FCONE); |
| UNPROTECT(1); |
| if (info) error(_("error [%d] from Lapack 'ztrcon()'"), info); |
| return val; |
| #else |
| error(_("Fortran complex functions are not available on this platform")); |
| return R_NilValue; /* -Wall */ |
| #endif |
| } |
| |
| /* Complex case of solve.default: see the comments in La_solve */ |
| static SEXP La_solve_cmplx(SEXP A, SEXP Bin) |
| { |
| #ifdef HAVE_FORTRAN_DOUBLE_COMPLEX |
| int n, p, info, *ipiv, *Adims, *Bdims; |
| Rcomplex *avals; |
| SEXP B, Adn, Bdn; |
| |
| if (!isMatrix(A)) error(_("'a' must be a complex matrix")); |
| Adims = INTEGER(coerceVector(getAttrib(A, R_DimSymbol), INTSXP)); |
| n = Adims[0]; |
| if(n == 0) error(_("'a' is 0-diml")); |
| int n2 = Adims[1]; |
| if(n2 != n) error(_("'a' (%d x %d) must be square"), n, n2); |
| Adn = getAttrib(A, R_DimNamesSymbol); |
| |
| if (isMatrix(Bin)) { |
| Bdims = INTEGER(coerceVector(getAttrib(Bin, R_DimSymbol), INTSXP)); |
| p = Bdims[1]; |
| if(p == 0) error(_("no right-hand side in 'b'")); |
| int p2 = Bdims[0]; |
| if(p2 != n) |
| error(_("'b' (%d x %d) must be compatible with 'a' (%d x %d)"), |
| p2, p, n, n); |
| PROTECT(B = allocMatrix(CPLXSXP, n, p)); |
| SEXP Bindn = getAttrib(Bin, R_DimNamesSymbol); |
| if (!isNull(Adn) || !isNull(Bindn)) { |
| Bdn = allocVector(VECSXP, 2); |
| if (!isNull(Adn)) SET_VECTOR_ELT(Bdn, 0, VECTOR_ELT(Adn, 1)); |
| if (!isNull(Bindn)) SET_VECTOR_ELT(Bdn, 1, VECTOR_ELT(Bindn, 1)); |
| if (!isNull(VECTOR_ELT(Bdn, 0)) && !isNull(VECTOR_ELT(Bdn, 1))) |
| setAttrib(B, R_DimNamesSymbol, Bdn); |
| } |
| } else { |
| p = 1; |
| if(length(Bin) != n) |
| error(_("'b' (%d x %d) must be compatible with 'a' (%d x %d)"), |
| length(Bin), p, n, n); |
| PROTECT(B = allocVector(CPLXSXP, n)); |
| if (!isNull(Adn)) setAttrib(B, R_NamesSymbol, VECTOR_ELT(Adn, 1)); |
| } |
| Bin = PROTECT(coerceVector(Bin, CPLXSXP)); |
| Memcpy(COMPLEX(B), COMPLEX(Bin), (size_t)n * p); |
| |
| ipiv = (int *) R_alloc(n, sizeof(int)); |
| |
| /* work on a copy of A */ |
| if(TYPEOF(A) != CPLXSXP) { |
| A = coerceVector(A, CPLXSXP); |
| avals = COMPLEX(A); |
| } else { |
| avals = (Rcomplex *) R_alloc((size_t)n * n, sizeof(Rcomplex)); |
| Memcpy(avals, COMPLEX(A), (size_t) n * n); |
| } |
| PROTECT(A); |
| F77_CALL(zgesv)(&n, &p, avals, &n, ipiv, COMPLEX(B), &n, &info); |
| if (info < 0) |
| error(_("argument %d of Lapack routine %s had invalid value"), |
| -info, "zgesv"); |
| if (info > 0) |
| error(("Lapack routine zgesv: system is exactly singular")); |
| UNPROTECT(3); /* B, Bin, A */ |
| return B; |
| #else |
| error(_("Fortran complex functions are not available on this platform")); |
| return R_NilValue; /* -Wall */ |
| #endif |
| } |
| |
| /* Complex case of qr.default */ |
| static SEXP La_qr_cmplx(SEXP Ain) |
| { |
| #ifdef HAVE_FORTRAN_DOUBLE_COMPLEX |
| int i, m, n, *Adims, info, lwork; |
| Rcomplex *work, tmp; |
| double *rwork; |
| |
| if (!(isMatrix(Ain) && isComplex(Ain))) |
| error(_("'a' must be a complex matrix")); |
| SEXP Adn = getAttrib(Ain, R_DimNamesSymbol); |
| Adims = INTEGER(coerceVector(getAttrib(Ain, R_DimSymbol), INTSXP)); |
| m = Adims[0]; n = Adims[1]; |
| SEXP A = PROTECT(allocMatrix(CPLXSXP, m, n)); |
| Memcpy(COMPLEX(A), COMPLEX(Ain), (size_t)m * n); |
| rwork = (double *) R_alloc(2*(size_t)n, sizeof(double)); |
| |
| SEXP jpvt = PROTECT(allocVector(INTSXP, n)); |
| for (i = 0; i < n; i++) INTEGER(jpvt)[i] = 0; |
| SEXP tau = PROTECT(allocVector(CPLXSXP, m < n ? m : n)); |
| lwork = -1; |
| F77_CALL(zgeqp3)(&m, &n, COMPLEX(A), &m, INTEGER(jpvt), COMPLEX(tau), |
| &tmp, &lwork, rwork, &info); |
| if (info != 0) |
| error(_("error code %d from Lapack routine '%s'"), info, "zgeqp3"); |
| lwork = (int) tmp.r; |
| work = (Rcomplex *) R_alloc(lwork, sizeof(Rcomplex)); |
| F77_CALL(zgeqp3)(&m, &n, COMPLEX(A), &m, INTEGER(jpvt), COMPLEX(tau), |
| work, &lwork, rwork, &info); |
| if (info != 0) |
| error(_("error code %d from Lapack routine '%s'"), info, "zgeqp3"); |
| SEXP val = PROTECT(allocVector(VECSXP, 4)); |
| SEXP nm = PROTECT(allocVector(STRSXP, 4)); |
| SET_STRING_ELT(nm, 0, mkChar("qr")); |
| SET_STRING_ELT(nm, 1, mkChar("rank")); |
| SET_STRING_ELT(nm, 2, mkChar("qraux")); |
| SET_STRING_ELT(nm, 3, mkChar("pivot")); |
| setAttrib(val, R_NamesSymbol, nm); |
| // Fix up dimnames(A) |
| if(!isNull(Adn)) { |
| SEXP Adn2 = duplicate(Adn); |
| SEXP cn = VECTOR_ELT(Adn, 1), cn2 = VECTOR_ELT(Adn2, 1); |
| if(!isNull(cn)) { // pivot them |
| for (int j = 0; j < n; j++) |
| SET_STRING_ELT(cn2, j, STRING_ELT(cn, INTEGER(jpvt)[j]-1)); |
| } |
| setAttrib(A, R_DimNamesSymbol, Adn2); |
| } |
| SET_VECTOR_ELT(val, 0, A); |
| SET_VECTOR_ELT(val, 1, ScalarInteger(m < n ? m : n)); |
| SET_VECTOR_ELT(val, 2, tau); |
| SET_VECTOR_ELT(val, 3, jpvt); |
| UNPROTECT(5); |
| return val; |
| #else |
| error(_("Fortran complex functions are not available on this platform")); |
| return R_NilValue; /* -Wall */ |
| #endif |
| } |
| |
| /* Complex case of qr_coef */ |
| static SEXP qr_coef_cmplx(SEXP Q, SEXP Bin) |
| { |
| #ifdef HAVE_FORTRAN_DOUBLE_COMPLEX |
| int n, nrhs, lwork, info, k, *Bdims; |
| SEXP B, qr = VECTOR_ELT(Q, 0), tau = VECTOR_ELT(Q, 2); |
| Rcomplex *work, tmp; |
| |
| k = LENGTH(tau); |
| if (!isMatrix(Bin)) error(_("'b' must be a complex matrix")); |
| |
| if (!isComplex(Bin)) B = PROTECT(coerceVector(Bin, CPLXSXP)); |
| else B = PROTECT(duplicate(Bin)); |
| |
| n = INTEGER(coerceVector(getAttrib(qr, R_DimSymbol), INTSXP))[0]; |
| Bdims = INTEGER(coerceVector(getAttrib(Bin, R_DimSymbol), INTSXP)); |
| if(Bdims[0] != n) |
| error(_("right-hand side should have %d not %d rows"), n, Bdims[0]); |
| nrhs = Bdims[1]; |
| lwork = -1; |
| F77_CALL(zunmqr)("L", "C", &n, &nrhs, &k, |
| COMPLEX(qr), &n, COMPLEX(tau), COMPLEX(B), &n, |
| &tmp, &lwork, &info FCONE FCONE); |
| if (info != 0) |
| error(_("error code %d from Lapack routine '%s'"), info, "zunmqr"); |
| lwork = (int) tmp.r; |
| work = (Rcomplex *) R_alloc(lwork, sizeof(Rcomplex)); |
| F77_CALL(zunmqr)("L", "C", &n, &nrhs, &k, |
| COMPLEX(qr), &n, COMPLEX(tau), COMPLEX(B), &n, |
| work, &lwork, &info FCONE FCONE); |
| if (info != 0) |
| error(_("error code %d from Lapack routine '%s'"), info, "zunmqr"); |
| F77_CALL(ztrtrs)("U", "N", "N", &k, &nrhs, |
| COMPLEX(qr), &n, COMPLEX(B), &n, &info |
| FCONE FCONE FCONE); |
| if (info != 0) |
| error(_("error code %d from Lapack routine '%s'"), info, "ztrtrs"); |
| UNPROTECT(1); |
| return B; |
| #else |
| error(_("Fortran complex functions are not available on this platform")); |
| return R_NilValue; /* -Wall */ |
| #endif |
| } |
| |
| /* Complex case of qr.qy and qr.qty */ |
| static SEXP qr_qy_cmplx(SEXP Q, SEXP Bin, SEXP trans) |
| { |
| #ifdef HAVE_FORTRAN_DOUBLE_COMPLEX |
| int n, nrhs, lwork, info, k, *Bdims, tr; |
| SEXP B, qr = VECTOR_ELT(Q, 0), tau = VECTOR_ELT(Q, 2); |
| Rcomplex *work, tmp; |
| |
| k = LENGTH(tau); |
| if (!(isMatrix(Bin) && isComplex(Bin))) |
| error(_("'b' must be a complex matrix")); |
| tr = asLogical(trans); |
| if(tr == NA_LOGICAL) error(_("invalid '%s' argument"), "trans"); |
| |
| if (!isReal(Bin)) B = PROTECT(coerceVector(Bin, CPLXSXP)); |
| else B = PROTECT(duplicate(Bin)); |
| n = INTEGER(coerceVector(getAttrib(qr, R_DimSymbol), INTSXP))[0]; |
| Bdims = INTEGER(coerceVector(getAttrib(B, R_DimSymbol), INTSXP)); |
| if(Bdims[0] != n) |
| error(_("right-hand side should have %d not %d rows"), n, Bdims[0]); |
| nrhs = Bdims[1]; |
| lwork = -1; |
| F77_CALL(zunmqr)("L", tr ? "C" : "N", &n, &nrhs, &k, |
| COMPLEX(qr), &n, COMPLEX(tau), COMPLEX(B), &n, |
| &tmp, &lwork, &info FCONE FCONE); |
| if (info != 0) |
| error(_("error code %d from Lapack routine '%s'"), info, "zunmqr"); |
| lwork = (int) tmp.r; |
| work = (Rcomplex *) R_alloc(lwork, sizeof(Rcomplex)); |
| F77_CALL(zunmqr)("L", tr ? "C" : "N", &n, &nrhs, &k, |
| COMPLEX(qr), &n, COMPLEX(tau), COMPLEX(B), &n, |
| work, &lwork, &info FCONE FCONE); |
| if (info != 0) |
| error(_("error code %d from Lapack routine '%s'"), info, "zunmqr"); |
| UNPROTECT(1); |
| return B; |
| #else |
| error(_("Fortran complex functions are not available on this platform")); |
| return R_NilValue; /* -Wall */ |
| #endif |
| } |
| |
| static SEXP La_svd_cmplx(SEXP jobu, SEXP x, SEXP s, SEXP u, SEXP v) |
| { |
| #ifdef HAVE_FORTRAN_DOUBLE_COMPLEX |
| if (!(isString(jobu))) |
| error(_("'jobu' must be a character string")); |
| int *xdims = INTEGER(coerceVector(getAttrib(x, R_DimSymbol), INTSXP)); |
| int n = xdims[0], p = xdims[1]; |
| const char *jz = CHAR(STRING_ELT(jobu, 0)); |
| |
| /* The underlying LAPACK, specifically ZLARF, does not work with |
| * long arrays */ |
| if ((double)n * (double)p > INT_MAX) |
| error(_("matrices of 2^31 or more elements are not supported")); |
| |
| /* work on a copy of x */ |
| Rcomplex *xvals = (Rcomplex *) R_alloc(n * (size_t) p, sizeof(Rcomplex)); |
| Memcpy(xvals, COMPLEX(x), n * (size_t) p); |
| |
| int *iwork= (int *) R_alloc(8*(size_t)(n < p ? n : p), sizeof(int)); |
| size_t mn0 = (n < p ? n : p), mn1 = (n > p ? n : p), lrwork; |
| if (strcmp(jz, "N")) { |
| size_t f1 = 5 * mn1 + 7, f2 = 2 * mn1 + 2 * mn0 + 1; |
| lrwork = (f1 > f2 ? f1 : f2) * mn0; |
| // 7 replaces 5: bug 111 in http://www.netlib.org/lapack/Errata/index2.html |
| } else lrwork = 7 * mn0; |
| double *rwork = (double *) R_alloc(lrwork, sizeof(double)); |
| /* ask for optimal size of lwork array */ |
| int lwork = -1, info; |
| Rcomplex tmp; |
| int ldu, ldv; |
| SEXP dims = getAttrib(u, R_DimSymbol); |
| if (TYPEOF(dims) != INTSXP) error("non-integer dims"); |
| ldu = INTEGER(dims)[0]; |
| dims = getAttrib(v, R_DimSymbol); |
| if (TYPEOF(dims) != INTSXP) error("non-integer dims"); |
| ldv = INTEGER(dims)[0]; |
| F77_CALL(zgesdd)(jz, &n, &p, xvals, &n, REAL(s), |
| COMPLEX(u), &ldu, COMPLEX(v), &ldv, |
| &tmp, &lwork, rwork, iwork, &info FCONE); |
| if (info != 0) |
| error(_("error code %d from Lapack routine '%s'"), info, "zgesdd"); |
| lwork = (int) tmp.r; |
| Rcomplex *work = (Rcomplex *) R_alloc(lwork, sizeof(Rcomplex)); |
| F77_CALL(zgesdd)(jz, &n, &p, xvals, &n, REAL(s), |
| COMPLEX(u), &ldu, COMPLEX(v), &ldv, |
| work, &lwork, rwork, iwork, &info FCONE); |
| if (info != 0) |
| error(_("error code %d from Lapack routine '%s'"), info, "zgesdd"); |
| |
| SEXP val = PROTECT(allocVector(VECSXP, 3)); |
| SEXP nm = PROTECT(allocVector(STRSXP, 3)); |
| SET_STRING_ELT(nm, 0, mkChar("d")); |
| SET_STRING_ELT(nm, 1, mkChar("u")); |
| SET_STRING_ELT(nm, 2, mkChar("vt")); |
| setAttrib(val, R_NamesSymbol, nm); |
| SET_VECTOR_ELT(val, 0, s); |
| SET_VECTOR_ELT(val, 1, u); |
| SET_VECTOR_ELT(val, 2, v); |
| UNPROTECT(2); |
| return val; |
| #else |
| error(_("Fortran complex functions are not available on this platform")); |
| return R_NilValue; /* -Wall */ |
| #endif |
| } |
| |
| /* Complex, symmetric case of eigen */ |
| static SEXP La_rs_cmplx(SEXP xin, SEXP only_values) |
| { |
| #ifdef HAVE_FORTRAN_DOUBLE_COMPLEX |
| int *xdims, n, lwork, info, ov; |
| char jobv[2] = "N", uplo[2] = "L"; |
| Rcomplex *work, *rx, tmp; |
| double *rwork, *rvalues; |
| |
| xdims = INTEGER(coerceVector(getAttrib(xin, R_DimSymbol), INTSXP)); |
| n = xdims[0]; |
| if (n != xdims[1]) error(_("'x' must be a square complex matrix")); |
| ov = asLogical(only_values); |
| if (ov == NA_LOGICAL) error(_("invalid '%s' argument"), "only.values"); |
| if (ov) jobv[0] = 'N'; else jobv[0] = 'V'; |
| |
| SEXP x = PROTECT(allocMatrix(CPLXSXP, n, n)); |
| rx = COMPLEX(x); |
| Memcpy(rx, COMPLEX(xin), (size_t) n * n); |
| SEXP values = PROTECT(allocVector(REALSXP, n)); |
| rvalues = REAL(values); |
| |
| rwork = (double *) R_alloc((3*(size_t)n-2) > 1 ? 3*(size_t)n-2 : 1, |
| sizeof(double)); |
| /* ask for optimal size of work array */ |
| lwork = -1; |
| F77_CALL(zheev)(jobv, uplo, &n, rx, &n, rvalues, &tmp, &lwork, rwork, |
| &info FCONE FCONE); |
| if (info != 0) |
| error(_("error code %d from Lapack routine '%s'"), info, "zheev"); |
| lwork = (int) tmp.r; |
| work = (Rcomplex *) R_alloc(lwork, sizeof(Rcomplex)); |
| F77_CALL(zheev)(jobv, uplo, &n, rx, &n, rvalues, work, &lwork, rwork, |
| &info FCONE FCONE); |
| if (info != 0) |
| error(_("error code %d from Lapack routine '%s'"), info, "zheev"); |
| SEXP ret, nm; |
| if (!ov) { |
| ret = PROTECT(allocVector(VECSXP, 2)); |
| nm = PROTECT(allocVector(STRSXP, 2)); |
| SET_STRING_ELT(nm, 1, mkChar("vectors")); |
| SET_VECTOR_ELT(ret, 1, x); |
| } else { |
| ret = PROTECT(allocVector(VECSXP, 1)); |
| nm = PROTECT(allocVector(STRSXP, 1)); |
| } |
| SET_STRING_ELT(nm, 0, mkChar("values")); |
| setAttrib(ret, R_NamesSymbol, nm); |
| SET_VECTOR_ELT(ret, 0, values); |
| UNPROTECT(4); |
| return ret; |
| #else |
| error(_("Fortran complex functions are not available on this platform")); |
| return R_NilValue; /* -Wall */ |
| #endif |
| } |
| |
| /* Complex, general case of eigen */ |
| static SEXP La_rg_cmplx(SEXP x, SEXP only_values) |
| { |
| #ifdef HAVE_FORTRAN_DOUBLE_COMPLEX |
| int n, lwork, info, *xdims, ov; |
| Rcomplex *work, *left, *right, *xvals, tmp; |
| double *rwork; |
| char jobVL[2] = "N", jobVR[2] = "N"; |
| SEXP ret, nm, values, val = R_NilValue; |
| |
| xdims = INTEGER(coerceVector(getAttrib(x, R_DimSymbol), INTSXP)); |
| n = xdims[0]; |
| if (n != xdims[1]) error(_("'x' must be a square numeric matrix")); |
| |
| /* work on a copy of x */ |
| xvals = (Rcomplex *) R_alloc((size_t)n * n, sizeof(Rcomplex)); |
| Memcpy(xvals, COMPLEX(x), (size_t) n * n); |
| ov = asLogical(only_values); |
| if (ov == NA_LOGICAL) error(_("invalid '%s' argument"), "only.values"); |
| left = right = (Rcomplex *) 0; |
| if (!ov) { |
| jobVR[0] = 'V'; |
| PROTECT(val = allocMatrix(CPLXSXP, n, n)); |
| right = COMPLEX(val); |
| } |
| PROTECT(values = allocVector(CPLXSXP, n)); |
| rwork = (double *) R_alloc(2*(size_t)n, sizeof(double)); |
| /* ask for optimal size of work array */ |
| lwork = -1; |
| F77_CALL(zgeev)(jobVL, jobVR, &n, xvals, &n, COMPLEX(values), |
| left, &n, right, &n, &tmp, &lwork, rwork, &info |
| FCONE FCONE); |
| if (info != 0) |
| error(_("error code %d from Lapack routine '%s'"), info, "zgeev"); |
| lwork = (int) tmp.r; |
| work = (Rcomplex *) R_alloc(lwork, sizeof(Rcomplex)); |
| F77_CALL(zgeev)(jobVL, jobVR, &n, xvals, &n, COMPLEX(values), |
| left, &n, right, &n, work, &lwork, rwork, &info |
| FCONE FCONE); |
| if (info != 0) |
| error(_("error code %d from Lapack routine '%s'"), info, "zgeev"); |
| |
| if(!ov){ |
| ret = PROTECT(allocVector(VECSXP, 2)); |
| nm = PROTECT(allocVector(STRSXP, 2)); |
| SET_STRING_ELT(nm, 1, mkChar("vectors")); |
| SET_VECTOR_ELT(ret, 1, val); |
| } else { |
| ret = PROTECT(allocVector(VECSXP, 1)); |
| nm = PROTECT(allocVector(STRSXP, 1)); |
| } |
| SET_STRING_ELT(nm, 0, mkChar("values")); |
| SET_VECTOR_ELT(ret, 0, values); |
| setAttrib(ret, R_NamesSymbol, nm); |
| UNPROTECT(ov ? 3 : 4); |
| return ret; |
| #else |
| error(_("Fortran complex functions are not available on this platform")); |
| return R_NilValue; /* -Wall */ |
| #endif |
| } |
| |
| /* ------------------------------------------------------------ */ |
| |
| static SEXP La_chol(SEXP A, SEXP pivot, SEXP stol) |
| { |
| if (!isMatrix(A)) error(_("'a' must be a numeric matrix")); |
| |
| SEXP ans = PROTECT(isReal(A) ? duplicate(A): coerceVector(A, REALSXP)); |
| SEXP adims = getAttrib(A, R_DimSymbol); |
| if (TYPEOF(adims) != INTSXP) error("non-integer dims"); |
| int m = INTEGER(adims)[0], n = INTEGER(adims)[1]; |
| |
| if (m != n) error(_("'a' must be a square matrix")); |
| if (m <= 0) error(_("'a' must have dims > 0")); |
| size_t N = n; |
| for (int j = 0; j < n; j++) /* zero the lower triangle */ |
| for (int i = j+1; i < n; i++) REAL(ans)[i + N * j] = 0.; |
| |
| int piv = asInteger(pivot); |
| if (piv != 0 && piv != 1) error("invalid '%s' value", "pivot"); |
| if(!piv) { |
| int info; |
| F77_CALL(dpotrf)("U", &m, REAL(ans), &m, &info FCONE); |
| if (info != 0) { |
| if (info > 0) |
| error(_("the leading minor of order %d is not positive definite"), |
| info); |
| error(_("argument %d of Lapack routine %s had invalid value"), |
| -info, "dpotrf"); |
| } |
| } else { |
| double tol = asReal(stol); |
| SEXP piv = PROTECT(allocVector(INTSXP, m)); |
| int *ip = INTEGER(piv); |
| double *work = (double *) R_alloc(2 * (size_t)m, sizeof(double)); |
| int rank, info; |
| F77_CALL(dpstrf)("U", &m, REAL(ans), &m, ip, &rank, &tol, work, &info |
| FCONE); |
| if (info != 0) { |
| if (info > 0) |
| warning(_("the matrix is either rank-deficient or indefinite")); |
| else |
| error(_("argument %d of Lapack routine %s had invalid value"), |
| -info, "dpstrf"); |
| } |
| setAttrib(ans, install("pivot"), piv); |
| SEXP s_rank = install("rank"); |
| setAttrib(ans, s_rank, ScalarInteger(rank)); |
| SEXP cn, dn = getAttrib(ans, R_DimNamesSymbol); |
| if (!isNull(dn) && !isNull(cn = VECTOR_ELT(dn, 1))) { |
| // need to pivot the colnames |
| SEXP dn2 = PROTECT(duplicate(dn)); |
| SEXP cn2 = VECTOR_ELT(dn2, 1); |
| for(int i = 0; i < m; i++) |
| SET_STRING_ELT(cn2, i, STRING_ELT(cn, ip[i] - 1)); // base 1 |
| setAttrib(ans, R_DimNamesSymbol, dn2); |
| UNPROTECT(1); |
| } |
| UNPROTECT(1); // piv |
| } |
| UNPROTECT(1); // ans |
| return ans; |
| } |
| |
| static SEXP La_chol2inv(SEXP A, SEXP size) |
| { |
| int sz = asInteger(size); |
| if (sz == NA_INTEGER || sz < 1) { |
| error(_("'size' argument must be a positive integer")); |
| return R_NilValue; /* -Wall */ |
| } else { |
| SEXP ans, Amat = A; /* -Wall: we initialize here as for the 1x1 case */ |
| int m = 1, n = 1, nprot = 0; |
| |
| if (sz == 1 && !isMatrix(A) && isReal(A)) { |
| /* nothing to do; m = n = 1; ... */ |
| } else if (isMatrix(A)) { |
| SEXP adims = getAttrib(A, R_DimSymbol); |
| if (TYPEOF(adims) != INTSXP) error("non-integer dims"); |
| Amat = PROTECT(coerceVector(A, REALSXP)); nprot++; |
| m = INTEGER(adims)[0]; n = INTEGER(adims)[1]; |
| } else error(_("'a' must be a numeric matrix")); |
| |
| if (sz > n) { UNPROTECT(nprot); error(_("'size' cannot exceed ncol(x) = %d"), n); } |
| if (sz > m) { UNPROTECT(nprot); error(_("'size' cannot exceed nrow(x) = %d"), m); } |
| ans = PROTECT(allocMatrix(REALSXP, sz, sz)); nprot++; |
| size_t M = m, SZ = sz; |
| for (int j = 0; j < sz; j++) { |
| for (int i = 0; i <= j; i++) |
| REAL(ans)[i + j * SZ] = REAL(Amat)[i + j * M]; |
| } |
| int info; |
| F77_CALL(dpotri)("U", &sz, REAL(ans), &sz, &info FCONE); |
| if (info != 0) { |
| UNPROTECT(nprot); |
| if (info > 0) |
| error(_("element (%d, %d) is zero, so the inverse cannot be computed"), |
| info, info); |
| error(_("argument %d of Lapack routine %s had invalid value"), |
| -info, "dpotri"); |
| } |
| for (int j = 0; j < sz; j++) |
| for (int i = j+1; i < sz; i++) |
| REAL(ans)[i + j * SZ] = REAL(ans)[j + i * SZ]; |
| UNPROTECT(nprot); |
| return ans; |
| } |
| } |
| |
| /* ------------------------------------------------------------ */ |
| |
| /* Real case of solve.default */ |
| static SEXP La_solve(SEXP A, SEXP Bin, SEXP tolin) |
| { |
| int n, p; |
| double *avals, anorm, rcond, tol = asReal(tolin), *work; |
| SEXP B, Adn, Bdn; |
| |
| if (!(isMatrix(A) && |
| (TYPEOF(A) == REALSXP || TYPEOF(A) == INTSXP || TYPEOF(A) == LGLSXP))) |
| error(_("'a' must be a numeric matrix")); |
| int *Adims = INTEGER(coerceVector(getAttrib(A, R_DimSymbol), INTSXP)); |
| n = Adims[0]; |
| if(n == 0) error(_("'a' is 0-diml")); |
| int n2 = Adims[1]; |
| if(n2 != n) error(_("'a' (%d x %d) must be square"), n, n2); |
| Adn = getAttrib(A, R_DimNamesSymbol); |
| |
| if (isMatrix(Bin)) { |
| int *Bdims = INTEGER(coerceVector(getAttrib(Bin, R_DimSymbol), INTSXP)); |
| p = Bdims[1]; |
| if(p == 0) error(_("no right-hand side in 'b'")); |
| int p2 = Bdims[0]; |
| if(p2 != n) |
| error(_("'b' (%d x %d) must be compatible with 'a' (%d x %d)"), |
| p2, p, n, n); |
| PROTECT(B = allocMatrix(REALSXP, n, p)); |
| SEXP Bindn = getAttrib(Bin, R_DimNamesSymbol); |
| // This is somewhat odd, but Matrix relies on dropping NULL dimnames |
| if (!isNull(Adn) || !isNull(Bindn)) { |
| // rownames(ans) = colnames(A), colnames(ans) = colnames(Bin) |
| Bdn = allocVector(VECSXP, 2); |
| if (!isNull(Adn)) SET_VECTOR_ELT(Bdn, 0, VECTOR_ELT(Adn, 1)); |
| if (!isNull(Bindn)) SET_VECTOR_ELT(Bdn, 1, VECTOR_ELT(Bindn, 1)); |
| if (!isNull(VECTOR_ELT(Bdn, 0)) || !isNull(VECTOR_ELT(Bdn, 1))) |
| setAttrib(B, R_DimNamesSymbol, Bdn); |
| } |
| } else { |
| p = 1; |
| if(length(Bin) != n) |
| error(_("'b' (%d x %d) must be compatible with 'a' (%d x %d)"), |
| length(Bin), p, n, n); |
| PROTECT(B = allocVector(REALSXP, n)); |
| if (!isNull(Adn)) setAttrib(B, R_NamesSymbol, VECTOR_ELT(Adn, 1)); |
| } |
| PROTECT(Bin = coerceVector(Bin, REALSXP)); |
| Memcpy(REAL(B), REAL(Bin), (size_t)n * p); |
| |
| int *ipiv = (int *) R_alloc(n, sizeof(int)); |
| |
| /* work on a copy of A */ |
| if (!isReal(A)) { |
| A = coerceVector(A, REALSXP); |
| avals = REAL(A); |
| } else { |
| avals = (double *) R_alloc((size_t)n * n, sizeof(double)); |
| Memcpy(avals, REAL(A), (size_t)n * n); |
| } |
| PROTECT(A); |
| int info; |
| F77_CALL(dgesv)(&n, &p, avals, &n, ipiv, REAL(B), &n, &info); |
| if (info < 0) |
| error(_("argument %d of Lapack routine %s had invalid value"), |
| -info, "dgesv"); |
| if (info > 0) |
| error(_("Lapack routine %s: system is exactly singular: U[%d,%d] = 0"), |
| "dgesv", info, info); |
| if(tol > 0) { |
| char one[2] = "1"; |
| anorm = F77_CALL(dlange)(one, &n, &n, REAL(A), &n, |
| (double*) NULL FCONE); |
| work = (double *) R_alloc(4*(size_t)n, sizeof(double)); |
| F77_CALL(dgecon)(one, &n, avals, &n, &anorm, &rcond, work, ipiv, |
| &info FCONE); |
| if (rcond < tol) |
| error(_("system is computationally singular: reciprocal condition number = %g"), |
| rcond); |
| } |
| UNPROTECT(3); /* B, Bin, A */ |
| return B; |
| } |
| |
| /* Real case of qr.default */ |
| static SEXP La_qr(SEXP Ain) |
| { |
| int m, n; |
| |
| if (!isMatrix(Ain)) error(_("'a' must be a numeric matrix")); |
| SEXP Adn = getAttrib(Ain, R_DimNamesSymbol); |
| int *Adims = INTEGER(coerceVector(getAttrib(Ain, R_DimSymbol), INTSXP)); |
| m = Adims[0]; n = Adims[1]; |
| SEXP A; |
| if (!isReal(Ain)) { |
| A = PROTECT(coerceVector(Ain, REALSXP)); |
| } else { |
| A = PROTECT(allocMatrix(REALSXP, m, n)); |
| Memcpy(REAL(A), REAL(Ain), (size_t)m * n); |
| } |
| |
| SEXP jpvt = PROTECT(allocVector(INTSXP, n)); |
| for (int i = 0; i < n; i++) INTEGER(jpvt)[i] = 0; |
| SEXP tau = PROTECT(allocVector(REALSXP, m < n ? m : n)); |
| int info, lwork = -1; |
| double tmp; |
| F77_CALL(dgeqp3)(&m, &n, REAL(A), &m, INTEGER(jpvt), REAL(tau), |
| &tmp, &lwork, &info); |
| if (info < 0) |
| error(_("error code %d from Lapack routine '%s'"), info, "dgeqp3"); |
| lwork = (int) tmp; |
| double *work = (double *) R_alloc(lwork, sizeof(double)); |
| F77_CALL(dgeqp3)(&m, &n, REAL(A), &m, INTEGER(jpvt), REAL(tau), |
| work, &lwork, &info); |
| if (info < 0) |
| error(_("error code %d from Lapack routine '%s'"), info, "dgeqp3"); |
| SEXP val = PROTECT(allocVector(VECSXP, 4)); |
| SEXP nm = PROTECT(allocVector(STRSXP, 4)); |
| SET_STRING_ELT(nm, 0, mkChar("qr")); |
| SET_STRING_ELT(nm, 1, mkChar("rank")); |
| SET_STRING_ELT(nm, 2, mkChar("qraux")); |
| SET_STRING_ELT(nm, 3, mkChar("pivot")); |
| setAttrib(val, R_NamesSymbol, nm); |
| // Fix up dimnames(A) |
| if(!isNull(Adn)) { |
| SEXP Adn2 = duplicate(Adn); |
| SEXP cn = VECTOR_ELT(Adn, 1), cn2 = VECTOR_ELT(Adn2, 1); |
| if(!isNull(cn)) { // pivot them |
| for (int j = 0; j < n; j++) |
| SET_STRING_ELT(cn2, j, STRING_ELT(cn, INTEGER(jpvt)[j]-1)); |
| } |
| setAttrib(A, R_DimNamesSymbol, Adn2); |
| } |
| SET_VECTOR_ELT(val, 0, A); |
| SET_VECTOR_ELT(val, 1, ScalarInteger(m < n ? m : n)); |
| SET_VECTOR_ELT(val, 2, tau); |
| SET_VECTOR_ELT(val, 3, jpvt); |
| UNPROTECT(5); |
| return val; |
| } |
| |
| /* Real case of qr.coef */ |
| static SEXP qr_coef_real(SEXP Q, SEXP Bin) |
| { |
| if (!isMatrix(Bin)) error(_("'b' must be a numeric matrix")); |
| |
| SEXP B = PROTECT(isReal(Bin) ? duplicate(Bin) : coerceVector(Bin, REALSXP)), |
| qr = VECTOR_ELT(Q, 0), // qr$qr |
| tau = VECTOR_ELT(Q, 2); // qr$qraux |
| int k = LENGTH(tau), |
| n = INTEGER(coerceVector(getAttrib(qr, R_DimSymbol), INTSXP))[0], |
| *Bdims = INTEGER(coerceVector(getAttrib(B, R_DimSymbol), INTSXP)); |
| if(Bdims[0] != n) |
| error(_("right-hand side should have %d not %d rows"), n, Bdims[0]); |
| int nrhs = Bdims[1], |
| lwork = -1, info; |
| double tmp; |
| F77_CALL(dormqr)("L", "T", &n, &nrhs, &k, |
| REAL(qr), &n, REAL(tau), REAL(B), &n, |
| &tmp, &lwork, &info FCONE FCONE); |
| if (info != 0) |
| error(_("error code %d from Lapack routine '%s'"), info, "dormqr [tmp]"); |
| lwork = (int) tmp; |
| double *work = (double *) R_alloc(lwork, sizeof(double)); |
| F77_CALL(dormqr)("L", "T", &n, &nrhs, &k, |
| REAL(qr), &n, REAL(tau), REAL(B), &n, |
| work, &lwork, &info FCONE FCONE); |
| if (info != 0) |
| error(_("error code %d from Lapack routine '%s'"), info, "dormqr [work]"); |
| F77_CALL(dtrtrs)("U", "N", "N", &k, &nrhs, |
| REAL(qr), &n, REAL(B), &n, &info |
| FCONE FCONE FCONE); |
| if (info != 0) |
| error(_("error code %d from Lapack routine '%s'"), info, "dtrtrs"); |
| UNPROTECT(1); |
| return B; |
| } |
| |
| /* Real case of qr.qy and qr.aty */ |
| static SEXP qr_qy_real(SEXP Q, SEXP Bin, SEXP trans) |
| { |
| int n, nrhs, k, *Bdims, tr; |
| SEXP B, qr = VECTOR_ELT(Q, 0), tau = VECTOR_ELT(Q, 2); |
| double *work, tmp; |
| |
| k = LENGTH(tau); |
| if (!isMatrix(Bin)) error(_("'b' must be a numeric matrix")); |
| tr = asLogical(trans); |
| if(tr == NA_LOGICAL) error(_("invalid '%s' argument"), "trans"); |
| |
| B = PROTECT(isReal(Bin) ? duplicate(Bin) : coerceVector(Bin, REALSXP)); |
| n = INTEGER(coerceVector(getAttrib(qr, R_DimSymbol), INTSXP))[0]; |
| Bdims = INTEGER(coerceVector(getAttrib(B, R_DimSymbol), INTSXP)); |
| if(Bdims[0] != n) |
| error(_("right-hand side should have %d not %d rows"), n, Bdims[0]); |
| nrhs = Bdims[1]; |
| int lwork = -1, info; |
| F77_CALL(dormqr)("L", tr ? "T" : "N", &n, &nrhs, &k, |
| REAL(qr), &n, REAL(tau), REAL(B), &n, |
| &tmp, &lwork, &info FCONE FCONE); |
| if (info != 0) |
| error(_("error code %d from Lapack routine '%s'"), info, "dormqr"); |
| lwork = (int) tmp; |
| work = (double *) R_alloc(lwork, sizeof(double)); |
| F77_CALL(dormqr)("L", tr ? "T" : "N", &n, &nrhs, &k, |
| REAL(qr), &n, REAL(tau), REAL(B), &n, |
| work, &lwork, &info FCONE FCONE); |
| if (info != 0) |
| error(_("error code %d from Lapack routine '%s'"), info, "dormqr"); |
| UNPROTECT(1); |
| return B; |
| } |
| |
| /* TODO : add a *complex* version, using LAPACK ZGETRF() */ |
| static SEXP det_ge_real(SEXP Ain, SEXP logarithm) |
| { |
| int info, sign = 1, useLog = asLogical(logarithm); |
| double modulus = 0.0; /* -Wall */ |
| |
| if (!isMatrix(Ain)) error(_("'a' must be a numeric matrix")); |
| if (useLog == NA_LOGICAL) error(_("argument 'logarithm' must be logical")); |
| SEXP A = PROTECT(isReal(Ain) ? duplicate(Ain): coerceVector(Ain, REALSXP)); |
| int *Adims = INTEGER(coerceVector(getAttrib(Ain, R_DimSymbol), INTSXP)); |
| int n = Adims[0]; |
| if (Adims[1] != n) error(_("'a' must be a square matrix")); |
| int *jpvt = (int *) R_alloc(n, sizeof(int)); |
| F77_CALL(dgetrf)(&n, &n, REAL(A), &n, jpvt, &info); |
| if (info < 0) |
| error(_("error code %d from Lapack routine '%s'"), info, "dgetrf"); |
| else if (info > 0) { /* Singular matrix: U[i,i] (i := info) is 0 */ |
| /*warning("Lapack dgetrf(): singular matrix: U[%d,%d]=0", info,info);*/ |
| modulus = |
| #ifdef _not_quite_/* unfortunately does not work -- FIXME */ |
| (ISNAN(REAL(A)[(info-1)*n + (info-1)])) /* pivot is NA/NaN */ |
| ? R_NaN : |
| #endif |
| (useLog ? R_NegInf : 0.); |
| } |
| else { |
| for (int i = 0; i < n; i++) if (jpvt[i] != (i + 1)) sign = -sign; |
| if (useLog) { |
| modulus = 0.0; |
| size_t N1 = n+1; |
| for (int i = 0; i < n; i++) { |
| double dii = REAL(A)[i * N1]; /* ith diagonal element */ |
| modulus += log(dii < 0 ? -dii : dii); |
| if (dii < 0) sign = -sign; |
| } |
| } else { |
| modulus = 1.0; |
| size_t N1 = n+1; |
| for (int i = 0; i < n; i++) modulus *= REAL(A)[i * N1]; |
| if (modulus < 0) { |
| modulus = -modulus; |
| sign = -sign; |
| } |
| } |
| } |
| SEXP val = PROTECT(allocVector(VECSXP, 2)); |
| SEXP nm = PROTECT(allocVector(STRSXP, 2)); |
| SET_STRING_ELT(nm, 0, mkChar("modulus")); |
| SET_STRING_ELT(nm, 1, mkChar("sign")); |
| setAttrib(val, R_NamesSymbol, nm); |
| SET_VECTOR_ELT(val, 0, ScalarReal(modulus)); |
| SEXP s_logarithm = install("logarithm"); |
| setAttrib(VECTOR_ELT(val, 0), s_logarithm, ScalarLogical(useLog)); |
| SET_VECTOR_ELT(val, 1, ScalarInteger(sign)); |
| setAttrib(val, R_ClassSymbol, ScalarString(mkChar("det"))); |
| UNPROTECT(3); |
| return val; |
| } |
| |
| static SEXP mod_do_lapack(SEXP call, SEXP op, SEXP args, SEXP env) |
| { |
| SEXP ans = R_NilValue; |
| |
| switch(PRIMVAL(op)) { |
| case 0: ans = La_qr_cmplx(CAR(args)); break; |
| case 1: ans = La_rs(CAR(args), CADR(args)); break; |
| case 2: ans = La_rs_cmplx(CAR(args), CADR(args)); break; |
| case 3: ans = La_rg(CAR(args), CADR(args)); break; |
| case 41: ans = La_rg_cmplx(CAR(args), CADR(args)); break; |
| case 5: ans = La_rs(CAR(args), CADR(args)); break; |
| case 51: ans = La_rs_cmplx(CAR(args), CADR(args)); break; |
| case 6: ans = La_dlange(CAR(args), CADR(args)); break; |
| case 7: ans = La_dgecon(CAR(args), CADR(args)); break; |
| case 8: ans = La_dtrcon(CAR(args), CADR(args)); break; |
| case 9: ans = La_zgecon(CAR(args), CADR(args)); break; |
| case 10: ans = La_ztrcon(CAR(args), CADR(args)); break; |
| case 11: ans = La_solve_cmplx(CAR(args), CADR(args)); break; |
| |
| case 100: ans = La_solve(CAR(args), CADR(args), CADDR(args)); break; |
| case 101: ans = La_qr(CAR(args)); break; |
| |
| case 200: ans = La_chol(CAR(args), CADR(args), CADDR(args)); break; |
| case 201: ans = La_chol2inv(CAR(args), CADR(args)); break; |
| |
| case 300: ans = qr_coef_real(CAR(args), CADR(args)); break; |
| case 301: ans = qr_qy_real(CAR(args), CADR(args), CADDR(args)); break; |
| case 302: ans = det_ge_real(CAR(args), CADR(args)); break; |
| case 303: ans = qr_coef_cmplx(CAR(args), CADR(args)); break; |
| case 304: ans = qr_qy_cmplx(CAR(args), CADR(args), CADDR(args)); break; |
| |
| case 400: |
| { |
| SEXP a1, a2, a3, a4; |
| a1 = CAR(args); args = CDR(args); |
| a2 = CAR(args); args = CDR(args); |
| a3 = CAR(args); args = CDR(args); |
| a4 = CAR(args); args = CDR(args); |
| ans = La_svd(a1, a2, a3, a4, CAR(args)); |
| break; |
| } |
| case 401: |
| { |
| SEXP a1, a2, a3, a4; |
| a1 = CAR(args); args = CDR(args); |
| a2 = CAR(args); args = CDR(args); |
| a3 = CAR(args); args = CDR(args); |
| a4 = CAR(args); args = CDR(args); |
| ans = La_svd_cmplx(a1, a2, a3, a4, CAR(args)); |
| break; |
| } |
| case 1000: |
| { |
| int major, minor, patch; |
| char str[20]; |
| F77_CALL(ilaver)(&major, &minor, &patch); |
| snprintf(str, 20, "%d.%d.%d", major, minor, patch); |
| ans = mkString(str); |
| break; |
| } |
| case 1001: |
| { |
| #if defined(HAVE_DLADDR) && defined(HAVE_REALPATH) |
| Dl_info dl_info; |
| // dladdr has first arg void * on Solaris. This is not POSIX. |
| if (dladdr((void *) F77_NAME(ilaver), &dl_info)) { |
| char buf[PATH_MAX+1]; |
| char *res = realpath(dl_info.dli_fname, buf); |
| if (res) { |
| ans = mkString(res); |
| break; |
| } |
| } |
| #endif |
| ans = mkString(""); /* LAPACK library not known */ |
| break; |
| } |
| |
| }// switch |
| |
| return ans; |
| } |
| |
| |
| /* ------------------------------------------------------------ */ |
| |
| #include <Rmodules/Rlapack.h> |
| #include <R_ext/Rdynload.h> |
| |
| void |
| R_init_lapack(DllInfo *info) |
| { |
| R_LapackRoutines *tmp; |
| tmp = Calloc(1, R_LapackRoutines); |
| |
| tmp->do_lapack = mod_do_lapack; |
| R_setLapackRoutines(tmp); |
| } |