diff -Nru orig/PRIMME/src/Makevars patched/PRIMME/src/Makevars --- orig/PRIMME/src/Makevars 2023-05-08 18:05:53.000000000 +0200 +++ patched/PRIMME/src/Makevars 2024-02-08 10:45:35.959046600 +0100 @@ -1,5 +1,9 @@ PKG_CXXFLAGS = -I../inst/include -DPRIMME_INT_SIZE=0 -DF77UNDERSCORE -DUSE_XHEEV -DUSE_ZGESV -DUSE_XHEGV -DPRIMME_INT_SIZE=0 -DPRIMME_WITHOUT_FLOAT -DPRIMME_BLAS_RCOMPLEX -# Linker will discard primmeext if R provides a full LAPACK +# On some platorms, linker will discard primmeext if R provides a full LAPACK. +# +# Note ld.lld on Windows, however, will fail when a symbol from Rlapack (DLL) +# is later found in libprimmeext.a. Functions in Rlapack shouldn't hence be +# in libprimmeext.a PKG_LIBS = primme/libprimme.a $(LAPACK_LIBS) primme/libprimmeext.a $(BLAS_LIBS) $(FLIBS) $(SHLIB): primme/libprimme.a primme/libprimmeext.a diff -Nru orig/PRIMME/src/primme/zhegv.c patched/PRIMME/src/primme/zhegv.c --- orig/PRIMME/src/primme/zhegv.c 2023-07-29 14:58:47.000000000 +0200 +++ patched/PRIMME/src/primme/zhegv.c 2024-02-08 10:29:05.191430300 +0100 @@ -397,258 +397,6 @@ static real c_b1094 = 0.f; static real c_b1095 = 1.f; -/* Subroutine */ int zheev_(char *jobz, char *uplo, integer *n, doublecomplex - *a, integer *lda, doublereal *w, doublecomplex *work, integer *lwork, - doublereal *rwork, integer *info) -{ - /* System generated locals */ - integer a_dim1, a_offset, i__1, i__2; - doublereal d__1; - - /* Local variables */ - integer nb; - doublereal eps; - integer inde; - doublereal anrm; - integer imax; - doublereal rmin, rmax; - extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *, - integer *); - doublereal sigma; - extern logical lsame_(char *, char *); - integer iinfo; - logical lower, wantz; - extern doublereal dlamch_(char *); - integer iscale; - doublereal safmin; - extern integer ilaenv_(integer *, char *, char *, integer *, integer *, - integer *, integer *, ftnlen, ftnlen); - extern /* Subroutine */ int xerbla_(char *, integer *); - doublereal bignum; - extern doublereal zlanhe_(char *, char *, integer *, doublecomplex *, - integer *, doublereal *); - integer indtau; - extern /* Subroutine */ int dsterf_(integer *, doublereal *, doublereal *, - integer *), zlascl_(char *, integer *, integer *, doublereal *, - doublereal *, integer *, integer *, doublecomplex *, integer *, - integer *); - integer indwrk; - extern /* Subroutine */ int zhetrd_(char *, integer *, doublecomplex *, - integer *, doublereal *, doublereal *, doublecomplex *, - doublecomplex *, integer *, integer *); - integer llwork; - doublereal smlnum; - integer lwkopt; - logical lquery; - extern /* Subroutine */ int zsteqr_(char *, integer *, doublereal *, - doublereal *, doublecomplex *, integer *, doublereal *, integer *), zungtr_(char *, integer *, doublecomplex *, integer *, - doublecomplex *, doublecomplex *, integer *, integer *); - - -/* -- LAPACK driver routine (version 3.2) -- */ -/* -- LAPACK is a software package provided by Univ. of Tennessee, -- */ -/* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */ -/* November 2006 */ - - -/* Purpose */ -/* ======= */ - -/* ZHEEV computes all eigenvalues and, optionally, eigenvectors of a */ -/* complex Hermitian matrix A. */ - -/* Arguments */ -/* ========= */ - -/* JOBZ (input) CHARACTER*1 */ -/* = 'N': Compute eigenvalues only; */ -/* = 'V': Compute eigenvalues and eigenvectors. */ - -/* UPLO (input) CHARACTER*1 */ -/* = 'U': Upper triangle of A is stored; */ -/* = 'L': Lower triangle of A is stored. */ - -/* N (input) INTEGER */ -/* The order of the matrix A. N >= 0. */ - -/* A (input/output) COMPLEX*16 array, dimension (LDA, N) */ -/* On entry, the Hermitian matrix A. If UPLO = 'U', the */ -/* leading N-by-N upper triangular part of A contains the */ -/* upper triangular part of the matrix A. If UPLO = 'L', */ -/* the leading N-by-N lower triangular part of A contains */ -/* the lower triangular part of the matrix A. */ -/* On exit, if JOBZ = 'V', then if INFO = 0, A contains the */ -/* orthonormal eigenvectors of the matrix A. */ -/* If JOBZ = 'N', then on exit the lower triangle (if UPLO='L') */ -/* or the upper triangle (if UPLO='U') of A, including the */ -/* diagonal, is destroyed. */ - -/* LDA (input) INTEGER */ -/* The leading dimension of the array A. LDA >= f2cmax(1,N). */ - -/* W (output) DOUBLE PRECISION array, dimension (N) */ -/* If INFO = 0, the eigenvalues in ascending order. */ - -/* WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK)) */ -/* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */ - -/* LWORK (input) INTEGER */ -/* The length of the array WORK. LWORK >= f2cmax(1,2*N-1). */ -/* For optimal efficiency, LWORK >= (NB+1)*N, */ -/* where NB is the blocksize for ZHETRD returned by ILAENV. */ - -/* If LWORK = -1, then a workspace query is assumed; the routine */ -/* only calculates the optimal size of the WORK array, returns */ -/* this value as the first entry of the WORK array, and no error */ -/* message related to LWORK is issued by XERBLA. */ - -/* RWORK (workspace) DOUBLE PRECISION array, dimension (f2cmax(1, 3*N-2)) */ - -/* INFO (output) INTEGER */ -/* = 0: successful exit */ -/* < 0: if INFO = -i, the i-th argument had an illegal value */ -/* > 0: if INFO = i, the algorithm failed to converge; i */ -/* off-diagonal elements of an intermediate tridiagonal */ -/* form did not converge to zero. */ - -/* ===================================================================== */ - - -/* Test the input parameters. */ - - /* Parameter adjustments */ - a_dim1 = *lda; - a_offset = 1 + a_dim1; - a -= a_offset; - --w; - --work; - --rwork; - - /* Function Body */ - wantz = lsame_(jobz, "V"); - lower = lsame_(uplo, "L"); - lquery = *lwork == -1; - - *info = 0; - if (! (wantz || lsame_(jobz, "N"))) { - *info = -1; - } else if (! (lower || lsame_(uplo, "U"))) { - *info = -2; - } else if (*n < 0) { - *info = -3; - } else if (*lda < f2cmax(1,*n)) { - *info = -5; - } - - if (*info == 0) { - nb = ilaenv_(&c__1, "ZHETRD", uplo, n, &c_n1, &c_n1, &c_n1, (ftnlen)6, - (ftnlen)1); -/* Computing MAX */ - i__1 = 1, i__2 = (nb + 1) * *n; - lwkopt = f2cmax(i__1,i__2); - work[1].r = (doublereal) lwkopt, work[1].i = 0.; - -/* Computing MAX */ - i__1 = 1, i__2 = (*n << 1) - 1; - if (*lwork < f2cmax(i__1,i__2) && ! lquery) { - *info = -8; - } - } - - if (*info != 0) { - i__1 = -(*info); - xerbla_("ZHEEV ", &i__1); - return 0; - } else if (lquery) { - return 0; - } - -/* Quick return if possible */ - - if (*n == 0) { - return 0; - } - - if (*n == 1) { - i__1 = a_dim1 + 1; - w[1] = a[i__1].r; - work[1].r = 1., work[1].i = 0.; - if (wantz) { - i__1 = a_dim1 + 1; - a[i__1].r = 1., a[i__1].i = 0.; - } - return 0; - } - -/* Get machine constants. */ - - safmin = dlamch_("Safe minimum"); - eps = dlamch_("Precision"); - smlnum = safmin / eps; - bignum = 1. / smlnum; - rmin = sqrt(smlnum); - rmax = sqrt(bignum); - -/* Scale matrix to allowable range, if necessary. */ - - anrm = zlanhe_("M", uplo, n, &a[a_offset], lda, &rwork[1]); - iscale = 0; - if (anrm > 0. && anrm < rmin) { - iscale = 1; - sigma = rmin / anrm; - } else if (anrm > rmax) { - iscale = 1; - sigma = rmax / anrm; - } - if (iscale == 1) { - zlascl_(uplo, &c__0, &c__0, &c_b18, &sigma, n, n, &a[a_offset], lda, - info); - } - -/* Call ZHETRD to reduce Hermitian matrix to tridiagonal form. */ - - inde = 1; - indtau = 1; - indwrk = indtau + *n; - llwork = *lwork - indwrk + 1; - zhetrd_(uplo, n, &a[a_offset], lda, &w[1], &rwork[inde], &work[indtau], & - work[indwrk], &llwork, &iinfo); - -/* For eigenvalues only, call DSTERF. For eigenvectors, first call */ -/* ZUNGTR to generate the unitary matrix, then call ZSTEQR. */ - - if (! wantz) { - dsterf_(n, &w[1], &rwork[inde], info); - } else { - zungtr_(uplo, n, &a[a_offset], lda, &work[indtau], &work[indwrk], & - llwork, &iinfo); - indwrk = inde + *n; - zsteqr_(jobz, n, &w[1], &rwork[inde], &a[a_offset], lda, &rwork[ - indwrk], info); - } - -/* If matrix was scaled, then rescale eigenvalues appropriately. */ - - if (iscale == 1) { - if (*info == 0) { - imax = *n; - } else { - imax = *info - 1; - } - d__1 = 1. / sigma; - dscal_(&imax, &d__1, &w[1], &c__1); - } - -/* Set WORK(1) to optimal complex workspace size. */ - - work[1].r = (doublereal) lwkopt, work[1].i = 0.; - - return 0; - -/* End of ZHEEV */ - -} /* zheev_ */ - /* Subroutine */ int zhegs2_(integer *itype, char *uplo, integer *n, doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb, integer *info) @@ -5803,219 +5551,6 @@ } /* zpotf2_ */ -/* Subroutine */ int zpotrf_(char *uplo, integer *n, doublecomplex *a, - integer *lda, integer *info) -{ - /* System generated locals */ - integer a_dim1, a_offset, i__1, i__2, i__3, i__4; - doublecomplex z__1; - - /* Local variables */ - integer j, jb, nb; - extern logical lsame_(char *, char *); - extern /* Subroutine */ int zgemm_(char *, char *, integer *, integer *, - integer *, doublecomplex *, doublecomplex *, integer *, - doublecomplex *, integer *, doublecomplex *, doublecomplex *, - integer *), zherk_(char *, char *, integer *, - integer *, doublereal *, doublecomplex *, integer *, doublereal *, - doublecomplex *, integer *); - logical upper; - extern /* Subroutine */ int ztrsm_(char *, char *, char *, char *, - integer *, integer *, doublecomplex *, doublecomplex *, integer *, - doublecomplex *, integer *), - zpotf2_(char *, integer *, doublecomplex *, integer *, integer *), xerbla_(char *, integer *); - extern integer ilaenv_(integer *, char *, char *, integer *, integer *, - integer *, integer *, ftnlen, ftnlen); - - -/* -- LAPACK routine (version 3.3.1) -- */ -/* -- LAPACK is a software package provided by Univ. of Tennessee, -- */ -/* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */ -/* -- April 2011 -- */ - - -/* Purpose */ -/* ======= */ - -/* ZPOTRF computes the Cholesky factorization of a complex Hermitian */ -/* positive definite matrix A. */ - -/* The factorization has the form */ -/* A = U**H * U, if UPLO = 'U', or */ -/* A = L * L**H, if UPLO = 'L', */ -/* where U is an upper triangular matrix and L is lower triangular. */ - -/* This is the block version of the algorithm, calling Level 3 BLAS. */ - -/* Arguments */ -/* ========= */ - -/* UPLO (input) CHARACTER*1 */ -/* = 'U': Upper triangle of A is stored; */ -/* = 'L': Lower triangle of A is stored. */ - -/* N (input) INTEGER */ -/* The order of the matrix A. N >= 0. */ - -/* A (input/output) COMPLEX*16 array, dimension (LDA,N) */ -/* On entry, the Hermitian matrix A. If UPLO = 'U', the leading */ -/* N-by-N upper triangular part of A contains the upper */ -/* triangular part of the matrix A, and the strictly lower */ -/* triangular part of A is not referenced. If UPLO = 'L', the */ -/* leading N-by-N lower triangular part of A contains the lower */ -/* triangular part of the matrix A, and the strictly upper */ -/* triangular part of A is not referenced. */ - -/* On exit, if INFO = 0, the factor U or L from the Cholesky */ -/* factorization A = U**H *U or A = L*L**H. */ - -/* LDA (input) INTEGER */ -/* The leading dimension of the array A. LDA >= f2cmax(1,N). */ - -/* INFO (output) INTEGER */ -/* = 0: successful exit */ -/* < 0: if INFO = -i, the i-th argument had an illegal value */ -/* > 0: if INFO = i, the leading minor of order i is not */ -/* positive definite, and the factorization could not be */ -/* completed. */ - -/* ===================================================================== */ - - -/* Test the input parameters. */ - - /* Parameter adjustments */ - a_dim1 = *lda; - a_offset = 1 + a_dim1; - a -= a_offset; - - /* Function Body */ - *info = 0; - upper = lsame_(uplo, "U"); - if (! upper && ! lsame_(uplo, "L")) { - *info = -1; - } else if (*n < 0) { - *info = -2; - } else if (*lda < f2cmax(1,*n)) { - *info = -4; - } - if (*info != 0) { - i__1 = -(*info); - xerbla_("ZPOTRF", &i__1); - return 0; - } - -/* Quick return if possible */ - - if (*n == 0) { - return 0; - } - -/* Determine the block size for this environment. */ - - nb = ilaenv_(&c__1, "ZPOTRF", uplo, n, &c_n1, &c_n1, &c_n1, (ftnlen)6, ( - ftnlen)1); - if (nb <= 1 || nb >= *n) { - -/* Use unblocked code. */ - - zpotf2_(uplo, n, &a[a_offset], lda, info); - } else { - -/* Use blocked code. */ - - if (upper) { - -/* Compute the Cholesky factorization A = U**H *U. */ - - i__1 = *n; - i__2 = nb; - for (j = 1; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) { - -/* Update and factorize the current diagonal block and test */ -/* for non-positive-definiteness. */ - -/* Computing MIN */ - i__3 = nb, i__4 = *n - j + 1; - jb = f2cmin(i__3,i__4); - i__3 = j - 1; - zherk_("Upper", "Conjugate transpose", &jb, &i__3, &c_b613, & - a[j * a_dim1 + 1], lda, &c_b18, &a[j + j * a_dim1], - lda); - zpotf2_("Upper", &jb, &a[j + j * a_dim1], lda, info); - if (*info != 0) { - goto L30; - } - if (j + jb <= *n) { - -/* Compute the current block row. */ - - i__3 = *n - j - jb + 1; - i__4 = j - 1; - z__1.r = -1., z__1.i = -0.; - zgemm_("Conjugate transpose", "No transpose", &jb, &i__3, - &i__4, &z__1, &a[j * a_dim1 + 1], lda, &a[(j + jb) - * a_dim1 + 1], lda, &c_b1, &a[j + (j + jb) * - a_dim1], lda); - i__3 = *n - j - jb + 1; - ztrsm_("Left", "Upper", "Conjugate transpose", "Non-unit", - &jb, &i__3, &c_b1, &a[j + j * a_dim1], lda, &a[j - + (j + jb) * a_dim1], lda); - } -/* L10: */ - } - - } else { - -/* Compute the Cholesky factorization A = L*L**H. */ - - i__2 = *n; - i__1 = nb; - for (j = 1; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) { - -/* Update and factorize the current diagonal block and test */ -/* for non-positive-definiteness. */ - -/* Computing MIN */ - i__3 = nb, i__4 = *n - j + 1; - jb = f2cmin(i__3,i__4); - i__3 = j - 1; - zherk_("Lower", "No transpose", &jb, &i__3, &c_b613, &a[j + - a_dim1], lda, &c_b18, &a[j + j * a_dim1], lda); - zpotf2_("Lower", &jb, &a[j + j * a_dim1], lda, info); - if (*info != 0) { - goto L30; - } - if (j + jb <= *n) { - -/* Compute the current block column. */ - - i__3 = *n - j - jb + 1; - i__4 = j - 1; - z__1.r = -1., z__1.i = -0.; - zgemm_("No transpose", "Conjugate transpose", &i__3, &jb, - &i__4, &z__1, &a[j + jb + a_dim1], lda, &a[j + - a_dim1], lda, &c_b1, &a[j + jb + j * a_dim1], lda); - i__3 = *n - j - jb + 1; - ztrsm_("Right", "Lower", "Conjugate transpose", "Non-unit" - , &i__3, &jb, &c_b1, &a[j + j * a_dim1], lda, &a[ - j + jb + j * a_dim1], lda); - } -/* L20: */ - } - } - } - goto L40; - -L30: - *info = *info + j - 1; - -L40: - return 0; - -/* End of ZPOTRF */ - -} /* zpotrf_ */ /* Subroutine */ int zsteqr_(char *compz, integer *n, doublereal *d__, doublereal *e, doublecomplex *z__, integer *ldz, doublereal *work,