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

update linalg to version 1.4 corresponding to LAPACK 3.12.1

Browse Source
This commit is contained in:
Axel Kohlmeyer
2025-01-08 15:52:30 -05:00
parent 0abb371fbe
commit 8ec9f37611
100 changed files with 1424 additions and 2031 deletions
Download Patch File Download Diff File Expand all files Collapse all files

48
lib/linalg/dlanv2.cpp
View File

@ -2,16 +2,28 @@
extern "C" {
#endif
#include "lmp_f2c.h"
static doublereal c_b3 = 1.;
static doublereal c_b6 = 1.;
int dlanv2_(doublereal *a, doublereal *b, doublereal *c__, doublereal *d__, doublereal *rt1r,
doublereal *rt1i, doublereal *rt2r, doublereal *rt2i, doublereal *cs, doublereal *sn)
{
integer i__1;
doublereal d__1, d__2;
double d_lmp_sign(doublereal *, doublereal *), sqrt(doublereal);
double log(doublereal), pow_lmp_di(doublereal *, integer *), d_lmp_sign(doublereal *, doublereal *),
sqrt(doublereal);
doublereal p, z__, aa, bb, cc, dd, cs1, sn1, sab, sac, eps, tau, temp, scale, bcmax, bcmis,
sigma;
extern doublereal dlapy2_(doublereal *, doublereal *), dlamch_(char *, ftnlen);
integer count;
doublereal safmn2;
extern doublereal dlapy2_(doublereal *, doublereal *);
doublereal safmx2;
extern doublereal dlamch_(char *, ftnlen);
doublereal safmin;
safmin = dlamch_((char *)"S", (ftnlen)1);
eps = dlamch_((char *)"P", (ftnlen)1);
d__1 = dlamch_((char *)"B", (ftnlen)1);
i__1 = (integer)(log(safmin / eps) / log(dlamch_((char *)"B", (ftnlen)1)) / 2.);
safmn2 = pow_lmp_di(&d__1, &i__1);
safmx2 = 1. / safmn2;
if (*c__ == 0.) {
*cs = 1.;
*sn = 0.;
@ -23,7 +35,7 @@ int dlanv2_(doublereal *a, doublereal *b, doublereal *c__, doublereal *d__, doub
*a = temp;
*b = -(*c__);
*c__ = 0.;
} else if (*a - *d__ == 0. && d_lmp_sign(&c_b3, b) != d_lmp_sign(&c_b3, c__)) {
} else if (*a - *d__ == 0. && d_lmp_sign(&c_b6, b) != d_lmp_sign(&c_b6, c__)) {
*cs = 1.;
*sn = 0.;
} else {
@ -32,7 +44,7 @@ int dlanv2_(doublereal *a, doublereal *b, doublereal *c__, doublereal *d__, doub
d__1 = abs(*b), d__2 = abs(*c__);
bcmax = max(d__1, d__2);
d__1 = abs(*b), d__2 = abs(*c__);
bcmis = min(d__1, d__2) * d_lmp_sign(&c_b3, b) * d_lmp_sign(&c_b3, c__);
bcmis = min(d__1, d__2) * d_lmp_sign(&c_b6, b) * d_lmp_sign(&c_b6, c__);
d__1 = abs(p);
scale = max(d__1, bcmax);
z__ = p / scale * p + bcmax / scale * bcmis;
@ -47,24 +59,44 @@ int dlanv2_(doublereal *a, doublereal *b, doublereal *c__, doublereal *d__, doub
*b -= *c__;
*c__ = 0.;
} else {
count = 0;
sigma = *b + *c__;
L10:
++count;
d__1 = abs(temp), d__2 = abs(sigma);
scale = max(d__1, d__2);
if (scale >= safmx2) {
sigma *= safmn2;
temp *= safmn2;
if (count <= 20) {
goto L10;
}
}
if (scale <= safmn2) {
sigma *= safmx2;
temp *= safmx2;
if (count <= 20) {
goto L10;
}
}
p = temp * .5;
tau = dlapy2_(&sigma, &temp);
*cs = sqrt((abs(sigma) / tau + 1.) * .5);
*sn = -(p / (tau * *cs)) * d_lmp_sign(&c_b3, &sigma);
*sn = -(p / (tau * *cs)) * d_lmp_sign(&c_b6, &sigma);
aa = *a * *cs + *b * *sn;
bb = -(*a) * *sn + *b * *cs;
cc = *c__ * *cs + *d__ * *sn;
dd = -(*c__) * *sn + *d__ * *cs;
*a = aa * *cs + cc * *sn;
*b = bb * *cs + dd * *sn;
*c__ = -aa * *sn + cc * *cs;
*c__ = -(aa * *sn) + cc * *cs;
*d__ = -bb * *sn + dd * *cs;
temp = (*a + *d__) * .5;
*a = temp;
*d__ = temp;
if (*c__ != 0.) {
if (*b != 0.) {
if (d_lmp_sign(&c_b3, b) == d_lmp_sign(&c_b3, c__)) {
if (d_lmp_sign(&c_b6, b) == d_lmp_sign(&c_b6, c__)) {
sab = sqrt((abs(*b)));
sac = sqrt((abs(*c__)));
d__1 = sab * sac;