ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #3992 from akohlmey/linalg-lapack-3.12

Browse Source 
Update internal linear algebra library to LAPACK 3.12
This commit is contained in:
Axel Kohlmeyer
2023-11-24 19:16:02 -05:00
committed by GitHub
parent 239407c4eb 2a1af662d4
commit 10f207b986
15 changed files with 109 additions and 103 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
lib/linalg/README
View File

@ -11,7 +11,7 @@ resulting library will follow the Fortran binary conventions.
Note that this is an *incomplete* subset of full BLAS/LAPACK.
The files correspond to LAPACK version 3.11.0.
The files correspond to LAPACK version 3.12.0.
You should only need to build and use the library in this directory if
you want to build LAMMPS with one of the listed packages AND you do not

16
lib/linalg/dbdsqr.cpp
View File

@ -40,7 +40,7 @@ int dbdsqr_(char *uplo, integer *n, integer *ncvt, integer *nru, integer *ncc, d
integer oldll;
doublereal shift, sigmn, oldsn;
extern int dswap_(integer *, doublereal *, integer *, doublereal *, integer *);
doublereal sminl, sigmx;
doublereal sigmx;
logical lower;
extern int dlasq1_(integer *, doublereal *, doublereal *, doublereal *, integer *),
dlasv2_(doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *,
@ -141,7 +141,7 @@ int dbdsqr_(char *uplo, integer *n, integer *ncvt, integer *nru, integer *ncc, d
d__2 = smax, d__3 = (d__1 = e[i__], abs(d__1));
smax = max(d__2, d__3);
}
sminl = 0.;
smin = 0.;
if (tol >= 0.) {
sminoa = abs(d__[1]);
if (sminoa == 0.) {
@ -185,7 +185,6 @@ L60:
d__[m] = 0.;
}
smax = (d__1 = d__[m], abs(d__1));
smin = smax;
i__1 = m - 1;
for (lll = 1; lll <= i__1; ++lll) {
ll = m - lll;
@ -197,7 +196,6 @@ L60:
if (abse <= thresh) {
goto L80;
}
smin = min(smin, abss);
d__1 = max(smax, abss);
smax = max(d__1, abse);
}
@ -243,7 +241,7 @@ L90:
}
if (tol >= 0.) {
mu = (d__1 = d__[ll], abs(d__1));
sminl = mu;
smin = mu;
i__1 = m - 1;
for (lll = ll; lll <= i__1; ++lll) {
if ((d__1 = e[lll], abs(d__1)) <= tol * mu) {
@ -251,7 +249,7 @@ L90:
goto L60;
}
mu = (d__2 = d__[lll + 1], abs(d__2)) * (mu / (mu + (d__1 = e[lll], abs(d__1))));
sminl = min(sminl, mu);
smin = min(smin, mu);
}
}
} else {
@ -262,7 +260,7 @@ L90:
}
if (tol >= 0.) {
mu = (d__1 = d__[m], abs(d__1));
sminl = mu;
smin = mu;
i__1 = ll;
for (lll = m - 1; lll >= i__1; --lll) {
if ((d__1 = e[lll], abs(d__1)) <= tol * mu) {
@ -270,14 +268,14 @@ L90:
goto L60;
}
mu = (d__2 = d__[lll], abs(d__2)) * (mu / (mu + (d__1 = e[lll], abs(d__1))));
sminl = min(sminl, mu);
smin = min(smin, mu);
}
}
}
oldll = ll;
oldm = m;
d__1 = eps, d__2 = tol * .01;
if (tol >= 0. && *n * tol * (sminl / smax) <= max(d__1, d__2)) {
if (tol >= 0. && *n * tol * (smin / smax) <= max(d__1, d__2)) {
shift = 0.;
} else {
if (idir == 1) {

17
lib/linalg/dgecon.cpp
View File

@ -20,6 +20,7 @@ int dgecon_(char *norm, integer *n, doublereal *a, integer *lda, doublereal *ano
integer *);
extern doublereal dlamch_(char *, ftnlen);
extern integer idamax_(integer *, doublereal *, integer *);
extern logical disnan_(doublereal *);
extern int xerbla_(char *, integer *, ftnlen);
doublereal ainvnm;
extern int dlatrs_(char *, char *, char *, char *, integer *, doublereal *, integer *,
@ -27,12 +28,13 @@ int dgecon_(char *norm, integer *n, doublereal *a, integer *lda, doublereal *ano
ftnlen);
logical onenrm;
char normin[1];
doublereal smlnum;
doublereal smlnum, hugeval;
a_dim1 = *lda;
a_offset = 1 + a_dim1;
a -= a_offset;
--work;
--iwork;
hugeval = dlamch_((char *)"Overflow", (ftnlen)8);
*info = 0;
onenrm = *(unsigned char *)norm == '1' || lsame_(norm, (char *)"O", (ftnlen)1, (ftnlen)1);
if (!onenrm && !lsame_(norm, (char *)"I", (ftnlen)1, (ftnlen)1)) {
@ -55,6 +57,13 @@ int dgecon_(char *norm, integer *n, doublereal *a, integer *lda, doublereal *ano
return 0;
} else if (*anorm == 0.) {
return 0;
} else if (disnan_(anorm)) {
*rcond = *anorm;
*info = -5;
return 0;
} else if (*anorm > hugeval) {
*info = -5;
return 0;
}
smlnum = dlamch_((char *)"Safe minimum", (ftnlen)12);
ainvnm = 0.;
@ -92,6 +101,12 @@ L10:
}
if (ainvnm != 0.) {
*rcond = 1. / ainvnm / *anorm;
} else {
*info = 1;
return 0;
}
if (disnan_(rcond) || *rcond > hugeval) {
*info = 1;
}
L20:
return 0;

4
lib/linalg/dgelsd.cpp
View File

@ -19,8 +19,7 @@ int dgelsd_(integer *m, integer *n, integer *nrhs, doublereal *a, integer *lda,
integer itau, nlvl, iascl, ibscl;
doublereal sfmin;
integer minmn, maxmn, itaup, itauq, mnthr, nwork;
extern int dlabad_(doublereal *, doublereal *),
dgebrd_(integer *, integer *, doublereal *, integer *, doublereal *, doublereal *,
extern int dgebrd_(integer *, integer *, doublereal *, integer *, doublereal *, doublereal *,
doublereal *, doublereal *, doublereal *, integer *, integer *);
extern doublereal dlamch_(char *, ftnlen),
dlange_(char *, integer *, integer *, doublereal *, integer *, doublereal *, ftnlen);
@ -189,7 +188,6 @@ int dgelsd_(integer *m, integer *n, integer *nrhs, doublereal *a, integer *lda,
sfmin = dlamch_((char *)"S", (ftnlen)1);
smlnum = sfmin / eps;
bignum = 1. / smlnum;
dlabad_(&smlnum, &bignum);
anrm = dlange_((char *)"M", m, n, &a[a_offset], lda, &work[1], (ftnlen)1);
iascl = 0;
if (anrm > 0. && anrm < smlnum) {

10
lib/linalg/dgelss.cpp
View File

@ -30,8 +30,7 @@ int dgelss_(integer *m, integer *n, integer *nrhs, doublereal *a, integer *lda,
integer minmn;
extern int dcopy_(integer *, doublereal *, integer *, doublereal *, integer *);
integer maxmn, itaup, itauq, mnthr, iwork;
extern int dlabad_(doublereal *, doublereal *),
dgebrd_(integer *, integer *, doublereal *, integer *, doublereal *, doublereal *,
extern int dgebrd_(integer *, integer *, doublereal *, integer *, doublereal *, doublereal *,
doublereal *, doublereal *, doublereal *, integer *, integer *);
extern doublereal dlamch_(char *, ftnlen),
dlange_(char *, integer *, integer *, doublereal *, integer *, doublereal *, ftnlen);
@ -208,7 +207,6 @@ int dgelss_(integer *m, integer *n, integer *nrhs, doublereal *a, integer *lda,
sfmin = dlamch_((char *)"S", (ftnlen)1);
smlnum = sfmin / eps;
bignum = 1. / smlnum;
dlabad_(&smlnum, &bignum);
anrm = dlange_((char *)"M", m, n, &a[a_offset], lda, &work[1], (ftnlen)1);
iascl = 0;
if (anrm > 0. && anrm < smlnum) {
@ -300,7 +298,7 @@ int dgelss_(integer *m, integer *n, integer *nrhs, doublereal *a, integer *lda,
&c_b46, &work[1], n, (ftnlen)1, (ftnlen)1);
dlacpy_((char *)"G", n, &bl, &work[1], n, &b[i__ * b_dim1 + 1], ldb, (ftnlen)1);
}
} else {
} else if (*nrhs == 1) {
dgemv_((char *)"T", n, n, &c_b79, &a[a_offset], lda, &b[b_offset], &c__1, &c_b46, &work[1],
&c__1, (ftnlen)1);
dcopy_(n, &work[1], &c__1, &b[b_offset], &c__1);
@ -376,7 +374,7 @@ int dgelss_(integer *m, integer *n, integer *nrhs, doublereal *a, integer *lda,
ldb, &c_b46, &work[iwork], m, (ftnlen)1, (ftnlen)1);
dlacpy_((char *)"G", m, &bl, &work[iwork], m, &b[i__ * b_dim1 + 1], ldb, (ftnlen)1);
}
} else {
} else if (*nrhs == 1) {
dgemv_((char *)"T", m, m, &c_b79, &work[il], &ldwork, &b[b_dim1 + 1], &c__1, &c_b46,
&work[iwork], &c__1, (ftnlen)1);
dcopy_(m, &work[iwork], &c__1, &b[b_dim1 + 1], &c__1);
@ -438,7 +436,7 @@ int dgelss_(integer *m, integer *n, integer *nrhs, doublereal *a, integer *lda,
ldb, &c_b46, &work[1], n, (ftnlen)1, (ftnlen)1);
dlacpy_((char *)"F", n, &bl, &work[1], n, &b[i__ * b_dim1 + 1], ldb, (ftnlen)1);
}
} else {
} else if (*nrhs == 1) {
dgemv_((char *)"T", m, n, &c_b79, &a[a_offset], lda, &b[b_offset], &c__1, &c_b46, &work[1],
&c__1, (ftnlen)1);
dcopy_(n, &work[1], &c__1, &b[b_offset], &c__1);

5
lib/linalg/dlabad.cpp
View File

@ -4,11 +4,6 @@ extern "C" {
#include "lmp_f2c.h"
int dlabad_(doublereal *small, doublereal *large)
{
double d_lmp_lg10(doublereal *), sqrt(doublereal);
if (d_lmp_lg10(large) > 2e3) {
*small = sqrt(*small);
*large = sqrt(*large);
}
return 0;
}
#ifdef __cplusplus

16
lib/linalg/dlaed2.cpp
View File

@ -5,7 +5,7 @@ extern "C" {
static doublereal c_b3 = -1.;
static integer c__1 = 1;
int dlaed2_(integer *k, integer *n, integer *n1, doublereal *d__, doublereal *q, integer *ldq,
integer *indxq, doublereal *rho, doublereal *z__, doublereal *dlamda, doublereal *w,
integer *indxq, doublereal *rho, doublereal *z__, doublereal *dlambda, doublereal *w,
doublereal *q2, integer *indx, integer *indxc, integer *indxp, integer *coltyp,
integer *info)
{
@ -35,7 +35,7 @@ int dlaed2_(integer *k, integer *n, integer *n1, doublereal *d__, doublereal *q,
q -= q_offset;
--indxq;
--z__;
--dlamda;
--dlambda;
--w;
--q2;
--indx;
@ -75,9 +75,9 @@ int dlaed2_(integer *k, integer *n, integer *n1, doublereal *d__, doublereal *q,
}
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
dlamda[i__] = d__[indxq[i__]];
dlambda[i__] = d__[indxq[i__]];
}
dlamrg_(n1, &n2, &dlamda[1], &c__1, &c__1, &indxc[1]);
dlamrg_(n1, &n2, &dlambda[1], &c__1, &c__1, &indxc[1]);
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
indx[i__] = indxq[indxc[i__]];
@ -94,11 +94,11 @@ int dlaed2_(integer *k, integer *n, integer *n1, doublereal *d__, doublereal *q,
for (j = 1; j <= i__1; ++j) {
i__ = indx[j];
dcopy_(n, &q[i__ * q_dim1 + 1], &c__1, &q2[iq2], &c__1);
dlamda[j] = d__[i__];
dlambda[j] = d__[i__];
iq2 += *n;
}
dlacpy_((char *)"A", n, n, &q2[1], n, &q[q_offset], ldq, (ftnlen)1);
dcopy_(n, &dlamda[1], &c__1, &d__[1], &c__1);
dcopy_(n, &dlambda[1], &c__1, &d__[1], &c__1);
goto L190;
}
i__1 = *n1;
@ -176,7 +176,7 @@ L80:
pj = nj;
} else {
++(*k);
dlamda[*k] = d__[pj];
dlambda[*k] = d__[pj];
w[*k] = z__[pj];
indxp[*k] = pj;
pj = nj;
@ -185,7 +185,7 @@ L80:
goto L80;
L100:
++(*k);
dlamda[*k] = d__[pj];
dlambda[*k] = d__[pj];
w[*k] = z__[pj];
indxp[*k] = pj;
for (j = 1; j <= 4; ++j) {

33
lib/linalg/dlaed3.cpp
View File

@ -3,10 +3,10 @@ extern "C" {
#endif
#include "lmp_f2c.h"
static integer c__1 = 1;
static doublereal c_b22 = 1.;
static doublereal c_b23 = 0.;
static doublereal c_b21 = 1.;
static doublereal c_b22 = 0.;
int dlaed3_(integer *k, integer *n, integer *n1, doublereal *d__, doublereal *q, integer *ldq,
doublereal *rho, doublereal *dlamda, doublereal *q2, integer *indx, integer *ctot,
doublereal *rho, doublereal *dlambda, doublereal *q2, integer *indx, integer *ctot,
doublereal *w, doublereal *s, integer *info)
{
integer q_dim1, q_offset, i__1, i__2;
@ -20,10 +20,9 @@ int dlaed3_(integer *k, integer *n, integer *n1, doublereal *d__, doublereal *q,
ftnlen, ftnlen),
dcopy_(integer *, doublereal *, integer *, doublereal *, integer *),
dlaed4_(integer *, integer *, doublereal *, doublereal *, doublereal *, doublereal *,
doublereal *, integer *);
extern doublereal dlamc3_(doublereal *, doublereal *);
extern int dlacpy_(char *, integer *, integer *, doublereal *, integer *, doublereal *,
integer *, ftnlen),
doublereal *, integer *),
dlacpy_(char *, integer *, integer *, doublereal *, integer *, doublereal *, integer *,
ftnlen),
dlaset_(char *, integer *, integer *, doublereal *, doublereal *, doublereal *, integer *,
ftnlen),
xerbla_(char *, integer *, ftnlen);
@ -31,7 +30,7 @@ int dlaed3_(integer *k, integer *n, integer *n1, doublereal *d__, doublereal *q,
q_dim1 = *ldq;
q_offset = 1 + q_dim1;
q -= q_offset;
--dlamda;
--dlambda;
--q2;
--indx;
--ctot;
@ -54,12 +53,8 @@ int dlaed3_(integer *k, integer *n, integer *n1, doublereal *d__, doublereal *q,
return 0;
}
i__1 = *k;
for (i__ = 1; i__ <= i__1; ++i__) {
dlamda[i__] = dlamc3_(&dlamda[i__], &dlamda[i__]) - dlamda[i__];
}
i__1 = *k;
for (j = 1; j <= i__1; ++j) {
dlaed4_(k, &j, &dlamda[1], &w[1], &q[j * q_dim1 + 1], rho, &d__[j], info);
dlaed4_(k, &j, &dlambda[1], &w[1], &q[j * q_dim1 + 1], rho, &d__[j], info);
if (*info != 0) {
goto L120;
}
@ -86,11 +81,11 @@ int dlaed3_(integer *k, integer *n, integer *n1, doublereal *d__, doublereal *q,
for (j = 1; j <= i__1; ++j) {
i__2 = j - 1;
for (i__ = 1; i__ <= i__2; ++i__) {
w[i__] *= q[i__ + j * q_dim1] / (dlamda[i__] - dlamda[j]);
w[i__] *= q[i__ + j * q_dim1] / (dlambda[i__] - dlambda[j]);
}
i__2 = *k;
for (i__ = j + 1; i__ <= i__2; ++i__) {
w[i__] *= q[i__ + j * q_dim1] / (dlamda[i__] - dlamda[j]);
w[i__] *= q[i__ + j * q_dim1] / (dlambda[i__] - dlambda[j]);
}
}
i__1 = *k;
@ -118,17 +113,17 @@ L110:
dlacpy_((char *)"A", &n23, k, &q[ctot[1] + 1 + q_dim1], ldq, &s[1], &n23, (ftnlen)1);
iq2 = *n1 * n12 + 1;
if (n23 != 0) {
dgemm_((char *)"N", (char *)"N", &n2, k, &n23, &c_b22, &q2[iq2], &n2, &s[1], &n23, &c_b23,
dgemm_((char *)"N", (char *)"N", &n2, k, &n23, &c_b21, &q2[iq2], &n2, &s[1], &n23, &c_b22,
&q[*n1 + 1 + q_dim1], ldq, (ftnlen)1, (ftnlen)1);
} else {
dlaset_((char *)"A", &n2, k, &c_b23, &c_b23, &q[*n1 + 1 + q_dim1], ldq, (ftnlen)1);
dlaset_((char *)"A", &n2, k, &c_b22, &c_b22, &q[*n1 + 1 + q_dim1], ldq, (ftnlen)1);
}
dlacpy_((char *)"A", &n12, k, &q[q_offset], ldq, &s[1], &n12, (ftnlen)1);
if (n12 != 0) {
dgemm_((char *)"N", (char *)"N", n1, k, &n12, &c_b22, &q2[1], n1, &s[1], &n12, &c_b23, &q[q_offset], ldq,
dgemm_((char *)"N", (char *)"N", n1, k, &n12, &c_b21, &q2[1], n1, &s[1], &n12, &c_b22, &q[q_offset], ldq,
(ftnlen)1, (ftnlen)1);
} else {
dlaset_((char *)"A", n1, k, &c_b23, &c_b23, &q[q_dim1 + 1], ldq, (ftnlen)1);
dlaset_((char *)"A", n1, k, &c_b22, &c_b22, &q[q_dim1 + 1], ldq, (ftnlen)1);
}
L120:
return 0;

22
lib/linalg/dlaed8.cpp
View File

@ -6,7 +6,7 @@ static doublereal c_b3 = -1.;
static integer c__1 = 1;
int dlaed8_(integer *icompq, integer *k, integer *n, integer *qsiz, doublereal *d__, doublereal *q,
integer *ldq, integer *indxq, doublereal *rho, integer *cutpnt, doublereal *z__,
doublereal *dlamda, doublereal *q2, integer *ldq2, doublereal *w, integer *perm,
doublereal *dlambda, doublereal *q2, integer *ldq2, doublereal *w, integer *perm,
integer *givptr, integer *givcol, doublereal *givnum, integer *indxp, integer *indx,
integer *info)
{
@ -35,7 +35,7 @@ int dlaed8_(integer *icompq, integer *k, integer *n, integer *qsiz, doublereal *
q -= q_offset;
--indxq;
--z__;
--dlamda;
--dlambda;
q2_dim1 = *ldq2;
q2_offset = 1 + q2_dim1;
q2 -= q2_offset;
@ -87,15 +87,15 @@ int dlaed8_(integer *icompq, integer *k, integer *n, integer *qsiz, doublereal *
}
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
dlamda[i__] = d__[indxq[i__]];
dlambda[i__] = d__[indxq[i__]];
w[i__] = z__[indxq[i__]];
}
i__ = 1;
j = *cutpnt + 1;
dlamrg_(&n1, &n2, &dlamda[1], &c__1, &c__1, &indx[1]);
dlamrg_(&n1, &n2, &dlambda[1], &c__1, &c__1, &indx[1]);
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
d__[i__] = dlamda[indx[i__]];
d__[i__] = dlambda[indx[i__]];
z__[i__] = w[indx[i__]];
}
imax = idamax_(n, &z__[1], &c__1);
@ -183,7 +183,7 @@ L80:
} else {
++(*k);
w[*k] = z__[jlam];
dlamda[*k] = d__[jlam];
dlambda[*k] = d__[jlam];
indxp[*k] = jlam;
jlam = j;
}
@ -192,21 +192,21 @@ L80:
L100:
++(*k);
w[*k] = z__[jlam];
dlamda[*k] = d__[jlam];
dlambda[*k] = d__[jlam];
indxp[*k] = jlam;
L110:
if (*icompq == 0) {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
jp = indxp[j];
dlamda[j] = d__[jp];
dlambda[j] = d__[jp];
perm[j] = indxq[indx[jp]];
}
} else {
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
jp = indxp[j];
dlamda[j] = d__[jp];
dlambda[j] = d__[jp];
perm[j] = indxq[indx[jp]];
dcopy_(qsiz, &q[perm[j] * q_dim1 + 1], &c__1, &q2[j * q2_dim1 + 1], &c__1);
}
@ -214,10 +214,10 @@ L110:
if (*k < *n) {
if (*icompq == 0) {
i__1 = *n - *k;
dcopy_(&i__1, &dlamda[*k + 1], &c__1, &d__[*k + 1], &c__1);
dcopy_(&i__1, &dlambda[*k + 1], &c__1, &d__[*k + 1], &c__1);
} else {
i__1 = *n - *k;
dcopy_(&i__1, &dlamda[*k + 1], &c__1, &d__[*k + 1], &c__1);
dcopy_(&i__1, &dlambda[*k + 1], &c__1, &d__[*k + 1], &c__1);
i__1 = *n - *k;
dlacpy_((char *)"A", qsiz, &i__1, &q2[(*k + 1) * q2_dim1 + 1], ldq2, &q[(*k + 1) * q_dim1 + 1],
ldq, (ftnlen)1);

19
lib/linalg/dlaed9.cpp
View File

@ -4,7 +4,7 @@ extern "C" {
#include "lmp_f2c.h"
static integer c__1 = 1;
int dlaed9_(integer *k, integer *kstart, integer *kstop, integer *n, doublereal *d__, doublereal *q,
integer *ldq, doublereal *rho, doublereal *dlamda, doublereal *w, doublereal *s,
integer *ldq, doublereal *rho, doublereal *dlambda, doublereal *w, doublereal *s,
integer *lds, integer *info)
{
integer q_dim1, q_offset, s_dim1, s_offset, i__1, i__2;
@ -15,14 +15,13 @@ int dlaed9_(integer *k, integer *kstart, integer *kstop, integer *n, doublereal
extern doublereal dnrm2_(integer *, doublereal *, integer *);
extern int dcopy_(integer *, doublereal *, integer *, doublereal *, integer *),
dlaed4_(integer *, integer *, doublereal *, doublereal *, doublereal *, doublereal *,
doublereal *, integer *);
extern doublereal dlamc3_(doublereal *, doublereal *);
extern int xerbla_(char *, integer *, ftnlen);
doublereal *, integer *),
xerbla_(char *, integer *, ftnlen);
--d__;
q_dim1 = *ldq;
q_offset = 1 + q_dim1;
q -= q_offset;
--dlamda;
--dlambda;
--w;
s_dim1 = *lds;
s_offset = 1 + s_dim1;
@ -49,13 +48,9 @@ int dlaed9_(integer *k, integer *kstart, integer *kstop, integer *n, doublereal
if (*k == 0) {
return 0;
}
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
dlamda[i__] = dlamc3_(&dlamda[i__], &dlamda[i__]) - dlamda[i__];
}
i__1 = *kstop;
for (j = *kstart; j <= i__1; ++j) {
dlaed4_(k, &j, &dlamda[1], &w[1], &q[j * q_dim1 + 1], rho, &d__[j], info);
dlaed4_(k, &j, &dlambda[1], &w[1], &q[j * q_dim1 + 1], rho, &d__[j], info);
if (*info != 0) {
goto L120;
}
@ -77,11 +72,11 @@ int dlaed9_(integer *k, integer *kstart, integer *kstop, integer *n, doublereal
for (j = 1; j <= i__1; ++j) {
i__2 = j - 1;
for (i__ = 1; i__ <= i__2; ++i__) {
w[i__] *= q[i__ + j * q_dim1] / (dlamda[i__] - dlamda[j]);
w[i__] *= q[i__ + j * q_dim1] / (dlambda[i__] - dlambda[j]);
}
i__2 = *k;
for (i__ = j + 1; i__ <= i__2; ++i__) {
w[i__] *= q[i__ + j * q_dim1] / (dlamda[i__] - dlamda[j]);
w[i__] *= q[i__ + j * q_dim1] / (dlambda[i__] - dlambda[j]);
}
}
i__1 = *k;

10
lib/linalg/dlasd8.cpp
View File

@ -4,7 +4,7 @@ extern "C" {
#include "lmp_f2c.h"
static integer c__1 = 1;
static integer c__0 = 0;
static doublereal c_b8 = 1.;
static doublereal c_b7 = 1.;
int dlasd8_(integer *icompq, integer *k, doublereal *d__, doublereal *z__, doublereal *vf,
doublereal *vl, doublereal *difl, doublereal *difr, integer *lddifr, doublereal *dsigma,
doublereal *work, integer *info)
@ -62,19 +62,15 @@ int dlasd8_(integer *icompq, integer *k, doublereal *d__, doublereal *z__, doubl
}
return 0;
}
i__1 = *k;
for (i__ = 1; i__ <= i__1; ++i__) {
dsigma[i__] = dlamc3_(&dsigma[i__], &dsigma[i__]) - dsigma[i__];
}
iwk1 = 1;
iwk2 = iwk1 + *k;
iwk3 = iwk2 + *k;
iwk2i = iwk2 - 1;
iwk3i = iwk3 - 1;
rho = dnrm2_(k, &z__[1], &c__1);
dlascl_((char *)"G", &c__0, &c__0, &rho, &c_b8, k, &c__1, &z__[1], k, info, (ftnlen)1);
dlascl_((char *)"G", &c__0, &c__0, &rho, &c_b7, k, &c__1, &z__[1], k, info, (ftnlen)1);
rho *= rho;
dlaset_((char *)"A", k, &c__1, &c_b8, &c_b8, &work[iwk3], k, (ftnlen)1);
dlaset_((char *)"A", k, &c__1, &c_b7, &c_b7, &work[iwk3], k, (ftnlen)1);
i__1 = *k;
for (j = 1; j <= i__1; ++j) {
dlasd4_(k, &j, &dsigma[1], &z__[1], &work[iwk1], &rho, &d__[j], &work[iwk2], info);

6
lib/linalg/dlatrs.cpp
View File

@ -16,7 +16,7 @@ int dlatrs_(char *uplo, char *trans, char *diag, char *normin, integer *n, doubl
extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *, integer *);
doublereal xbnd;
integer imax;
doublereal tmax, tjjs, xmax, grow, sumj;
doublereal tmax, tjjs, xmax, grow, sumj, work[1];
extern int dscal_(integer *, doublereal *, doublereal *, integer *);
extern logical lsame_(char *, char *, ftnlen, ftnlen);
doublereal tscal, uscal;
@ -100,7 +100,7 @@ int dlatrs_(char *uplo, char *trans, char *diag, char *normin, integer *n, doubl
i__1 = *n;
for (j = 2; j <= i__1; ++j) {
i__2 = j - 1;
d__1 = dlange_((char *)"M", &i__2, &c__1, &a[j * a_dim1 + 1], &c__1, &sumj, (ftnlen)1);
d__1 = dlange_((char *)"M", &i__2, &c__1, &a[j * a_dim1 + 1], &c__1, work, (ftnlen)1);
tmax = max(d__1, tmax);
}
} else {
@ -108,7 +108,7 @@ int dlatrs_(char *uplo, char *trans, char *diag, char *normin, integer *n, doubl
for (j = 1; j <= i__1; ++j) {
i__2 = *n - j;
d__1 =
dlange_((char *)"M", &i__2, &c__1, &a[j + 1 + j * a_dim1], &c__1, &sumj, (ftnlen)1);
dlange_((char *)"M", &i__2, &c__1, &a[j + 1 + j * a_dim1], &c__1, work, (ftnlen)1);
tmax = max(d__1, tmax);
}
}

4
lib/linalg/drscl.cpp
View File

@ -7,8 +7,7 @@ int drscl_(integer *n, doublereal *sa, doublereal *sx, integer *incx)
doublereal mul, cden;
logical done;
doublereal cnum, cden1, cnum1;
extern int dscal_(integer *, doublereal *, doublereal *, integer *),
dlabad_(doublereal *, doublereal *);
extern int dscal_(integer *, doublereal *, doublereal *, integer *);
extern doublereal dlamch_(char *, ftnlen);
doublereal bignum, smlnum;
--sx;
@ -17,7 +16,6 @@ int drscl_(integer *n, doublereal *sa, doublereal *sx, integer *incx)
}
smlnum = dlamch_((char *)"S", (ftnlen)1);
bignum = 1. / smlnum;
dlabad_(&smlnum, &bignum);
cden = *sa;
cnum = 1.;
L10:

26
lib/linalg/ilaenv.cpp
View File

@ -3,8 +3,8 @@ extern "C" {
#endif
#include "lmp_f2c.h"
static integer c__1 = 1;
static real c_b176 = (float)0.;
static real c_b177 = (float)1.;
static real c_b179 = (float)0.;
static real c_b180 = (float)1.;
static integer c__0 = 0;
integer ilaenv_(integer *ispec, char *name__, char *opts, integer *n1, integer *n2, integer *n3,
integer *n4, ftnlen name_len, ftnlen opts_len)
@ -201,6 +201,12 @@ L50:
} else {
nb = 64;
}
} else if (s_lmp_cmp(subnam + 3, (char *)"QP3RK", (ftnlen)4, (ftnlen)5) == 0) {
if (sname) {
nb = 32;
} else {
nb = 32;
}
}
} else if (s_lmp_cmp(c2, (char *)"PO", (ftnlen)2, (ftnlen)2) == 0) {
if (s_lmp_cmp(c3, (char *)"TRF", (ftnlen)3, (ftnlen)3) == 0) {
@ -402,6 +408,12 @@ L60:
} else {
nbmin = 2;
}
} else if (s_lmp_cmp(subnam + 3, (char *)"QP3RK", (ftnlen)4, (ftnlen)5) == 0) {
if (sname) {
nbmin = 2;
} else {
nbmin = 2;
}
}
} else if (s_lmp_cmp(c2, (char *)"SY", (ftnlen)2, (ftnlen)2) == 0) {
if (s_lmp_cmp(c3, (char *)"TRF", (ftnlen)3, (ftnlen)3) == 0) {
@ -493,6 +505,12 @@ L70:
} else {
nx = 128;
}
} else if (s_lmp_cmp(subnam + 3, (char *)"QP3RK", (ftnlen)4, (ftnlen)5) == 0) {
if (sname) {
nx = 128;
} else {
nx = 128;
}
}
} else if (s_lmp_cmp(c2, (char *)"SY", (ftnlen)2, (ftnlen)2) == 0) {
if (sname && s_lmp_cmp(c3, (char *)"TRD", (ftnlen)3, (ftnlen)3) == 0) {
@ -555,13 +573,13 @@ L130:
L140:
ret_val = 1;
if (ret_val == 1) {
ret_val = ieeeck_(&c__1, &c_b176, &c_b177);
ret_val = ieeeck_(&c__1, &c_b179, &c_b180);
}
return ret_val;
L150:
ret_val = 1;
if (ret_val == 1) {
ret_val = ieeeck_(&c__0, &c_b176, &c_b177);
ret_val = ieeeck_(&c__0, &c_b179, &c_b180);
}
return ret_val;
L160:

18
lib/linalg/zlaed8.cpp
View File

@ -5,7 +5,7 @@ extern "C" {
static doublereal c_b3 = -1.;
static integer c__1 = 1;
int zlaed8_(integer *k, integer *n, integer *qsiz, doublecomplex *q, integer *ldq, doublereal *d__,
doublereal *rho, integer *cutpnt, doublereal *z__, doublereal *dlamda,
doublereal *rho, integer *cutpnt, doublereal *z__, doublereal *dlambda,
doublecomplex *q2, integer *ldq2, doublereal *w, integer *indxp, integer *indx,
integer *indxq, integer *perm, integer *givptr, integer *givcol, doublereal *givnum,
integer *info)
@ -35,7 +35,7 @@ int zlaed8_(integer *k, integer *n, integer *qsiz, doublecomplex *q, integer *ld
q -= q_offset;
--d__;
--z__;
--dlamda;
--dlambda;
q2_dim1 = *ldq2;
q2_offset = 1 + q2_dim1;
q2 -= q2_offset;
@ -86,15 +86,15 @@ int zlaed8_(integer *k, integer *n, integer *qsiz, doublecomplex *q, integer *ld
}
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
dlamda[i__] = d__[indxq[i__]];
dlambda[i__] = d__[indxq[i__]];
w[i__] = z__[indxq[i__]];
}
i__ = 1;
j = *cutpnt + 1;
dlamrg_(&n1, &n2, &dlamda[1], &c__1, &c__1, &indx[1]);
dlamrg_(&n1, &n2, &dlambda[1], &c__1, &c__1, &indx[1]);
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
d__[i__] = dlamda[indx[i__]];
d__[i__] = dlambda[indx[i__]];
z__[i__] = w[indx[i__]];
}
imax = idamax_(n, &z__[1], &c__1);
@ -173,7 +173,7 @@ L70:
} else {
++(*k);
w[*k] = z__[jlam];
dlamda[*k] = d__[jlam];
dlambda[*k] = d__[jlam];
indxp[*k] = jlam;
jlam = j;
}
@ -182,19 +182,19 @@ L70:
L90:
++(*k);
w[*k] = z__[jlam];
dlamda[*k] = d__[jlam];
dlambda[*k] = d__[jlam];
indxp[*k] = jlam;
L100:
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
jp = indxp[j];
dlamda[j] = d__[jp];
dlambda[j] = d__[jp];
perm[j] = indxq[indx[jp]];
zcopy_(qsiz, &q[perm[j] * q_dim1 + 1], &c__1, &q2[j * q2_dim1 + 1], &c__1);
}
if (*k < *n) {
i__1 = *n - *k;
dcopy_(&i__1, &dlamda[*k + 1], &c__1, &d__[*k + 1], &c__1);
dcopy_(&i__1, &dlambda[*k + 1], &c__1, &d__[*k + 1], &c__1);
i__1 = *n - *k;
zlacpy_((char *)"A", qsiz, &i__1, &q2[(*k + 1) * q2_dim1 + 1], ldq2, &q[(*k + 1) * q_dim1 + 1], ldq,
(ftnlen)1);