1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
|
dnl
dnl Check whether ARPACK works (does not crash)
dnl
dnl Using a pure Fortran program doesn't seem to crash when linked
dnl with the buggy ARPACK library but the C++ program does. Maybe
dnl it is the memory allocation that exposes the bug and using statically
dnl allocated arrays in Fortran does not?
dnl
dnl This code is not used by arpack-ng itself.
dnl This is a macro for applications using arpack to detect that the version
dnl of arpack behave correctly (ie not arpack-ng)
dnl This is the work of Rik <rik@octave.org>
dnl
dnl This code is released under the same license as arpack
dnl
AC_DEFUN([CHECK_ARPACK_OK], [
AC_LANG_PUSH(C++)
AC_CACHE_CHECK([whether the arpack library works],
[cv_lib_arpack_ok], [
AC_RUN_IFELSE([AC_LANG_PROGRAM([[
// External functions from ARPACK library
extern "C" int
F77_FUNC (dnaupd, DNAUPD) (int&, const char *, const int&, const char *,
int&, const double&, double*, const int&,
double*, const int&, int*, int*, double*,
double*, const int&, int&, long int, long int);
extern "C" int
F77_FUNC (dneupd, DNEUPD) (const int&, const char *, int*, double*,
double*, double*, const int&,
const double&, const double&, double*,
const char*, const int&, const char *,
int&, const double&, double*, const int&,
double*, const int&, int*, int*, double*,
double*, const int&, int&, long int,
long int, long int);
extern "C" int
F77_FUNC (dgemv, DGEMV) (const char *, const int&, const int&,
const double&, const double*, const int&,
const double*, const int&, const double&,
double*, const int&, long int);
#include <cfloat>
void
doit (void)
{
// Based on the octave function EigsRealNonSymmetricMatrix from
// liboctave/eigs-base.cc.
// Problem matrix. See bug #31479
int n = 4;
double *m = new double [n * n];
m[0] = 1, m[4] = 0, m[8] = 0, m[12] = -1;
m[1] = 0, m[5] = 1, m[9] = 0, m[13] = 0;
m[2] = 0, m[6] = 0, m[10] = 1, m[14] = 0;
m[3] = 0, m[7] = 0, m[11] = 2, m[15] = 1;
double *resid = new double [4];
resid[0] = 0.960966;
resid[1] = 0.741195;
resid[2] = 0.150143;
resid[3] = 0.868067;
int *ip = new int [11];
ip[0] = 1; // ishift
ip[1] = 0; // ip[1] not referenced
ip[2] = 300; // mxiter, maximum number of iterations
ip[3] = 1; // NB blocksize in recurrence
ip[4] = 0; // nconv, number of Ritz values that satisfy convergence
ip[5] = 0; // ip[5] not referenced
ip[6] = 1; // mode
ip[7] = 0; // ip[7] to ip[10] are return values
ip[8] = 0;
ip[9] = 0;
ip[10] = 0;
int *ipntr = new int [14];
int k = 1;
int p = 3;
int lwork = 3 * p * (p + 2);
double *v = new double [n * (p + 1)];
double *workl = new double [lwork + 1];
double *workd = new double [3 * n + 1];
int ido = 0;
int info = 0;
double tol = DBL_EPSILON;
do
{
F77_FUNC (dnaupd, DNAUPD) (ido, "I", n, "LM", k, tol, resid, p,
v, n, ip, ipntr, workd, workl, lwork,
info, 1L, 2L);
if (ido == -1 || ido == 1 || ido == 2)
{
double *x = workd + ipntr[0] - 1;
double *y = workd + ipntr[1] - 1;
F77_FUNC (dgemv, DGEMV) ("N", n, n, 1.0, m, n, x, 1, 0.0,
y, 1, 1L);
}
else
{
if (info < 0)
{
return; // Error
}
break;
}
}
while (1);
int *sel = new int [p];
// In Octave, the dimensions of dr and di are k+1, but k+2 avoids segfault
double *dr = new double [k + 1];
double *di = new double [k + 1];
double *workev = new double [3 * p];
for (int i = 0; i < k + 1; i++)
dr[i] = di[i] = 0.;
int rvec = 1;
double sigmar = 0.0;
double sigmai = 0.0;
// In Octave, this is n*(k+1), but k+2 avoids segfault
double *z = new double [n * (k + 1)];
F77_FUNC (dneupd, DNEUPD) (rvec, "A", sel, dr, di, z, n, sigmar,
sigmai, workev, "I", n, "LM", k, tol,
resid, p, v, n, ip, ipntr, workd,
workl, lwork, info, 1L, 1L, 2L);
}
]], [[
for (int i = 0; i < 10; i++)
doit ();
]])],
[cv_lib_arpack_ok=yes],
[cv_lib_arpack_ok=no],
[cv_lib_arpack_ok=yes])])
AC_LANG_POP(C++)
if test "$cv_lib_arpack_ok" = "yes"; then
$1
else
$2
fi
])
|
