← Back to branch summary

~ubuntu-branches/debian/sid/arpack/sid

« back to all changes in this revision

Viewing changes to detect_arpack_bug.m4

  • Committer: Package Import Robot
  • Author(s): Sylvestre Ledru
  • Date: 2012-02-22 11:05:03 UTC
  • mfrom: (1.2.5)
  • Revision ID: package-import@ubuntu.com-20120222110503-h6ux3f5ilm5q76w0
Tags: 3.1.0-1
New upstream release

Show diffs side-by-side

added added

removed removed

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