/* * R : A Computer Language for Statistical Data Analysis * Copyright (C) 2001--2023 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 #endif #include #include /* for toupper */ #include /* for PATH_MAX */ #include /* for realpath */ #include /* for strcpy */ #ifdef HAVE_UNISTD_H # include /* for realpath on some systems */ #endif #ifdef HAVE_DLFCN_H # include /* 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' */ } /* Lapack valid 'uplo' (upper/lower) for triangular/symmetric matrices: */ static char La_valid_uplo(const char *uplostr) { if (strlen(uplostr) != 1) error(_("argument type[1]='%s' must be a character string of string length 1"), uplostr); char uplo = (char) toupper(*uplostr); if (uplo != 'U' && uplo != 'L') error(_("argument type[1]='%s' must be 'U' or 'L'"), uplostr); return uplo; /* 'U' or 'L' */ } /* La.svd, called from svd */ static SEXP La_svd(SEXP jobu, SEXP x, SEXP s, SEXP u, SEXP vt) { if (!isString(jobu)) error(_("'%s' must be a character string"), "jobu"); 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(_("'%s' must be a numeric matrix"), "A"); if (!isString(type)) error(_("'%s' must be a character string"), "type"); 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(_("'%s' must be a numeric matrix"), "A"); if (!isString(norm)) error(_("'%s' must be a character string"), "norm"); 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 code %d from Lapack routine '%s'"), info, "dgetrf()"); } 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 code %d from Lapack routine '%s'"), info, "dgecon()"); return val; } /* Real case of kappa.tri (and also rcond for a triangular matrix), both * uplo = "U" or uplo = "L" */ static SEXP La_dtrcon3(SEXP A, SEXP norm, SEXP uplo) { int *xdims, n, nprot = 0, info; char typNorm[] = {'\0', '\0'}; char uploC[] = {'\0', '\0'}; if (!isMatrix(A)) error(_("'%s' must be a numeric matrix"), "A"); if (!isString(norm)) error(_("'%s' must be a character string"), "norm"); if (!isString(uplo)) error(_("'%s' must be a character string"), "uplo"); 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))); uploC [0] = La_valid_uplo(CHAR(asChar(uplo))); nprot++; SEXP val = PROTECT(allocVector(REALSXP, 1)); F77_CALL(dtrcon)(typNorm, uploC, "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; } /* Real case of kappa.tri (and also rcond for a triangular matrix) */ /* hardcoded uplo = "U" -- using only upper triangular entries */ static SEXP La_dtrcon(SEXP A, SEXP norm) { int *xdims, n, nprot = 0, info; char typNorm[] = {'\0', '\0'}; if (!isMatrix(A)) error(_("'%s' must be a numeric matrix"), "A"); if (!isString(norm)) error(_("'%s' must be a character string"), "norm"); 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 code %d from Lapack routine '%s'"), info, "dtrcon()"); return val; } /* norm() except for 2-norm */ static SEXP La_zlange(SEXP A, SEXP type) { #ifdef HAVE_FORTRAN_DOUBLE_COMPLEX int *xdims, m, n; double *work = NULL; char typNorm[] = {'\0', '\0'}; if (!(isMatrix(A) && isComplex(A))) error(_("'%s' must be a complex matrix"), "A"); if (!isString(type)) error(_("'%s' must be a character string"), "type"); 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((size_t)m, sizeof(Rcomplex)); REAL(val)[0] = F77_CALL(zlange)(typNorm, &m, &n, COMPLEX(A), &m, work FCONE); UNPROTECT(1); return val; #else error(_("Fortran complex functions are not available on this platform")); return R_NilValue; /* -Wall */ #endif } /* 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(_("'%s' must be a character string"), "norm"); if (!(isMatrix(A) && isComplex(A))) error(_("'%s' must be a complex matrix"), "A"); 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 code %d from Lapack routine '%s'"), info, "zgetrf()"); } else { REAL(val)[0] = 0.; /* rcond = 0 <==> singularity */ UNPROTECT(1); return val; } } 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 code %d from Lapack routine '%s'"), info, "zgecon()"); 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) */ /* hardcoded uplo = "U" -- using only upper triangular entries */ 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(_("'%s' must be a character string"), "norm"); if (!(isMatrix(A) && isComplex(A))) error(_("'%s' must be a complex matrix"), "A"); 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 code %d from Lapack routine '%s'"), info, "ztrcon()"); return val; #else error(_("Fortran complex functions are not available on this platform")); return R_NilValue; /* -Wall */ #endif } // La_ztrcon() /* Complex case of kappa.tri (and also rcond for a triangular matrix) */ static SEXP La_ztrcon3(SEXP A, SEXP norm, SEXP uplo) { #ifdef HAVE_FORTRAN_DOUBLE_COMPLEX SEXP val; int *dims, n, info; char typNorm[] = {'\0', '\0'}; char uploC[] = {'\0', '\0'}; if (!(isMatrix(A) && isComplex(A))) error(_("'%s' must be a complex matrix"), "A"); if (!isString(norm)) error(_("'%s' must be a character string"), "norm"); if (!isString(uplo)) error(_("'%s' must be a character string"), "uplo"); 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))); uploC [0] = La_valid_uplo(CHAR(asChar(uplo))); val = PROTECT(allocVector(REALSXP, 1)); F77_CALL(ztrcon)(typNorm, uploC, "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, SEXP tolin) { #ifdef HAVE_FORTRAN_DOUBLE_COMPLEX int n, p, info, *ipiv, *Adims, *Bdims; Rcomplex *avals; SEXP B, Adn, Bdn; if (!isMatrix(A)) error(_("'%s' must be a complex matrix"), "a"); Adims = INTEGER(coerceVector(getAttrib(A, R_DimSymbol), INTSXP)); n = Adims[0]; if(n == 0) error(_("'a' is 0-diml")); size_t nl = n; 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), nl * 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(nl * nl, sizeof(Rcomplex)); Memcpy(avals, COMPLEX(A), nl * nl); } 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 %s: system is exactly singular: U[%d,%d] = 0"), "zgesv", info, info); int OK = 1; for (size_t i = 0; i < nl*nl; i++) if (!isfinite(avals[i].r) || !isfinite(avals[i].i)) {OK = 0; break;} double tol = asReal(tolin); if(OK == 1 && tol > 0) { char one[2] = "1"; double anorm = F77_CALL(zlange)(one, &n, &n, COMPLEX(A), &n, (double*) NULL FCONE); Rcomplex *work = (Rcomplex *) R_alloc(2*nl, sizeof(Rcomplex)); double *rwork = (double *) R_alloc(2*nl, sizeof(double)); double rcond; F77_CALL(zgecon)(one, &n, avals, &n, &anorm, &rcond, work, rwork, &info FCONE); if (rcond < tol) error(_("system is computationally singular: reciprocal condition number = %g"), rcond); } 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(_("'%s' must be a complex matrix"), "a"); 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(_("'%s' must be a complex matrix"), "b"); 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(_("'%s' must be a complex matrix"), "b"); 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(_("'%s' must be a character string"), "jobu"); 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(_("'%s' must be a numeric matrix"), "a"); 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"), 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 not positive definite")); 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(_("'%s' must be a numeric matrix"), "a"); 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, tol = asReal(tolin); SEXP B, Adn, Bdn; if (!(isMatrix(A) && (TYPEOF(A) == REALSXP || TYPEOF(A) == INTSXP || TYPEOF(A) == LGLSXP))) error(_("'%s' must be a numeric matrix"), "a"); int *Adims = INTEGER(coerceVector(getAttrib(A, R_DimSymbol), INTSXP)); n = Adims[0]; if(n == 0) error(_("'a' is 0-diml")); size_t nl = n; 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), nl * 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(nl * nl, sizeof(double)); Memcpy(avals, REAL(A), nl * nl); } 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); // LAPACK 3.11.0 fails here if A contains NaNs int OK = 1; for (size_t i = 0; i < nl*nl; i++) if (!isfinite(avals[i])) {OK = 0; break;} if(OK == 1 && tol > 0) { char one[2] = "1"; double rcond; double anorm = F77_CALL(dlange)(one, &n, &n, REAL(A), &n, (double*) NULL FCONE); double *work = (double *) R_alloc(4*nl, 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(_("'%s' must be a numeric matrix"), "a"); 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)); // qraux 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(_("'%s' must be a numeric matrix"), "b"); 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(_("'%s' must be a numeric matrix"), "b"); 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(_("'%s' must be a numeric matrix"), "a"); 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 61: ans = La_zlange(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 81: ans = La_dtrcon3(CAR(args), CADR(args), CADDR(args)); break; case 9: ans = La_zgecon(CAR(args), CADR(args)); break; case 10: ans = La_ztrcon(CAR(args), CADR(args)); break; case 13: ans = La_ztrcon3(CAR(args), CADR(args), CADDR(args)); break; case 11: ans = La_solve_cmplx(CAR(args), CADR(args), CADDR(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: // La_version() { 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: // La_library() { #if defined(HAVE_DLADDR) && defined(HAVE_REALPATH) Dl_info dl_info; /* the call to dladdr() converts a function pointer to an object pointer, which is not allowed by ISO C, but there is no compliant alternative to using dladdr() */ // dladdr has first arg void * on Solaris. This is not POSIX. #ifdef __clang__ # pragma clang diagnostic push # pragma clang diagnostic ignored "-Wpedantic" #elif defined __GNUC__ # pragma GCC diagnostic push # pragma GCC diagnostic ignored "-Wpedantic" #endif if (dladdr((void *) F77_NAME(ilaver), &dl_info)) { char buf[PATH_MAX+1]; char *res = realpath(dl_info.dli_fname, buf); if (res) { SEXP nfo = R_NilValue; if (strstr(res, "flexiblas")) nfo = R_flexiblas_info(); if (isNull(nfo)) nfo = mkChar(res); ans = ScalarString(nfo); break; } } #ifdef __clang__ # pragma clang diagnostic pop #elif defined __GNUC__ # pragma GCC diagnostic pop #endif #endif ans = mkString(""); /* LAPACK library not known */ break; } }// switch return ans; } /* ------------------------------------------------------------ */ #include #include void R_init_lapack(DllInfo *info) { R_LapackRoutines *tmp; tmp = R_Calloc(1, R_LapackRoutines); tmp->do_lapack = mod_do_lapack; R_setLapackRoutines(tmp); }