Mercurial > octave-nkf

comparison m4/acinclude.m4 @ 14147:71e28fda7be9 stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
use C++ program to test ARPACK * acinclude.m4 (OCTAVE_CHECK_ARPACK_OK): Use C++ instead of Fortran for the test program.
author John W. Eaton <jwe@octave.org>
date Thu, 05 Jan 2012 00:10:19 -0500
parents 834df9f10963
children 99428221b4e1
comparison
equal deleted inserted replaced
14146:2f05228cb959 14147:71e28fda7be9
1015 fi 1015 fi
1016 ]) 1016 ])
1017 dnl 1017 dnl
1018 dnl Check whether ARPACK works (does not crash) 1018 dnl Check whether ARPACK works (does not crash)
1019 dnl 1019 dnl
1020 dnl Using a pure Fortran program doesn't seem to crash when linked
1021 dnl with the buggy ARPACK library but the C++ program does. Maybe
1022 dnl it is the memory allocation that exposes the bug and using statically
1023 dnl allocated arrays in Fortran does not?
1024 dnl
1020 AC_DEFUN([OCTAVE_CHECK_ARPACK_OK], [ 1025 AC_DEFUN([OCTAVE_CHECK_ARPACK_OK], [
1021 AC_LANG_PUSH(Fortran 77) 1026 AC_LANG_PUSH(C++)
1022 AC_CACHE_CHECK([whether the arpack library works], 1027 AC_CACHE_CHECK([whether the arpack library works],
1023 [octave_cv_lib_arpack_ok], [ 1028 [octave_cv_lib_arpack_ok], [
1024 AC_RUN_IFELSE([AC_LANG_PROGRAM([], [[ 1029 AC_RUN_IFELSE([AC_LANG_PROGRAM([[
1025 * Based on function EigsRealNonSymmetricMatrix from liboctave/eigs-base.cc. 1030 // External functions from ARPACK library
1026 * Problem matrix. See bug #31479 1031 extern "C" int
1027 1032 F77_FUNC (dnaupd, DNAUPD) (int&, const char *, const int&, const char *,
1028 integer rvec, ido, info, ii, jj 1033 int&, const double&, double*, const int&,
1029 double precision tol, sigmar, sigmai 1034 double*, const int&, int*, int*, double*,
1030 1035 double*, const int&, int&, long int, long int);
1031 integer n, k, p, lwork 1036
1032 parameter (n = 4, k = 1, p = 3, lwork = 3*p*(p+2)) 1037 extern "C" int
1033 1038 F77_FUNC (dneupd, DNEUPD) (const int&, const char *, int*, double*,
1034 integer ip (11) 1039 double*, double*, const int&,
1035 integer ipntr (14) 1040 const double&, const double&, double*,
1036 integer sel (p) 1041 const char*, const int&, const char *,
1037 1042 int&, const double&, double*, const int&,
1038 double precision m (n, n) 1043 double*, const int&, int*, int*, double*,
1039 double precision resid (4) 1044 double*, const int&, int&, long int,
1040 double precision v (n*(p+1)); 1045 long int, long int);
1041 double precision workl (lwork+1) 1046
1042 double precision workd (3*n+1) 1047 extern "C" int
1043 * In Octave, the dimensions of dr and di are k+1, but k+2 avoids segfault 1048 F77_FUNC (dgemv, DGEMV) (const char *, const int&, const int&,
1044 double precision dr (k+1) 1049 const double&, const double*, const int&,
1045 double precision di (k+1) 1050 const double*, const int&, const double&,
1046 double precision workev (3*p) 1051 double*, const int&, long int);
1047 * In Octave, this is n*(k+1), but k+2 avoids segfault 1052
1048 double precision z (n*(k+1)) 1053 #include <cfloat>
1049 1054
1050 do 10 ii = 1, 100 1055 void
1051 1056 doit (void)
1052 m(1,1) = 1 1057 {
1053 m(2,1) = 0 1058 // Based on function EigsRealNonSymmetricMatrix from liboctave/eigs-base.cc.
1054 m(3,1) = 0 1059
1055 m(4,1) = 0 1060 // Problem matrix. See bug #31479
1056 1061 int n = 4;
1057 m(1,2) = 0 1062 double *m = new double [n * n];
1058 m(2,2) = 1 1063 m[0] = 1, m[4] = 0, m[8] = 0, m[12] = -1;
1059 m(3,2) = 0 1064 m[1] = 0, m[5] = 1, m[9] = 0, m[13] = 0;
1060 m(4,2) = 0 1065 m[2] = 0, m[6] = 0, m[10] = 1, m[14] = 0;
1061 1066 m[3] = 0, m[7] = 0, m[11] = 2, m[15] = 1;
1062 m(1,3) = 0 1067
1063 m(2,3) = 0 1068 double *resid = new double [4];
1064 m(3,3) = 1 1069
1065 m(4,3) = 2 1070 resid[0] = 0.960966;
1066 1071 resid[1] = 0.741195;
1067 m(1,4) = -1 1072 resid[2] = 0.150143;
1068 m(2,4) = 0 1073 resid[3] = 0.868067;
1069 m(3,4) = 0 1074
1070 m(4,4) = 1 1075 int *ip = new int [11];
1071 1076
1072 resid(1) = 0.960966 1077 ip[0] = 1; // ishift
1073 resid(1) = 0.741195 1078 ip[1] = 0; // ip[1] not referenced
1074 resid(1) = 0.150143 1079 ip[2] = 300; // mxiter, maximum number of iterations
1075 resid(1) = 0.868067 1080 ip[3] = 1; // NB blocksize in recurrence
1076 1081 ip[4] = 0; // nconv, number of Ritz values that satisfy convergence
1077 * ip(1) = ishift 1082 ip[5] = 0; // ip[5] not referenced
1078 * ip(2) is not referenced 1083 ip[6] = 1; // mode
1079 * ip(3) = maximum number of iterations 1084 ip[7] = 0; // ip[7] to ip[10] are return values
1080 * ip(4) = NB blocksize in recurrence 1085 ip[8] = 0;
1081 * ip(5) = nconv, number of Ritz values that satisfy convergence 1086 ip[9] = 0;
1082 * ip(6) is not referenced 1087 ip[10] = 0;
1083 * ip(7) = mode 1088
1084 * ip(8) to ip(11) are return values 1089 int *ipntr = new int [14];
1085 1090
1086 ip(1) = 1 1091 int k = 1;
1087 ip(2) = 0 1092 int p = 3;
1088 ip(3) = 300 1093 int lwork = 3 * p * (p + 2);
1089 ip(4) = 1 1094
1090 ip(5) = 0 1095 double *v = new double [n * (p + 1)];
1091 ip(6) = 0 1096 double *workl = new double [lwork + 1];
1092 ip(7) = 1 1097 double *workd = new double [3 * n + 1];
1093 ip(8) = 0 1098
1094 ip(9) = 0 1099 int ido = 0;
1095 ip(10) = 0 1100 int info = 0;
1096 ip(11) = 0 1101
1097 1102 double tol = DBL_EPSILON;
1098 ido = 0 1103
1099 info = 0 1104 do
1100 rvec = 1 1105 {
1101 sigmar = 0 1106 F77_FUNC (dnaupd, DNAUPD) (ido, "I", n, "LM", k, tol, resid, p,
1102 sigmai = 0 1107 v, n, ip, ipntr, workd, workl, lwork,
1103 tol = 2.0d-15 1108 info, 1L, 2L);
1104 1109
1105 1 continue 1110 if (ido == -1 || ido == 1 || ido == 2)
1106 1111 {
1107 call dnaupd (ido, 'I', n, 'LM', k, tol, resid, p, v, n, ip, 1112 double *x = workd + ipntr[0] - 1;
1108 $ ipntr, workd, workl, lwork, info) 1113 double *y = workd + ipntr[1] - 1;
1109 1114
1110 if (ido .eq. -1 .or. ido .eq. 1 .or. ido .eq. 2) then 1115 F77_FUNC (dgemv, DGEMV) ("N", n, n, 1.0, m, n, x, 1, 0.0,
1111 call dgemv ('N', n, n, 1.0, m, n, workd(ipntr(1)), 1116 y, 1, 1L);
1112 $ 1, 0.0, workd(ipntr(2)), 1) 1117 }
1113 goto 1 1118 else
1114 else if (info .lt. 0) then 1119 {
1115 * error 1120 if (info < 0)
1116 goto 9999 1121 {
1117 endif 1122 return; // Error
1118 1123 }
1119 do 10 jj = 1, k+1 1124
1120 dr(jj) = 0 1125 break;
1121 di(jj) = 0 1126 }
1122 10 continue 1127 }
1123 1128 while (1);
1124 call dneupd (rvec, "A", sel, dr, di, z, n, sigmar, sigmai, 1129
1125 $ workev,"I", n, "LM", k, tol, resid, p, v, n, ip, ipntr, 1130 int *sel = new int [p];
1126 $ workd,workl, lwork, info) 1131
1127 1132 // In Octave, the dimensions of dr and di are k+1, but k+2 avoids segfault
1128 100 continue 1133 double *dr = new double [k + 1];
1129 1134 double *di = new double [k + 1];
1130 9999 continue 1135 double *workev = new double [3 * p];
1136
1137 for (int i = 0; i < k + 1; i++)
1138 dr[i] = di[i] = 0.;
1139
1140 int rvec = 1;
1141
1142 double sigmar = 0.0;
1143 double sigmai = 0.0;
1144
1145 // In Octave, this is n*(k+1), but k+2 avoids segfault
1146 double *z = new double [n * (k + 1)];
1147
1148 F77_FUNC (dneupd, DNEUPD) (rvec, "A", sel, dr, di, z, n, sigmar,
1149 sigmai, workev, "I", n, "LM", k, tol,
1150 resid, p, v, n, ip, ipntr, workd,
1151 workl, lwork, info, 1L, 1L, 2L);
1152 }
1153 ]], [[
1154 for (int i = 0; i < 10; i++)
1155 doit ();
1131 ]])], 1156 ]])],
1132 [octave_cv_lib_arpack_ok=yes], 1157 [octave_cv_lib_arpack_ok=yes],
1133 [octave_cv_lib_arpack_ok=no], 1158 [octave_cv_lib_arpack_ok=no],
1134 [octave_cv_lib_arpack_ok=yes])]) 1159 [octave_cv_lib_arpack_ok=yes])])
1135 AC_LANG_POP(Fortran 77) 1160 AC_LANG_POP(C++)
1136 if test "$octave_cv_lib_arpack_ok" = "yes"; then 1161 if test "$octave_cv_lib_arpack_ok" = "yes"; then
1137 $1 1162 $1
1138 else 1163 else
1139 $2 1164 $2
1140 fi 1165 fi