Skip to content

Commit ddd4db1

Browse files
hjmjohnsonhjmjohnson
committed
COMP: Make FEM ITPACK tqlrat self-contained; drop dead determ_
itk::fem's ITPACK solver (dsrc2c.c eigvns_) called v3p_netlib's EISPACK tqlrat, which v3p_netlib removes under ITK_FUTURE_LEGACY_REMOVE, breaking the FEM link. Vendor self-contained, renamed copies of tqlrat/pythag/epslon (itkfem_tqlrat/pythag/epslon.c, symbols itkfem_*_) so FEM depends on no v3p_netlib eispack symbol in any configuration. The copies carry their own f2c typedefs and d_sign; integer is long, matching what dsrc2c.c passes. Also remove determ_, an unused ITPACK routine orphaned when its only caller eqrt1s was previously removed for copyright.
1 parent defccad commit ddd4db1

7 files changed

Lines changed: 572 additions & 46 deletions

File tree

‎Modules/Numerics/FEM/include/itpack.h‎

Lines changed: 0 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -525,16 +525,6 @@ extern "C"
525525
extern doublereal
526526
ddot_(integer * n, doublereal * dx, integer * incx, doublereal * dy, integer * incy);
527527

528-
/**
529-
* Subroutine that computes the determinant of a symmetric tridiagonal matrix
530-
* given by tri. det(tri - xlmda*i) = 0
531-
* \param n order of tridiagonal system
532-
* \param tri symmetric tridiagonal matrix of order n
533-
* \param xlmda argument for characteristic equation
534-
*/
535-
extern doublereal
536-
determ_(integer * n, doublereal * tri, doublereal * xlmda);
537-
538528
/**
539529
* Obtain default parameters
540530
* \param iparm array of 12 integer parameters

‎Modules/Numerics/FEM/src/CMakeLists.txt‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,9 @@
11
set(
22
ITKFEM_SRCS
33
dsrc2c.c
4+
itkfem_tqlrat.c
5+
itkfem_pythag.c
6+
itkfem_epslon.c
47
itkFEMElement2DC0LinearLine.cxx
58
itkFEMElement2DC0LinearLineStress.cxx
69
itkFEMElement2DC0LinearQuadrilateral.cxx

‎Modules/Numerics/FEM/src/dsrc2c.c‎

Lines changed: 2 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -26,6 +26,7 @@
2626

2727
#define V3P_NETLIB_SRC
2828
#include "v3p_netlib.h"
29+
#include "itkfem_eispack.h"
2930

3031
/* Modified by Peter Vanroose, Oct 2003: manual optimisation and clean-up */
3132

@@ -38,8 +39,6 @@ time(long * timer); /* #include <time.h> */
3839
extern doublereal
3940
cheby_(doublereal *, doublereal *, doublereal *, integer *, doublereal *, doublereal *);
4041
extern doublereal
41-
determ_(integer *, doublereal *, doublereal *);
42-
extern doublereal
4342
eigvns_(integer *, doublereal *, doublereal *, doublereal *, integer *);
4443
extern doublereal
4544
itpackddot_(integer *, doublereal *, integer *, doublereal *, integer *);
@@ -4433,39 +4432,6 @@ itpackddot_(integer * n, doublereal * dx, integer * incx, doublereal * dy, integ
44334432
return ret_val;
44344433
} /* itpackddot_ */
44354434

4436-
doublereal
4437-
determ_(integer * n, doublereal * tri, doublereal * xlmda)
4438-
{
4439-
/* Local variables */
4440-
static integer l;
4441-
static doublereal d1, d2, d3;
4442-
static integer icnt;
4443-
4444-
/* THIS SUBROUTINE COMPUTES THE DETERMINANT OF A SYMMETRIC */
4445-
/* TRIDIAGONAL MATRIX GIVEN BY TRI. DET(TRI - XLMDA*I) = 0 */
4446-
4447-
/* ... PARAMETER LIST */
4448-
4449-
/* N ORDER OF TRIDIAGONAL SYSTEM */
4450-
/* TRI SYMMETRIC TRIDIAGONAL MATRIX OF ORDER N */
4451-
/* XLMDA ARGUMENT FOR CHARACTERISTIC EQUATION */
4452-
4453-
d2 = tri[(*n << 1) - 2] - *xlmda;
4454-
d1 = d2 * (tri[(*n << 1) - 4] - *xlmda) - tri[(*n << 1) - 1];
4455-
if (*n == 2)
4456-
return d1;
4457-
4458-
for (icnt = 2; icnt < *n; ++icnt)
4459-
{
4460-
l = *n - icnt + 1;
4461-
d3 = d2;
4462-
d2 = d1;
4463-
d1 = (tri[((l - 1) << 1) - 2] - *xlmda) * d2 - d3 * tri[(l << 1) - 1];
4464-
}
4465-
4466-
return d1;
4467-
} /* determ_ */
4468-
44694435
/* Subroutine */
44704436
int
44714437
dfault_(integer * iparm, doublereal * rparm)
@@ -4782,7 +4748,7 @@ eigvns_(integer * n, doublereal * tri, doublereal * d, doublereal * e2, integer
47824748
}
47834749

47844750
// eqrt1s_(d, e2, n, &c__1, &c__0, ier);
4785-
v3p_netlib_tqlrat_(n, d, e2, ier);
4751+
itkfem_tqlrat_(n, d, e2, ier);
47864752

47874753
return -d[0];
47884754
} /* eigvns_ */

‎Modules/Numerics/FEM/src/itkfem_eispack.h‎

Lines changed: 42 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,42 @@
1+
/*=========================================================================
2+
*
3+
* Copyright NumFOCUS
4+
*
5+
* Licensed under the Apache License, Version 2.0 (the "License");
6+
* you may not use this file except in compliance with the License.
7+
* You may obtain a copy of the License at
8+
*
9+
* https://www.apache.org/licenses/LICENSE-2.0.txt
10+
*
11+
* Unless required by applicable law or agreed to in writing, software
12+
* distributed under the License is distributed on an "AS IS" BASIS,
13+
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14+
* See the License for the specific language governing permissions and
15+
* limitations under the License.
16+
*
17+
*=========================================================================*/
18+
#ifndef itkfem_eispack_h
19+
#define itkfem_eispack_h
20+
21+
/* FEM-private EISPACK eigenvalue helpers (renamed tqlrat/pythag/epslon).
22+
* v3p_netlib drops eispack under ITK_FUTURE_LEGACY_REMOVE; the itk::fem ITPACK
23+
* solver (dsrc2c.c eigvns_) needs tqlrat, so FEM owns self-contained copies and
24+
* no longer depends on the v3p_netlib eispack in any configuration. */
25+
26+
#ifdef __cplusplus
27+
extern "C"
28+
{
29+
#endif
30+
31+
int
32+
itkfem_tqlrat_(long * n, double * d, double * e2, long * ierr);
33+
double
34+
itkfem_pythag_(double * a, double * b);
35+
double
36+
itkfem_epslon_(double * x);
37+
38+
#ifdef __cplusplus
39+
}
40+
#endif
41+
42+
#endif

‎Modules/Numerics/FEM/src/itkfem_epslon.c‎

Lines changed: 94 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,94 @@
1+
/*=========================================================================
2+
*
3+
* Copyright NumFOCUS
4+
*
5+
* Licensed under the Apache License, Version 2.0 (the "License");
6+
* you may not use this file except in compliance with the License.
7+
* You may obtain a copy of the License at
8+
*
9+
* https://www.apache.org/licenses/LICENSE-2.0.txt
10+
*
11+
* Unless required by applicable law or agreed to in writing, software
12+
* distributed under the License is distributed on an "AS IS" BASIS,
13+
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14+
* See the License for the specific language governing permissions and
15+
* limitations under the License.
16+
*
17+
*=========================================================================*/
18+
19+
/* FEM-private copy of the EISPACK epslon routine, renamed to itkfem_epslon_.
20+
* v3p_netlib drops eispack under ITK_FUTURE_LEGACY_REMOVE, but the itk::fem
21+
* ITPACK solver (dsrc2c.c) still needs it, so FEM owns a self-contained copy
22+
* in every configuration. */
23+
24+
#include "itkfem_eispack.h"
25+
#ifdef __cplusplus
26+
extern "C"
27+
{
28+
#endif
29+
typedef double doublereal;
30+
#ifndef abs
31+
# define abs(x) ((x) >= 0 ? (x) : -(x))
32+
#endif
33+
34+
/*< double precision function epslon (x) >*/
35+
doublereal
36+
itkfem_epslon_(doublereal * x)
37+
{
38+
/* System generated locals */
39+
doublereal ret_val, d__1;
40+
41+
/* Local variables */
42+
doublereal a, b, c__, eps;
43+
44+
/*< double precision x >*/
45+
46+
/* estimate unit roundoff in quantities of size x. */
47+
48+
/*< double precision a,b,c,eps >*/
49+
50+
/* this program should function properly on all systems */
51+
/* satisfying the following two assumptions, */
52+
/* 1. the base used in representing floating point */
53+
/* numbers is not a power of three. */
54+
/* 2. the quantity a in statement 10 is represented to */
55+
/* the accuracy used in floating point variables */
56+
/* that are stored in memory. */
57+
/* the statement number 10 and the go to 10 are intended to */
58+
/* force optimizing compilers to generate code satisfying */
59+
/* assumption 2. */
60+
/* under these assumptions, it should be true that, */
61+
/* a is not exactly equal to four-thirds, */
62+
/* b has a zero for its last bit or digit, */
63+
/* c is not exactly equal to one, */
64+
/* eps measures the separation of 1.0 from */
65+
/* the next larger floating point number. */
66+
/* the developers of eispack would appreciate being informed */
67+
/* about any systems where these assumptions do not hold. */
68+
69+
/* this version dated 4/6/83. */
70+
71+
/*< a = 4.0d0/3.0d0 >*/
72+
a = 1.3333333333333333;
73+
/*< 10 b = a - 1.0d0 >*/
74+
L10:
75+
b = a - 1.;
76+
/*< c = b + b + b >*/
77+
c__ = b + b + b;
78+
/*< eps = dabs(c-1.0d0) >*/
79+
eps = (d__1 = c__ - 1., abs(d__1));
80+
/*< if (eps .eq. 0.0d0) go to 10 >*/
81+
if (eps == 0.)
82+
{
83+
goto L10;
84+
}
85+
/*< epslon = eps*dabs(x) >*/
86+
ret_val = eps * abs(*x);
87+
/*< return >*/
88+
return ret_val;
89+
/*< end >*/
90+
} /* itkfem_epslon_ */
91+
92+
#ifdef __cplusplus
93+
}
94+
#endif

‎Modules/Numerics/FEM/src/itkfem_pythag.c‎

Lines changed: 101 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,101 @@
1+
/*=========================================================================
2+
*
3+
* Copyright NumFOCUS
4+
*
5+
* Licensed under the Apache License, Version 2.0 (the "License");
6+
* you may not use this file except in compliance with the License.
7+
* You may obtain a copy of the License at
8+
*
9+
* https://www.apache.org/licenses/LICENSE-2.0.txt
10+
*
11+
* Unless required by applicable law or agreed to in writing, software
12+
* distributed under the License is distributed on an "AS IS" BASIS,
13+
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14+
* See the License for the specific language governing permissions and
15+
* limitations under the License.
16+
*
17+
*=========================================================================*/
18+
19+
/* FEM-private copy of the EISPACK pythag routine, renamed to itkfem_pythag_.
20+
* v3p_netlib drops eispack under ITK_FUTURE_LEGACY_REMOVE, but the itk::fem
21+
* ITPACK solver (dsrc2c.c) still needs it, so FEM owns a self-contained copy
22+
* in every configuration. */
23+
24+
#include "itkfem_eispack.h"
25+
#ifdef __cplusplus
26+
extern "C"
27+
{
28+
#endif
29+
typedef double doublereal;
30+
#ifndef abs
31+
# define abs(x) ((x) >= 0 ? (x) : -(x))
32+
#endif
33+
#ifndef max
34+
# define max(a, b) ((a) >= (b) ? (a) : (b))
35+
#endif
36+
#ifndef min
37+
# define min(a, b) ((a) <= (b) ? (a) : (b))
38+
#endif
39+
40+
/*< double precision function pythag(a,b) >*/
41+
doublereal
42+
itkfem_pythag_(doublereal * a, doublereal * b)
43+
{
44+
/* System generated locals */
45+
doublereal ret_val, d__1, d__2, d__3;
46+
47+
/* Local variables */
48+
doublereal p, r__, s, t, u;
49+
50+
/*< double precision a,b >*/
51+
52+
/* finds dsqrt(a**2+b**2) without overflow or destructive underflow */
53+
54+
/*< double precision p,r,s,t,u >*/
55+
/*< p = dmax1(dabs(a),dabs(b)) >*/
56+
/* Computing MAX */
57+
d__1 = abs(*a), d__2 = abs(*b);
58+
p = max(d__1, d__2);
59+
/*< if (p .eq. 0.0d0) go to 20 >*/
60+
if (p == 0.)
61+
{
62+
goto L20;
63+
}
64+
/*< r = (dmin1(dabs(a),dabs(b))/p)**2 >*/
65+
/* Computing MIN */
66+
d__2 = abs(*a), d__3 = abs(*b);
67+
/* Computing 2nd power */
68+
d__1 = min(d__2, d__3) / p;
69+
r__ = d__1 * d__1;
70+
/*< 10 continue >*/
71+
L10:
72+
/*< t = 4.0d0 + r >*/
73+
t = r__ + 4.;
74+
/*< if (t .eq. 4.0d0) go to 20 >*/
75+
if (t == 4.)
76+
{
77+
goto L20;
78+
}
79+
/*< s = r/t >*/
80+
s = r__ / t;
81+
/*< u = 1.0d0 + 2.0d0*s >*/
82+
u = s * 2. + 1.;
83+
/*< p = u*p >*/
84+
p = u * p;
85+
/*< r = (s/u)**2 * r >*/
86+
/* Computing 2nd power */
87+
d__1 = s / u;
88+
r__ = d__1 * d__1 * r__;
89+
/*< go to 10 >*/
90+
goto L10;
91+
/*< 20 pythag = p >*/
92+
L20:
93+
ret_val = p;
94+
/*< return >*/
95+
return ret_val;
96+
/*< end >*/
97+
} /* itkfem_pythag_ */
98+
99+
#ifdef __cplusplus
100+
}
101+
#endif

0 commit comments

Comments
 (0)