update of the release date

* Version numbers need to be kept up-to-date

* arpack.pc:  includedir should match the corresponding ./configure option

* Override pkgincludedir instead of creating a new directory variable

* Supplying arguments to AC_OUTPUT is obsolete behavior

* Distribute the entire source regardless of build options

Otherwise you won't be able to build ICB/MPI/etc. if the tarball wasn't
originally built with those options, because those files will be
missing.  Missing files leads to more breakage.

DIST_SUBDIRS is handled automatically.

* Clean up conditional recursion into subdirectories

Don't need to check for ICB inside a subdirectory if it won't be
descended into from the parent directory.

* Fix various include path problems

Should be using $(top_srcdir), not $(top_builddir).
-I flags go in CPPFLAGS, not CFLAGS nor FCFLAGS.

* Minor rearrangement of autoconf macros

Put all AC_CONFIG_*s together near the end.
Put all AM_CONDITIONALs together near the end.  These should never be
conditionally expanded.

* Run Fortran flag tests against the correct compiler

* Simplify --enable arguments

AC_HELP_STRING has been replaced by AS_HELP_STRING.
Setting enable_foo to $enableval is an unnecessary self-assignment.

* Fix missing files in the distribution

* Simplify creation of arpackdef.h

The ax_gen_arpack.m4 macro would add a bunch of extra garbage to
arpackdef.h, such as PACKAGE_NAME, VERSION, HAVE_BLAS, etc. that would
clash if any other project tried to #include it.  Shorten it to only the
one #define that we care about.

This file is placed in $(top_builddir) so add that path to CPPFLAGS.
-I$(top_builddir) should be first.  Note that autoheader would do this,
but now we don't need to use autoheader.

Move #include lines inside the include guards.

* Need PKG_PROG_PKG_CONFIG if first PKG_CHECK_MODULES is conditional

* Populate the version and URL fields of arpack.pc from CMake

Still not installing this yet because it doesn't contain full dependency
info for LAPACK/BLAS.

* Bump version to 3.7.0

* Update shared library version

Updating the minor number (libtool "age") due to the new dbgini/staini
symbols.  Note that CMake doesn't follow libtool's current:revision:age
numbering scheme, so care should be taken with any future changes to
make sure they generate files with the same numbers.

* Enable checks for eigen3 headers

.github/ISSUE_TEMPLATE.md

@@ -0,0 +1,39 @@
+## Expected behavior
+
+Do this.
+
+## Actual behavior
+
+Do that.
+
+## Where/how to reproduce the problem
+
+- arpack-ng: release or commit (SHA)
+- OS: debian ? ubuntu ? mac ?...
+- compiler: gcc ? intel ? clang ? version ?...
+- environment: export FFLAGS ?...
+- configure: ./configure --enable-stuff ? cmake -DSTUFFS=ON ?
+- input data ?
+
+## Steps to reproduce the problem
+
+-
+-
+-
+
+## Error message
+
+Bad stuff.
+
+## Traces
+
+~/arpack-ng> tail -n 300 test-suite.log
+
+~/arpack-ng/build> tail -n 300 Testing/Temporary/LastTest.log
+
+## Callstack
+
+(gdb) bt
+
+## Notes, remarks
+

.github/PULL_REQUEST_TEMPLATE.md

@@ -0,0 +1,9 @@
+## Pull request purpose
+
+fixing issue #xx ? enhancement ? new feature ?...
+
+## Detailed changes proposed in this pull request
+
+-
+-
+-

.gitignore

@@ -16,8 +16,10 @@ config.status
 libtool
 .deps/
 arpack.pc
+arpackdef.h
 
 # Generated by `make`
+.dirstamp
 *.o
 *.mod
 *.lo
@@ -40,6 +42,8 @@ TESTS/bug_58_double
 TESTS/bug_79_double_complex
 TESTS/icb_arpack_c
 TESTS/icb_arpack_cpp
+TESTS/bug_142
+TESTS/bug_142_gen
 EXAMPLES/BAND/[sd][sn]bdr[123456]
 EXAMPLES/BAND/[cz]nbdr[1234]
 EXAMPLES/COMPLEX/[cz]ndrv[1234]
@@ -48,6 +52,9 @@ EXAMPLES/SIMPLE/[sd]ssimp
 EXAMPLES/SIMPLE/[sdcz]nsimp
 EXAMPLES/SVD/[sd]svd
 EXAMPLES/SYM/[sd]sdrv[123456]
+EXAMPLES/MATRIX_MARKET/arpackmm
+EXAMPLES/MATRIX_MARKET/resid.out
+EXAMPLES/MATRIX_MARKET/v.out
 PARPACK/EXAMPLES/MPI/p[sd]ndrv[13]
 PARPACK/EXAMPLES/MPI/p[sd]sdrv1
 PARPACK/EXAMPLES/MPI/p[cz]ndrv1

.travis.yml

@@ -1,4 +1,4 @@
-sudo: false
+sudo: true
 
 language: c
 
@@ -18,15 +18,119 @@ addons:
       - autoconf
       - libblas-dev
       - liblapack-dev
+      - diffutils
+      - findutils
+      - libeigen3-dev
+
+services:
+  - docker
+
+stages:
+  # order stages
+  - name: precise
+  - name: trusty
+  - name: xenial
+  - name: coverage
+  - name: interface64
 
 jobs:
   include:
-    - script: mkdir -p build && cd build && cmake -D EXAMPLES=ON -D MPI=OFF -D ICB=OFF               .. && make VERBOSE=1 && make test;
-    - script: mkdir -p build && cd build && cmake -D EXAMPLES=ON -D MPI=ON  -D ICB=OFF               .. && make VERBOSE=1 && make test;
-    - script: mkdir -p build && cd build && cmake -D EXAMPLES=ON -D MPI=OFF -D ICB=ON                .. && make VERBOSE=1 && make test;
-    - script: mkdir -p build && cd build && cmake -D EXAMPLES=ON -D MPI=ON  -D ICB=ON                .. && make VERBOSE=1 && make test;
-    - script: mkdir -p build && cd build && cmake -D EXAMPLES=ON -D MPI=ON  -D ICB=ON -DCOVERALLS=ON .. && make VERBOSE=1 && make test;
-    - script: ./bootstrap && ./configure                                                                && make VERBOSE=1 && make check && make distcheck;
-    - script: ./bootstrap && ./configure --enable-mpi                                                   && make VERBOSE=1 && make check && make distcheck;
-    - script: ./bootstrap && ./configure              --enable-icb                                      && make VERBOSE=1 && make check && make distcheck;
-    - script: ./bootstrap && ./configure --enable-mpi --enable-icb                                      && make VERBOSE=1 && make check && make distcheck;
+  # precise <=> test "older" systems, without ICB, without cmake (too old to be supported)
+  - stage: precise
+    dist: precise
+    script: ./bootstrap && ./configure                                                  && make VERBOSE=1 && make check && make distcheck;
+  - stage: precise
+    dist: precise
+    script: ./bootstrap && ./configure --enable-mpi                                     && make VERBOSE=1 && make check && make distcheck;
+  # trusty <=> test "older" systems, without ICB
+  - stage: trusty
+    dist: trusty
+    script: mkdir -p build && cd build && cmake -D EXAMPLES=ON -D MPI=OFF -D ICB=OFF .. && make VERBOSE=1 && make test  && make package_source;
+  - stage: trusty
+    dist: trusty
+    script: mkdir -p build && cd build && cmake -D EXAMPLES=ON -D MPI=ON  -D ICB=OFF .. && make VERBOSE=1 && make test  && make package_source;
+  - stage: trusty
+    dist: trusty
+    script: ./bootstrap && ./configure                                                  && make VERBOSE=1 && make check && make distcheck;
+  - stage: trusty
+    dist: trusty
+    script: ./bootstrap && ./configure --enable-mpi                                     && make VERBOSE=1 && make check && make distcheck;
+  # xenial <=> test "recent" systems, with and without ICB
+  - stage: xenial
+    dist: xenial
+    script: mkdir -p build && cd build && cmake -D EXAMPLES=ON -D MPI=OFF -D ICB=OFF .. && make VERBOSE=1 && make test  && make package_source;
+  - stage: xenial
+    dist: xenial
+    script: mkdir -p build && cd build && cmake -D EXAMPLES=ON -D MPI=ON  -D ICB=OFF .. && make VERBOSE=1 && make test  && make package_source;
+  - stage: xenial
+    dist: xenial
+    script: mkdir -p build && cd build && cmake -D EXAMPLES=ON -D MPI=OFF -D ICB=ON  .. && make VERBOSE=1 && make test  && make package_source;
+  - stage: xenial
+    dist: xenial
+    script: mkdir -p build && cd build && cmake -DEXAMPLES=ON -DMPI=OFF -DICBEXMM=ON .. && make VERBOSE=1 && make test  && make package_source;
+  - stage: xenial
+    dist: xenial
+    script: mkdir -p build && cd build && cmake -D EXAMPLES=ON -D MPI=ON  -D ICB=ON  .. && make VERBOSE=1 && make test  && make package_source;
+  - stage: xenial
+    dist: xenial
+    script: ./bootstrap && ./configure                                                  && make VERBOSE=1 && make check && make distcheck;
+  - stage: xenial
+    dist: xenial
+    script: ./bootstrap && ./configure --enable-mpi                                     && make VERBOSE=1 && make check && make distcheck;
+  - stage: xenial
+    dist: xenial
+    script: ./bootstrap && ./configure              --enable-icb                        && make VERBOSE=1 && make check && make distcheck;
+  - stage: xenial
+    dist: xenial
+    script: ./bootstrap && ./configure              --enable-icb-exmm                   && make VERBOSE=1 && make check && make distcheck;
+  - stage: xenial
+    dist: xenial
+    script: ./bootstrap && ./configure --enable-mpi --enable-icb                        && make VERBOSE=1 && make check && make distcheck;
+  # xenial <=> coverage: "recent" systems with ICB
+  - stage: coverage
+    dist: xenial
+    script: mkdir -p build && cd build && cmake -DEXAMPLES=ON -DMPI=ON -DICBEXMM=ON -DCOVERALLS=ON .. && make && make test coveralls > /dev/null;
+  # xenial <=> interface64: "recent" systems with ICB + MKL + ILP64 (need debian/testing to get MKL-ILP64)
+  - stage: interface64
+    allow_failure:
+    dist: xenial
+    script: |
+      sudo docker pull debian                                                                                           \
+      &&                                                                                                                \
+      sudo docker run debian /bin/bash -c                                                                               \
+      "cat /etc/os-release                                                                                           && \
+       more /etc/apt/sources.list                                                                                    && \
+       sed -e 's/stretch/testing/'            -i /etc/apt/sources.list                                               && \
+       sed -e 's/main/main non-free contrib/' -i /etc/apt/sources.list                                               && \
+       more /etc/apt/sources.list                                                                                    && \
+       export DEBIAN_FRONTEND=noninteractive                                                                         && \
+       apt-get -y                                                                 update                             && \
+       apt-get -y --allow-unauthenticated -o Dpkg::Options::=--force-confdef      upgrade                            && \
+       apt-get -y --allow-unauthenticated -o Dpkg::Options::=--force-confdef dist-upgrade                            && \
+       cat /etc/os-release                                                                                           && \
+       apt-get -y install dialog apt-utils                                                                           && \
+       echo yes | apt-get -y install intel-mkl libmkl-dev                                                            && \
+       apt-get -y install build-essential                                                                            && \
+       apt-get -y install git gfortran gcc g++ openmpi-bin libopenmpi-dev automake autoconf libeigen3-dev            && \
+       cd /tmp                                                                                                       && \
+       git clone https://github.com/opencollab/arpack-ng                                                             && \
+       cd arpack-ng                                                                                                  && \
+       git checkout master                                                                                           && \
+       sed -e 's/LOG_FLAGS = /LOG_FLAGS = --allow-run-as-root --oversubscribe /' -i PARPACK/EXAMPLES/MPI/Makefile.am && \
+       sed -e 's/LOG_FLAGS = /LOG_FLAGS = --allow-run-as-root --oversubscribe /' -i PARPACK/TESTS/MPI/Makefile.am    && \
+       ./bootstrap                                                                                                   && \
+       export FFLAGS='-I/usr/include/mkl'                                                                            && \
+       export FCFLAGS='-I/usr/include/mkl'                                                                           && \
+       export LIBS='-Wl,--no-as-needed -lmkl_sequential -lmkl_core -lpthread -lm -ldl'                               && \
+       export INTERFACE64=1                                                                                          && \
+       ./configure --with-blas=mkl_gf_ilp64 --with-lapack=mkl_gf_ilp64 -enable-icb-exmm --enable-mpi                    \
+                   --disable-dependency-tracking                                                                     && \
+       export VERBOSE=1                                                                                              && \
+       make all                                                                                                      && \
+       make check"
+
+after_failure:
+
+  # show build error or test log to know what is wrong if errors occured.
+  - if [[ -f $TRAVIS_BUILD_DIR/build/Testing/Temporary/LastTest.log ]]; then tail -n 300 $TRAVIS_BUILD_DIR/build/Testing/Temporary/LastTest.log; fi
+  - find . -name test-suite.log | xargs tail -n 300

CHANGES

@@ -1,3 +1,38 @@
+arpack-ng - 3.7.0
+
+[ Franck Houssen ]
+ * [BUG FIX] ICB: missing workev for *[ds]neupd (real+not-sym) => API/ABI change for *[ds]neupd_c.
+ * [BUG FIX] autotools - make distcheck: fix circular dependencies.
+ * arpackmm: utility to test arpack with matrix market files.
+ * ICB: add ILP64 support.
+    The idea is:
+    - cmake      generates arpackdef.h from arpackdef.cmake.h.in
+    - autoheader generates arpackdef.h from arpackdef.autotools.h[.in]
+    - arpack includes arpackdef.h when/where needed:
+      - #define provide c_int/a_int according to architecture.
+      - all f90 who need to include "arpackdef.h" must be moved to F90.
+        Note: by convention, F90 are preprocessed (f90 are not).
+    - MPI does not support ILP64: integer*4 must be imposed in all
+      calls involving MPI (f90 example/test code).
+    To enable ILP64 users to compile/link, arpackdef.h is added in
+    the arpack installation (make install).
+
+ [ Kyle Guinn ]
+ * Autoconf/Automake simplifications and fixes.
+ * Simplify the generation of arpackdef.h.
+
+ -- Sylvestre Ledru <sylvestre@debian.org> Sat, 12 Jan 2019 16:24:00 +0100
+
+arpack-ng - 3.6.3
+
+[ Franck Houssen ]
+ * Add Fortran common initialization (block data).
+
+ [ Marco Caliari ]
+ * Give up forcing the initial residual to be in the range of the operator OP after a restart (Closes: #142).
+
+ -- Sylvestre Ledru <sylvestre@debian.org> Wed, 19 Sep 2018 09:59:59 +0200
+
 arpack-ng - 3.6.2
 
  * Remove all trailing whitespaces
@@ -6,7 +41,7 @@ arpack-ng - 3.6.2
  * Install: move headers into a dedicated directory (local/include/arpack).
    (Closes #126)
  * Add configuration summary.
- * Improve the flag detection. Hopefully fix the ppc64el and other archs 
+ * Improve the flag detection. Hopefully fix the ppc64el and other archs
    issues in Debian
 
  -- Sylvestre Ledru <sylvestre@debian.org> Sat, 23 Jun 2018 14:56:54 +0200

CMakeLists.txt

@@ -17,6 +17,7 @@ endif ()
 
 option(MPI "Enable parallel support" OFF)
 option(ICB "Enable support for *[ae]upd_c with ISO_C_BINDING" OFF)
+option(ICBEXMM "Enable support for matrix market example based on ICB" OFF)
 #option(SYSTEM_BLAS "Use system BLAS" ON)
 #option(SYSTEM_LAPACK "Use system LAPACK" ON)
 option(EXAMPLES "Compile ARPACK examples" OFF)
@@ -31,6 +32,9 @@ set(CMAKE_LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib)
 if (COVERALLS)
     include(Coveralls)
     set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -g -O0 -fprofile-arcs -ftest-coverage")
+    set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -g -O0 -fprofile-arcs -ftest-coverage")
+    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -O0 -fprofile-arcs -ftest-coverage")
+    # The no space is by design: issue in cmake. See CMP0004.
     set(EXTRA_LDLAGS "${EXTRA_LDLAGS}-lgcov")
 endif()
 
@@ -59,6 +63,35 @@ function(pexamples list_name)
     endforeach()
 endfunction(pexamples)
 
+if (ICBEXMM)
+    find_package(Eigen3)
+    if (NOT Eigen3_FOUND) # If not found, piggy-back pkg-config files.
+        message(WARNING "CMake didn't find the Eigen3 package. Try to look for pkg-config file...")
+        find_package(PkgConfig 3.2 REQUIRED)
+        pkg_check_modules(EIGEN3 REQUIRED eigen3>=3.2)
+    endif ()
+
+    # Look for headers.
+    find_path(EIGEN3_SPARSE_DIR NAMES Sparse PATHS ${EIGEN3_INCLUDE_DIRS} PATH_SUFFIXES Eigen)
+    if (NOT EIGEN3_SPARSE_DIR)
+        message(FATAL_ERROR "-- Eigen/Sparse header not found.")
+    endif ()
+    find_path(EIGEN3_ITERATIVE_SOLVER_DIR NAMES IterativeLinearSolvers PATHS ${EIGEN3_INCLUDE_DIRS} PATH_SUFFIXES Eigen)
+    if (NOT EIGEN3_ITERATIVE_SOLVER_DIR)
+        message(FATAL_ERROR "-- Eigen/IterativeLinearSolvers header not found.")
+    endif ()
+    find_path(EIGEN3_SLU_SOLVER_DIR NAMES SparseLU PATHS ${EIGEN3_INCLUDE_DIRS} PATH_SUFFIXES Eigen)
+    if (NOT EIGEN3_SLU_SOLVER_DIR)
+        message(FATAL_ERROR "-- Eigen/SparseLU header not found.")
+    endif ()
+    find_path(EIGEN3_SQR_SOLVER_DIR NAMES SparseQR PATHS ${EIGEN3_INCLUDE_DIRS} PATH_SUFFIXES Eigen)
+    if (NOT EIGEN3_SQR_SOLVER_DIR)
+        message(FATAL_ERROR "-- Eigen/SparseQR header not found.")
+    endif ()
+
+    set(ICB "ON")
+endif ()
+
 # Enable language(s) before any find_package (in particular before MPI find_package).
 if (ICB)
     enable_language(C CXX) # For testing binding with c/c++.
@@ -174,13 +207,13 @@ endif()
 #    file(GLOB arpacklapack_STAT_SRCS lapack/*.f)
 #endif()
 #file(GLOB arpackutil_STAT_SRCS util/*.f)
-file(GLOB arpacksrc_STAT_SRCS ${arpack_SOURCE_DIR}/SRC/*.f)
+file(GLOB arpacksrc_STAT_SRCS dbgini.f staini.f ${arpack_SOURCE_DIR}/SRC/*.f)
 
 set(arpacksrc_ICB "")
 set(parpacksrc_ICB "")
 if(ICB)
-  file(GLOB arpacksrc_ICB SRC/icba*.f90 debug_init.f90 debug_icb.f90 stat_icb.f90)
-  file(GLOB parpacksrc_ICB PARPACK/SRC/MPI/icbp*.f90 debug_init.f90 debug_icb.f90 stat_icb.f90)
+  file(GLOB arpacksrc_ICB SRC/icba*.F90 ICB/debug_icb.F90 ICB/stat_icb.F90)
+  file(GLOB parpacksrc_ICB PARPACK/SRC/MPI/icbp*.F90 ICB/debug_icb.F90 ICB/stat_icb.F90)
 endif()
 
 set(arpackutil_STAT_SRCS
@@ -202,37 +235,37 @@ set(arpackutil_STAT_SRCS
 
 if (MPI)
     file(GLOB parpackutil_STAT_SRCS PARPACK/UTIL/MPI/*.f)
-    file(GLOB parpacksrc_STAT_SRCS PARPACK/SRC/MPI/*.f)
+    file(GLOB parpacksrc_STAT_SRCS dbgini.f staini.f PARPACK/SRC/MPI/*.f)
 endif()
 
 # use -DBUILD_SHARED_LIBS=ON|OFF to control static/shared
 add_library(arpack ${arpackutil_STAT_SRCS} ${arpacksrc_STAT_SRCS} ${arpacksrc_ICB})
 
 if(ICB)
-  install(FILES arpack.h DESTINATION include/arpack)
-  install(FILES arpack.hpp DESTINATION include/arpack)
+  install(FILES ICB/arpack.h DESTINATION include/arpack)
+  install(FILES ICB/arpack.hpp DESTINATION include/arpack)
   if (MPI)
-    install(FILES parpack.h DESTINATION include/arpack)
-    install(FILES parpack.hpp DESTINATION include/arpack)
+    install(FILES ICB/parpack.h DESTINATION include/arpack)
+    install(FILES ICB/parpack.hpp DESTINATION include/arpack)
   endif()
 endif()
 
 install(FILES debug.h DESTINATION include/arpack)
 if(ICB)
-  install(FILES debug_c.h DESTINATION include/arpack)
-  install(FILES debug_c.hpp DESTINATION include/arpack)
+  install(FILES ICB/debug_c.h DESTINATION include/arpack)
+  install(FILES ICB/debug_c.hpp DESTINATION include/arpack)
 endif()
 
 install(FILES stat.h DESTINATION include/arpack)
 if(ICB)
-  install(FILES stat_c.h DESTINATION include/arpack)
-  install(FILES stat_c.hpp DESTINATION include/arpack)
+  install(FILES ICB/stat_c.h DESTINATION include/arpack)
+  install(FILES ICB/stat_c.hpp DESTINATION include/arpack)
 endif()
 
 target_link_libraries(arpack ${BLAS_LIBRARIES})
 target_link_libraries(arpack ${LAPACK_LIBRARIES} ${EXTRA_LDLAGS})
 set_target_properties(arpack PROPERTIES OUTPUT_NAME arpack${LIBSUFFIX})
-set_target_properties(arpack PROPERTIES VERSION 2.0.0)
+set_target_properties(arpack PROPERTIES VERSION 2.1.0)
 set_target_properties(arpack PROPERTIES SOVERSION 2)
 
 if (MPI)
@@ -249,11 +282,12 @@ if (MPI)
     target_link_libraries(parpack ${MPI_Fortran_LIBRARIES})
     target_link_libraries(parpack arpack)
     set_target_properties(parpack PROPERTIES OUTPUT_NAME parpack${LIBSUFFIX})
-    set_target_properties(parpack PROPERTIES VERSION 2.0.0)
+    set_target_properties(parpack PROPERTIES VERSION 2.1.0)
     set_target_properties(parpack PROPERTIES SOVERSION 2)
 endif ()
 
 if (INTERFACE64)
+    set(INTERFACE64 1)
     if ("${CMAKE_Fortran_COMPILER_ID}" MATCHES "GNU")
       set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -fdefault-integer-8")
     elseif ("${CMAKE_Fortran_COMPILER_ID}" MATCHES "Intel")
@@ -261,10 +295,15 @@ if (INTERFACE64)
     else ()
       message(WARNING "build script does not know how to make your Fortran compiler use 64-bit integers: set it manually via FFLAGS.")
     endif ()
-
-    set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -DBLASINT=int64_t")
+else ()
+    set(INTERFACE64 0)
 endif ()
 
+message("-- Generating arpackdef.h")
+configure_file(arpackdef.h.in "${PROJECT_BINARY_DIR}/arpackdef.h" @ONLY)
+include_directories("${PROJECT_BINARY_DIR}") # Find arpackdef.h
+install(FILES "${PROJECT_BINARY_DIR}/arpackdef.h" DESTINATION include/arpack)
+
 set(CMAKE_INCLUDE_CURRENT_DIR ON)
 
 ############################
@@ -483,6 +522,14 @@ add_executable(bug_79_double_complex TESTS/bug_79_double_complex.f)
 target_link_libraries(bug_79_double_complex arpack ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES} ${EXTRA_LDLAGS})
 add_test(bug_79_double_complex_tst Tests/bug_79_double_complex)
 
+add_executable(bug_142 TESTS/bug_142.f)
+target_link_libraries(bug_142 arpack ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES} ${EXTRA_LDLAGS})
+add_test(bug_142_tst Tests/bug_142)
+
+add_executable(bug_142_gen TESTS/bug_142_gen.f)
+target_link_libraries(bug_142_gen arpack ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES} ${EXTRA_LDLAGS})
+add_test(bug_142_gen_tst Tests/bug_142_gen)
+
 if(MPI)
   add_executable(issue46 PARPACK/TESTS/MPI/issue46.f)
   target_link_libraries(issue46 parpack arpack ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES} ${EXTRA_LDLAGS})
@@ -491,23 +538,36 @@ endif()
 
 if(ICB)
   add_executable(icb_arpack_c TESTS/icb_arpack_c.c)
-  target_include_directories(icb_arpack_c PUBLIC ${PROJECT_SOURCE_DIR}) # Get arpack.h
+  target_include_directories(icb_arpack_c PUBLIC ${PROJECT_SOURCE_DIR}/ICB) # Get arpack.h
   target_link_libraries(icb_arpack_c arpack ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES} ${EXTRA_LDLAGS})
   add_test(icb_arpack_c_tst Tests/icb_arpack_c)
 
   add_executable(icb_arpack_cpp TESTS/icb_arpack_cpp.cpp)
-  target_include_directories(icb_arpack_cpp PUBLIC ${PROJECT_SOURCE_DIR}) # Get arpack.hpp
+  target_include_directories(icb_arpack_cpp PUBLIC ${PROJECT_SOURCE_DIR}/ICB) # Get arpack.hpp
   target_link_libraries(icb_arpack_cpp arpack ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES} ${EXTRA_LDLAGS})
   add_test(icb_arpack_cpp_tst Tests/icb_arpack_cpp)
 
+  if (ICBEXMM)
+    add_executable(arpackmm EXAMPLES/MATRIX_MARKET/arpackmm.cpp)
+    target_include_directories(arpackmm PUBLIC ${PROJECT_SOURCE_DIR}/ICB ${EIGEN3_INCLUDE_DIRS}) # Get arpack.h + eigen
+    target_link_libraries(arpackmm arpack ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES} ${EXTRA_LDLAGS})
+    configure_file(EXAMPLES/MATRIX_MARKET/As.mtx Tests/As.mtx)
+    configure_file(EXAMPLES/MATRIX_MARKET/An.mtx Tests/An.mtx)
+    configure_file(EXAMPLES/MATRIX_MARKET/Az.mtx Tests/Az.mtx)
+    configure_file(EXAMPLES/MATRIX_MARKET/B.mtx  Tests/B.mtx)
+    configure_file(EXAMPLES/MATRIX_MARKET/Bz.mtx Tests/Bz.mtx)
+    configure_file(EXAMPLES/MATRIX_MARKET/arpackmm.sh Tests/arpackmm.sh)
+    add_test(NAME arpackmm_tst WORKING_DIRECTORY ${arpack_BINARY_DIR}/Tests COMMAND arpackmm.sh)
+  endif()
+
   if (MPI)
     add_executable(icb_parpack_c PARPACK/TESTS/MPI/icb_parpack_c.c)
-    target_include_directories(icb_parpack_c PUBLIC ${PROJECT_SOURCE_DIR} ${MPI_C_INCLUDE_DIRS}) # Get parpack.h mpi.h
+    target_include_directories(icb_parpack_c PUBLIC ${PROJECT_SOURCE_DIR}/ICB ${MPI_C_INCLUDE_DIRS}) # Get parpack.h mpi.h
     target_link_libraries(icb_parpack_c parpack arpack ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES} ${EXTRA_LDLAGS} ${MPI_C_LIBRARIES})
     add_test(icb_parpack_c_tst mpirun -n 2 Tests/icb_parpack_c)
 
     add_executable(icb_parpack_cpp PARPACK/TESTS/MPI/icb_parpack_cpp.cpp)
-    target_include_directories(icb_parpack_cpp PUBLIC ${PROJECT_SOURCE_DIR} ${MPI_CXX_INCLUDE_DIRS}) # Get parpack.hpp mpi.h
+    target_include_directories(icb_parpack_cpp PUBLIC ${PROJECT_SOURCE_DIR}/ICB ${MPI_CXX_INCLUDE_DIRS}) # Get parpack.hpp mpi.h
     target_link_libraries(icb_parpack_cpp parpack arpack ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES} ${EXTRA_LDLAGS} ${MPI_CXX_LIBRARIES})
     add_test(icb_parpack_cpp_tst mpirun -n 2 Tests/icb_parpack_cpp)
   endif()
@@ -519,21 +579,6 @@ endif()
 # 'make install' to the correct location
 include(GNUInstallDirs)
 
-# Convert variable names to those expected by the .pc file.
-set(prefix ${CMAKE_INSTALL_PREFIX})
-set(exec_prefix \${prefix})
-set(libdir \${exec_prefix}/${CMAKE_INSTALL_LIBDIR})
-set(PACKAGE_NAME ${PROJECT_NAME})
-# TODO: Fill these in with something appropriate.
-#set(PACKAGE_VERSION)
-#set(PACKAGE_URL)
-#set(LAPACK_LIBS)
-#set(BLAS_LIBS)
-configure_file(arpack.pc.in arpack${LIBSUFFIX}.pc @ONLY)
-
-#install(FILES ${CMAKE_CURRENT_BINARY_DIR}/arpack${LIBSUFFIX}.pc
-#    DESTINATION ${CMAKE_INSTALL_LIBDIR}/pkgconfig)
-
 install(TARGETS arpack
     ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR}
     LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}
@@ -550,7 +595,7 @@ endif ()
 configure_file(arpack-ng-config.cmake.in "${PROJECT_BINARY_DIR}/arpack-ng-config.cmake" @ONLY)
 install(FILES "${PROJECT_BINARY_DIR}/arpack-ng-config.cmake" DESTINATION "${CMAKE_INSTALL_PREFIX}/lib/cmake") # find_package(arpack-ng)
 set(arpack_ng_MAJOR_VERSION 3)
-set(arpack_ng_MINOR_VERSION 6)
+set(arpack_ng_MINOR_VERSION 7)
 set(arpack_ng_PATCH_VERSION 0)
 set(arpack_ng_VERSION ${arpack_ng_MAJOR_VERSION}.${arpack_ng_MINOR_VERSION}.${arpack_ng_PATCH_VERSION})
 configure_file(arpack-ng-config-version.cmake.in "${PROJECT_BINARY_DIR}/arpack-ng-config-version.cmake" @ONLY)
@@ -564,6 +609,21 @@ set(CPACK_SOURCE_IGNORE_FILES "/BUILD/" "/Build/" "/build/" "/local/") # Do not
 set(CPACK_SOURCE_PACKAGE_FILE_NAME "arpack-ng-${CPACK_VERSION_MAJOR}.${CPACK_VERSION_MINOR}.${CPACK_VERSION_PATCH}")
 include(CPack)
 
+# Convert variable names to those expected by the .pc file.
+set(prefix ${CMAKE_INSTALL_PREFIX})
+set(exec_prefix \${prefix})
+set(libdir \${exec_prefix}/${CMAKE_INSTALL_LIBDIR})
+set(includedir \${prefix}/${CMAKE_INSTALL_INCLUDEDIR})
+set(PACKAGE_NAME ${PROJECT_NAME})
+set(PACKAGE_VERSION ${arpack_ng_VERSION})
+set(PACKAGE_URL "https://github.com/opencollab/arpack-ng/")
+# TODO: Fill these in with something appropriate.
+#set(LAPACK_LIBS)
+#set(BLAS_LIBS)
+configure_file(arpack.pc.in arpack${LIBSUFFIX}.pc @ONLY)
+#install(FILES ${CMAKE_CURRENT_BINARY_DIR}/arpack${LIBSUFFIX}.pc
+#    DESTINATION ${CMAKE_INSTALL_LIBDIR}/pkgconfig)
+
 
 # We don't want this to run on every build.
 option(COVERALLS "Generate coveralls data" OFF)
@@ -573,8 +633,9 @@ if (COVERALLS)
     set(COVERAGE_SRCS awesome.c code.c files.c)
 
     # Create the coveralls target.
+    # Also lists the c/cpp files for test purposes
     coveralls_setup(
-        "${arpackutil_STAT_SRCS} ${arpacksrc_STAT_SRCS} ${arpacksrc_ICB}" # The source files.
+        "${arpackutil_STAT_SRCS} ${arpacksrc_STAT_SRCS} ${arpacksrc_ICB} ${parpacksrc_DIR} ${parpackutil_DIR} ${PROJECT_SOURCE_DIR}/TESTS/icb_arpack_c.c ${PROJECT_SOURCE_DIR}/TESTS/icb_arpack_cpp.cpp" # The source files.
         ON                 # If we should upload.
         "${PROJECT_SOURCE_DIR}/cmake/") # (Optional) Alternate project cmake module path.
 endif()
@@ -609,6 +670,7 @@ message("-- Configuration summary for arpack-ng-${arpack_ng_VERSION}:")
 message("   -- prefix: ${CMAKE_INSTALL_PREFIX}")
 message("   -- MPI: ${MPI}")
 message("   -- ICB: ${ICB}")
+message("   -- INTERFACE64: ${INTERFACE64}")
 cprsummary("FC" "${CMAKE_Fortran_COMPILER}"
                 "${CMAKE_Fortran_FLAGS_DEBUG}"
                 "${CMAKE_Fortran_FLAGS_MINSIZEREL}"
@@ -638,3 +700,6 @@ if (MPI)
 endif()
 libsummary("BLAS" "" "${BLAS_LIBRARIES}")
 libsummary("LAPACK" "" "${LAPACK_LIBRARIES}")
+if (ICBEXMM)
+    libsummary("EIGEN3" "${EIGEN3_INCLUDE_DIRS}" "")
+endif()

EXAMPLES/MATRIX_MARKET/An.mtx

@@ -0,0 +1,79 @@
+%% MatrixMarket matrix coordinate double symmetric
+% This example is 0-based without (optional) nnz
+%
+% This is a 1D diffusive laplacian matrix (fixed at first end <=> invertible)
+%
+%              1                              1          
+%              .                             .--.        
+%             / \                            |  |        
+%            /   \ 0                       0 |  |        
+%  phi_i  --o  o  o--o--  =>  grad(phi_i)  --o  o  o--o--
+%              i  j                            i|  |j    
+%                                               |  |     
+%                                               .--.     
+%                                                -1.     
+%
+%                 1                              1       
+%                 .                             .--.     
+%                / \                            |  |     
+%               /   \ 0                        i|  |j    
+%  phi_j  --o--o  o  o--  =>  grad(phi_j)  --o--o  o  o--
+%              i  j                                |  | 0
+%                                                  |  |  
+%                                                  .--.  
+%                                                   -1.  
+%
+%                  i     j
+%              | l_ii  l_ij | i
+%  laplacian = |            |
+%              | l_ji  l_jj | j
+%
+%  distance(i, j) = d = 1.
+%
+%  l_ii = int_[i,j](grad(phi_i).grad(phi_i)) = d*(-1.)*(-1.) =  1.
+%  l_ij = int_[i,j](grad(phi_i).grad(phi_j)) = d*(-1.)*( 1.) = -1.
+%  l_ji = int_[i,j](grad(phi_j).grad(phi_i)) = d*( 1.)*(-1.) = -1.
+%  l_jj = int_[i,j](grad(phi_j).grad(phi_j)) = d*( 1.)*( 1.) =  1.
+%
+%                i     j
+%              | d_ii  d_ij | i
+%  diffusion = |            |
+%              | d_ji  d_jj | j
+%
+%  d_ii = int_[i,j](phi_i.grad(phi_i)) = int_[i,j]((1-x)*(-1.)) = -d*0.5  = -0.5
+%  d_ij = int_[i,j](phi_i.grad(phi_j)) = int_[i,j]((1-x)*( 1.)) =  d*0.5  =  0.5
+%  d_ji = int_[i,j](phi_j.grad(phi_i)) = int_[i,j](   x *(-1.)) = -d*0.5  = -0.5
+%  d_jj = int_[i,j](phi_j.grad(phi_j)) = int_[i,j](   x *( 1.)) =  d*0.5  =  0.5
+%
+% A <=> assembly of {kappa*laplacian + rho*diffusion}
+%       where kappa = 100 and rho = 2
+%
+% n m [nnz]
+% i j Aij
+
+8 8
+
+0  0    1.
+1  1  200.
+2  2  200.
+3  3  200.
+4  4  200.
+5  5  200.
+6  6  200.
+7  7  101.
+
+1  0     0.
+2  1  -101.
+3  2  -101.
+4  3  -101.
+5  4  -101.
+6  5  -101.
+7  6  -101.
+
+0  1    0.
+1  2  -99.
+2  3  -99.
+3  4  -99.
+4  5  -99.
+5  6  -99.
+6  7  -99.

EXAMPLES/MATRIX_MARKET/As.mtx

@@ -0,0 +1,69 @@
+%% MatrixMarket matrix coordinate double symmetric
+% This example is 0-based without (optional) nnz
+%
+% This is a 1D laplacian matrix (fixed at first end <=> invertible)
+%
+%              1                              1          
+%              .                             .--.        
+%             / \                            |  |        
+%            /   \ 0                       0 |  |        
+%  phi_i  --o  o  o--o--  =>  grad(phi_i)  --o  o  o--o--
+%              i  j                            i|  |j    
+%                                               |  |     
+%                                               .--.     
+%                                                -1.     
+%
+%                 1                              1       
+%                 .                             .--.     
+%                / \                            |  |     
+%               /   \ 0                        i|  |j    
+%  phi_j  --o--o  o  o--  =>  grad(phi_j)  --o--o  o  o--
+%              i  j                                |  | 0
+%                                                  |  |  
+%                                                  .--.  
+%                                                   -1.  
+%
+%                  i     j
+%              | l_ii  l_ij | i
+%  laplacian = |            |
+%              | l_ji  l_jj | j
+%
+%  distance(i, j) = d = 1.
+%
+%  l_ii = int_[i,j](grad(phi_i).grad(phi_i)) = d*(-1.)*(-1.) =  1.
+%  l_ij = int_[i,j](grad(phi_i).grad(phi_j)) = d*(-1.)*( 1.) = -1.
+%  l_ji = int_[i,j](grad(phi_j).grad(phi_i)) = d*( 1.)*(-1.) = -1.
+%  l_jj = int_[i,j](grad(phi_j).grad(phi_j)) = d*( 1.)*( 1.) =  1.
+%
+% A <=> assembly of {kappa*laplacian}
+%       where kappa = 100
+%
+% n m [nnz]
+% i j Aij
+
+8 8
+
+0  0    1.
+1  1  200.
+2  2  200.
+3  3  200.
+4  4  200.
+5  5  200.
+6  6  200.
+7  7  100.
+
+1  0     0.00
+2  1  -100.00
+3  2  -100.00
+4  3  -100.00
+5  4  -100.00
+6  5  -100.00
+7  6  -100.00
+
+0  1     0.00
+1  2  -100.00
+2  3  -100.00
+3  4  -100.00
+4  5  -100.00
+5  6  -100.00
+6  7  -100.00

EXAMPLES/MATRIX_MARKET/Az.mtx

@@ -0,0 +1,79 @@
+%% MatrixMarket matrix coordinate double symmetric
+% This example is 0-based without (optional) nnz
+%
+% This is a 1D diffusive laplacian matrix (fixed at first end <=> invertible)
+%
+%              1                              1          
+%              .                             .--.        
+%             / \                            |  |        
+%            /   \ 0                       0 |  |        
+%  phi_i  --o  o  o--o--  =>  grad(phi_i)  --o  o  o--o--
+%              i  j                            i|  |j    
+%                                               |  |     
+%                                               .--.     
+%                                                -1.     
+%
+%                 1                              1       
+%                 .                             .--.     
+%                / \                            |  |     
+%               /   \ 0                        i|  |j    
+%  phi_j  --o--o  o  o--  =>  grad(phi_j)  --o--o  o  o--
+%              i  j                                |  | 0
+%                                                  |  |  
+%                                                  .--.  
+%                                                   -1.  
+%
+%                  i     j
+%              | l_ii  l_ij | i
+%  laplacian = |            |
+%              | l_ji  l_jj | j
+%
+%  distance(i, j) = d = 1.
+%
+%  l_ii = int_[i,j](grad(phi_i).grad(phi_i)) = d*(-1.)*(-1.) =  1.
+%  l_ij = int_[i,j](grad(phi_i).grad(phi_j)) = d*(-1.)*( 1.) = -1.
+%  l_ji = int_[i,j](grad(phi_j).grad(phi_i)) = d*( 1.)*(-1.) = -1.
+%  l_jj = int_[i,j](grad(phi_j).grad(phi_j)) = d*( 1.)*( 1.) =  1.
+%
+%                i     j
+%              | d_ii  d_ij | i
+%  diffusion = |            |
+%              | d_ji  d_jj | j
+%
+%  d_ii = int_[i,j](phi_i.grad(phi_i)) = int_[i,j]((1-x)*(-1.)) = -d*0.5  = -0.5
+%  d_ij = int_[i,j](phi_i.grad(phi_j)) = int_[i,j]((1-x)*( 1.)) =  d*0.5  =  0.5
+%  d_ji = int_[i,j](phi_j.grad(phi_i)) = int_[i,j](   x *(-1.)) = -d*0.5  = -0.5
+%  d_jj = int_[i,j](phi_j.grad(phi_j)) = int_[i,j](   x *( 1.)) =  d*0.5  =  0.5
+%
+% A <=> assembly of {kappa*laplacian + rho*diffusion}
+%       where kappa = 100 and rho = 2
+%
+% n m [nnz]
+% i j Aij
+
+8 8
+
+0  0  (  1., 0.)
+1  1  (200., 0.)
+2  2  (200., 0.)
+3  3  (200., 0.)
+4  4  (200., 0.)
+5  5  (200., 0.)
+6  6  (200., 0.)
+7  7  (101., 0.)
+
+1  0  (   0., 0.)
+2  1  (-101., 0.)
+3  2  (-101., 0.)
+4  3  (-101., 0.)
+5  4  (-101., 0.)
+6  5  (-101., 0.)
+7  6  (-101., 0.)
+
+0  1  (  0., 0.)
+1  2  (-99., 0.)
+2  3  (-99., 0.)
+3  4  (-99., 0.)
+4  5  (-99., 0.)
+5  6  (-99., 0.)
+6  7  (-99., 0.)

EXAMPLES/MATRIX_MARKET/B.mtx

@@ -0,0 +1,62 @@
+%% MatrixMarket matrix coordinate double general
+% This example is 1-based with (optional) nnz
+%
+% This is a 1D mass matrix
+%
+%              1        
+%              .        
+%             / \       
+%            /   \ 0    
+%  phi_i  --o  o  o--o--
+%              i  j     
+%
+%                 1     
+%                 .     
+%                / \    
+%               /   \ 0 
+%  phi_j  --o--o  o  o--
+%              i  j     
+%
+%                i     j
+%            | i_ii  i_ij | i
+%  inertia = |            |
+%            | i_ji  i_jj | j
+%
+%  distance(i, j) = d = 1.
+%
+%  i_ii = int_[i,j](phi_i.phi_i) = int_[i,j]((1-x)(1-x)) = i_jj (area under the curve)
+%  i_ij = int_[i,j](phi_i.phi_j) = int_[i,j]((1-x)   x ) = i_ji
+%  i_ji = int_[i,j](phi_j.phi_i) = int_[i,j](   x (1-x)) = d*1./6.
+%  i_jj = int_[i,j](phi_j.phi_j) = int_[i,j](   x    x ) = d*1./3.
+%
+% B <=> mass assembly
+%
+% n m [nnz]
+% i j Bij
+
+8 8 22
+
+1  1  0.333
+2  2  0.333
+3  3  0.333
+4  4  0.333
+5  5  0.333
+6  6  0.333
+7  7  0.333
+8  8  0.333
+
+2  1  0.166
+3  2  0.166
+4  3  0.166
+5  4  0.166
+6  5  0.166
+7  6  0.166
+8  7  0.166
+
+1  2  0.166
+2  3  0.166
+3  4  0.166
+4  5  0.166
+5  6  0.166
+6  7  0.166
+7  8  0.166

EXAMPLES/MATRIX_MARKET/Bz.mtx

@@ -0,0 +1,62 @@
+%% MatrixMarket matrix coordinate double general
+% This example is 1-based with (optional) nnz
+%
+% This is a 1D mass matrix
+%
+%              1        
+%              .        
+%             / \       
+%            /   \ 0    
+%  phi_i  --o  o  o--o--
+%              i  j     
+%
+%                 1     
+%                 .     
+%                / \    
+%               /   \ 0 
+%  phi_j  --o--o  o  o--
+%              i  j     
+%
+%                i     j
+%            | i_ii  i_ij | i
+%  inertia = |            |
+%            | i_ji  i_jj | j
+%
+%  distance(i, j) = d = 1.
+%
+%  i_ii = int_[i,j](phi_i.phi_i) = int_[i,j]((1-x)(1-x)) = i_jj (area under the curve)
+%  i_ij = int_[i,j](phi_i.phi_j) = int_[i,j]((1-x)   x ) = i_ji
+%  i_ji = int_[i,j](phi_j.phi_i) = int_[i,j](   x (1-x)) = d*1./6.
+%  i_jj = int_[i,j](phi_j.phi_j) = int_[i,j](   x    x ) = d*1./3.
+%
+% B <=> mass assembly
+%
+% n m [nnz]
+% i j Bij
+
+8 8 22
+
+1  1  (0.333, 0.)
+2  2  (0.333, 0.)
+3  3  (0.333, 0.)
+4  4  (0.333, 0.)
+5  5  (0.333, 0.)
+6  6  (0.333, 0.)
+7  7  (0.333, 0.)
+8  8  (0.333, 0.)
+
+2  1  (0.166, 0.)
+3  2  (0.166, 0.)
+4  3  (0.166, 0.)
+5  4  (0.166, 0.)
+6  5  (0.166, 0.)
+7  6  (0.166, 0.)
+8  7  (0.166, 0.)
+
+1  2  (0.166, 0.)
+2  3  (0.166, 0.)
+3  4  (0.166, 0.)
+4  5  (0.166, 0.)
+5  6  (0.166, 0.)
+6  7  (0.166, 0.)
+7  8  (0.166, 0.)

EXAMPLES/MATRIX_MARKET/Makefile.am

@@ -0,0 +1,13 @@
+LDADD = $(top_builddir)/SRC/libarpack$(LIBSUFFIX).la
+AM_CPPFLAGS = -I$(top_builddir) -I$(top_srcdir)/ICB $(EIGEN3_CFLAGS)
+
+EXTRA_DIST = README
+
+check_PROGRAMS = arpackmm
+dist_check_SCRIPTS = arpackmm.sh
+dist_check_DATA = An.mtx As.mtx Az.mtx B.mtx Bz.mtx
+TESTS = arpackmm.sh
+
+arpackmm_SOURCES = arpackmm.cpp
+
+CLEANFILES = resid.out v.out

EXAMPLES/MATRIX_MARKET/README

@@ -0,0 +1,22 @@
+arpackmm: utility to test arpack with matrix market files.
+
+This is meant to run arpack with different options to find
+what are the ones that are the best for your particular
+problem.
+
+Typically: computing small eigen values may breakdown and/or
+be slow (arpack is good at finding large eigen values only).
+To compute small eigen values, it may be a better choice
+to look for large eigen values with invert or shift+invert.
+
+If solve breaks down, you may increase --nbCV, try --restart,
+play with --shiftReal/Imag and/or --invert.
+
+To build this utility, you need:
+- to use a fortran compiler which supports iso_c_binding.
+- to have installed eigen3 (to deal with the RCI).
+- to configure arpack-ng this way:
+  - autotools:
+    ~arpack-ng> ./configure --enable-icb-exmm; make all check
+  - cmake:
+    ~arpack-ng/build> cmake -D ICBEXMM=ON ..; make all test

EXAMPLES/MATRIX_MARKET/arpackmm.cpp

@@ -0,0 +1,963 @@
+// This code sample is meant for convenience (not performance):
+//   - test/run arpack (eigen values / vectors, timing).
+//   - play with modes: shift, invert, shift + invert.
+//   - use with user matrices (matrix market format).
+
+#include <iostream>
+#include <string>
+#include <sstream> // stringstream.
+#include <fstream> // [io]fstream.
+#include <vector>
+#include <complex>
+#include <algorithm> // max_element.
+#include <chrono>
+#include <limits> // epsilon.
+#include <cmath> // fabs.
+#include <iomanip> // setw.
+#include <cassert> // assert.
+#include <memory> // unique_ptr.
+#include "arpack.h"
+#include "debug_c.hpp"
+
+#include <Eigen/Sparse>
+#include <Eigen/IterativeLinearSolvers>
+#include <Eigen/SparseLU>
+#include <Eigen/SparseQR>
+
+using namespace std;
+
+typedef Eigen::SparseMatrix<        double>                     EigMatR;  // Real.
+typedef Eigen::Triplet     <        double>                     EigCooR;  // Real.
+typedef Eigen::SparseMatrix<complex<double>>                    EigMatC;  // Complex.
+typedef Eigen::Triplet     <complex<double>>                    EigCooC;  // Complex.
+typedef Eigen::Matrix      <        double,  Eigen::Dynamic, 1> EigVecR;  // Real.
+typedef Eigen::Map         <EigVecR>                            EigMpVR;  // Real.
+typedef Eigen::Matrix      <complex<double>, Eigen::Dynamic, 1> EigVecC;  // Complex.
+typedef Eigen::Map         <EigVecC>                            EigMpVC;  // Complex.
+typedef Eigen::BiCGSTAB         <EigMatR>                       EigBiCGR; // Real.
+typedef Eigen::ConjugateGradient<EigMatR>                       EigCGR;   // Real.
+typedef Eigen::SparseLU<EigMatR, Eigen::COLAMDOrdering<int>>    EigSLUR;  // Real.
+typedef Eigen::SparseQR<EigMatR, Eigen::COLAMDOrdering<int>>    EigSQRR;  // Real.
+typedef Eigen::BiCGSTAB         <EigMatC>                       EigBiCGC; // Complex.
+typedef Eigen::ConjugateGradient<EigMatC>                       EigCGC;   // Complex.
+typedef Eigen::SparseLU<EigMatC, Eigen::COLAMDOrdering<int>>    EigSLUC;  // Complex.
+typedef Eigen::SparseQR<EigMatC, Eigen::COLAMDOrdering<int>>    EigSQRC;  // Complex.
+
+class options {
+  public:
+    options() {
+      fileA = "A.mtx";
+      fileB = "N.A."; // Not available.
+      nbEV = 1;
+      nbCV = 2*nbEV + 1;
+      stdPb = true; // Standard or generalized (= not standard).
+      symPb = true;
+      cpxPb = false;
+      mag = string("LM"); // Large magnitude.
+      shiftReal = false; shiftImag = false;
+      sigmaReal = 0.; sigmaImag = 0.; // Eigen value translation: look for lambda+sigma instead of lambda.
+      invert = false; // Eigen value invertion: look for 1./lambda instead of lambda.
+      tol = 1.e-06;
+      maxIt = 100;
+      schur = false; // Compute Ritz vectors.
+      slv = "BiCG";
+      slvItrTol = nullptr;
+      slvItrMaxIt = nullptr;
+      slvDrtPvtThd = nullptr;
+      check = true;
+      verbose = 0;
+      debug = 0;
+      restart = false;
+    };
+
+    int readCmdLine(int argc, char ** argv) {
+      // Check for command line independent parameters.
+
+      for (int a = 1; argv && a < argc; a++) {
+        string clo = argv[a]; // Command line option.
+        if (clo == "--help") return usage(0);
+        if (clo == "--A") {
+          a++; if (a >= argc) {cerr << "Error: bad " << clo << " - need argument" << endl; return usage();}
+          fileA = argv[a];
+        }
+        if (clo == "--nbEV") {
+          a++; if (a >= argc) {cerr << "Error: bad " << clo << " - need argument" << endl; return usage();}
+          stringstream nEV(argv[a]);
+          nEV >> nbEV; if (!nEV) {cerr << "Error: bad " << clo << " - bad argument" << endl; return usage();}
+          nbCV = 2*nbEV + 1;
+        }
+        if (clo == "--genPb") {
+          stdPb = false;
+          fileB = "B.mtx";
+        }
+        if (clo == "--nonSymPb") symPb = false;
+        if (clo == "--cpxPb") {
+          symPb = false;
+          cpxPb = true;
+        }
+        if (clo == "--mag") {
+          a++; if (a >= argc) {cerr << "Error: bad " << clo << " - need argument" << endl; return usage();}
+          mag = argv[a]; // small mag (likely poor perf) <=> large mag + invert (likely good perf).
+          bool ok = (mag == "LM" || mag == "SM" || mag == "LR" || mag == "SR" || mag == "LI" || mag == "SI") ? true : false;
+          if (!ok) {cerr << "Error: bad " << clo << " - bad argument" << endl; return usage();}
+        }
+        if (clo == "--shiftReal") {
+          shiftReal = true;
+          a++; if (a >= argc) {cerr << "Error: bad " << clo << " - need argument" << endl; return usage();}
+          stringstream s(argv[a]);
+          s >> sigmaReal; if (!s) {cerr << "Error: bad " << clo << " - bad argument" << endl; return usage();}
+        }
+        if (clo == "--shiftImag") {
+          shiftImag = true;
+          a++; if (a >= argc) {cerr << "Error: bad " << clo << " - need argument" << endl; return usage();}
+          stringstream s(argv[a]);
+          s >> sigmaImag; if (!s) {cerr << "Error: bad " << clo << " - bad argument" << endl; return usage();}
+        }
+        if (clo == "--invert") invert = true;
+        if (clo == "--tol") {
+          a++; if (a >= argc) {cerr << "Error: bad " << clo << " - need argument" << endl; return usage();}
+          stringstream t(argv[a]);
+          t >> tol; if (!t) {cerr << "Error: bad " << clo << " - bad argument" << endl; return usage();}
+        }
+        if (clo == "--maxIt") {
+          a++; if (a >= argc) {cerr << "Error: bad " << clo << " - need argument" << endl; return usage();}
+          stringstream mi(argv[a]);
+          mi >> maxIt; if (!mi) {cerr << "Error: bad " << clo << " - bad argument" << endl; return usage();}
+        }
+        if (clo == "--slv") {
+          a++; if (a >= argc) {cerr << "Error: bad " << clo << " - need argument" << endl; return usage();}
+          slv = argv[a];
+        }
+        if (clo == "--slvItrTol") {
+          a++; if (a >= argc) {cerr << "Error: bad " << clo << " - need argument" << endl; return usage();}
+          stringstream t(argv[a]);
+          double tol = 0.;
+          t >> tol; if (!t) {cerr << "Error: bad " << clo << " - bad argument" << endl; return usage();}
+          slvItrTol = unique_ptr<double>(new double);
+          if (slvItrTol) *slvItrTol = tol;
+        }
+        if (clo == "--slvItrMaxIt") {
+          a++; if (a >= argc) {cerr << "Error: bad " << clo << " - need argument" << endl; return usage();}
+          stringstream mi(argv[a]);
+          int maxIt = 0;
+          mi >> maxIt; if (!mi) {cerr << "Error: bad " << clo << " - bad argument" << endl; return usage();}
+          slvItrMaxIt = unique_ptr<int>(new int);
+          if (slvItrMaxIt) *slvItrMaxIt = maxIt;
+        }
+        if (clo == "--slvDrtPvtThd") {
+          a++; if (a >= argc) {cerr << "Error: bad " << clo << " - need argument" << endl; return usage();}
+          stringstream t(argv[a]);
+          double thd = 0.;
+          t >> thd; if (!t) {cerr << "Error: bad " << clo << " - bad argument" << endl; return usage();}
+          slvDrtPvtThd = unique_ptr<double>(new double);
+          if (slvDrtPvtThd) *slvDrtPvtThd = thd;
+        }
+        if (clo == "--noCheck") check = false;
+        if (clo == "--verbose") {
+          a++; if (a >= argc) {cerr << "Error: bad " << clo << " - need argument" << endl; return usage();}
+          stringstream vb(argv[a]);
+          vb >> verbose; if (!vb) {cerr << "Error: bad " << clo << " - bad argument" << endl; return usage();}
+        }
+        if (clo == "--debug") {
+          a++; if (a >= argc) {cerr << "Error: bad " << clo << " - need argument" << endl; return usage();}
+          stringstream dbg(argv[a]);
+          dbg >> debug; if (!dbg) {cerr << "Error: bad " << clo << " - bad argument" << endl; return usage();}
+          if (debug > 3) debug = 3;
+          debug_c(6, -6, debug, debug, debug, debug, debug, debug, debug, debug, debug, debug, debug,
+                  debug, debug, debug, debug, debug, debug, debug, debug, debug, debug, debug);
+        }
+        if (clo == "--restart") restart = true;
+      }
+
+      // Check for command line dependent parameters.
+
+      for (int a = 1; argv && a < argc; a++) {
+        string clo = argv[a]; // Command line option.
+        if (clo == "--nbCV") {
+          a++; if (a >= argc) {cerr << "Error: bad " << clo << " - need argument" << endl; return usage();}
+          stringstream nCV(argv[a]);
+          nCV >> nbCV; if (!nCV) {cerr << "Error: bad " << clo << " - bad argument" << endl; return usage();}
+        }
+        if (clo == "--B") {
+          a++; if (a >= argc) {cerr << "Error: bad " << clo << " - need argument" << endl; return usage();}
+          fileB = argv[a];
+        }
+      }
+
+      return 0;
+    };
+
+    int usage(int rc = 1) {
+      cout << "Usage: running arpack to check for eigen values/vectors." << endl;
+      cout << endl;
+      cout << "  --A F:            file name of matrix A such that A X = lambda X. (standard)" << endl;
+      cout << "                    default: A.mtx" << endl;
+      cout << "  --B F:            file name of matrix B such that A X = lambda B X. (generalized)" << endl;
+      cout << "                    default: N.A. for standard problem, or, B.mtx for generalized problem" << endl;
+      cout << "  --nbEV:           number of eigen values/vectors to compute." << endl;
+      cout << "                    default: 1" << endl;
+      cout << "  --nbCV:           number of columns of the matrix V." << endl;
+      cout << "                    default: 2*nbEV+1" << endl;
+      cout << "  --genPb:          generalized problem." << endl;
+      cout << "                    default: standard problem" << endl;
+      cout << "  --nonSymPb:       non symmetric problem (<=> use dn[ae]upd)." << endl;
+      cout << "                    default: symmetric problem (<=> use ds[ae]upd)" << endl;
+      cout << "  --cpxPb:          complex (non symmetric) problem (<=> use zn[ae]upd)." << endl;
+      cout << "                    default: false (<=> use d*[ae]upd)" << endl;
+      cout << "  --mag M:          set magnitude of eigen values to look for (LM, SM, LR, SR, LI, SI)." << endl;
+      cout << "                    default: large magnitude (LM)" << endl;
+      cout << "  --shiftReal S:    real shift where sigma = S (look for lambda+S instead of lambda)." << endl;
+      cout << "                    default: no shift, S = 0." << endl;
+      cout << "  --shiftImag S:    imaginary shift where sigma = S (look for lambda+S instead of lambda)." << endl;
+      cout << "                    default: no shift, S = 0." << endl;
+      cout << "  --invert:         invert mode (look for 1./lambda instead of lambda)." << endl;
+      cout << "                    default: no invert" << endl;
+      cout << "  --tol T:          tolerance T." << endl;
+      cout << "                    default: 1.e-06" << endl;
+      cout << "  --maxIt M:        maximum iterations M." << endl;
+      cout << "                    default: 100" << endl;
+      cout << "  --schur:          compute Schur vectors." << endl;
+      cout << "                      the Schur decomposition is such that A = Q^H x T x Q where:" << endl;
+      cout << "                        - the H superscript refers to the Hermitian transpose: Q^H = (Q^t)^*." << endl;
+      cout << "                        - Q is unitary: Q is such that Q^H x Q = I." << endl;
+      cout << "                        - T is an upper-triangular matrix whose diagonal elements are the eigenvalues of A." << endl;
+      cout << "                      every square matrix has a Schur decomposition: columns of Q are the Schur vectors." << endl;
+      cout << "                      for a general matrix A, there is no relation between Schur vectors of A and eigenvectors of A." << endl;
+      cout << "                      if q_j is the i-th Schur vector, then A x q_j is a linear combination of q_1, ..., q_j." << endl;
+      cout << "                      Schur vectors q_1, q_2, ..., q_j span an invariant subspace of A." << endl;
+      cout << "                      the Schur vectors and eigenvectors of A are the same if A is a normal matrix." << endl;
+      cout << "                    default: compute Ritz vectors (approximations of eigen vectors)" << endl;
+      cout << "  --slv S:          solver (BiCG, CG, LU)" << endl;
+      cout << "                      BiCG: iterative method, any matrices" << endl;
+      cout << "                        CG: iterative method, sym matrices only" << endl;
+      cout << "                        LU:    direct method, any matrices" << endl;
+      cout << "                        QR:    direct method, any matrices" << endl;
+      cout << "                    default: BiCG" << endl;
+      cout << "  --slvItrTol T:    solver tolerance T (for iterative solvers)." << endl;
+      cout << "                    default: eigen default value" << endl;
+      cout << "  --slvItrMaxIt M:  solver maximum iterations M (for iterative solvers)." << endl;
+      cout << "                    default: eigen default value" << endl;
+      cout << "  --slvDrtPvtThd T: solver pivot threshold T (for direct solvers)." << endl;
+      cout << "                    default: eigen default value" << endl;
+      cout << "  --noCheck:        check arpack eigen values/vectors." << endl;
+      cout << "                    check will fail if Schur vectors are computed and A is NOT a normal matrix." << endl;
+      cout << "                    default: check" << endl;
+      cout << "  --verbose V:      verbosity level (up to 3)." << endl;
+      cout << "                    default: 0" << endl;
+      cout << "  --debug D:        debug level (up to 3)." << endl;
+      cout << "                    default: 0" << endl;
+      cout << "  --restart:        restart from previous run (which had produced resid.out and v.out)." << endl;
+      cout << "                    default: false" << endl;
+      if (rc == 0) exit(0);
+      return rc;
+    };
+
+    friend ostream & operator<< (ostream & ostr, options const & opt);
+
+    string fileA;
+    string fileB;
+    a_int nbEV;
+    a_int nbCV;
+    bool stdPb; // Standard or generalized (= not standard).
+    bool symPb;
+    bool cpxPb;
+    string mag; // Magnitude <=> "which" arpack parameter.
+    bool shiftReal, shiftImag;
+    double sigmaReal, sigmaImag; // Eigen value translation: look for lambda+sigma instead of lambda.
+    bool invert; // Eigen value invertion: look for 1./lambda instead of lambda.
+    double tol;
+    int maxIt;
+    bool schur;
+    string slv;
+    unique_ptr<double> slvItrTol;
+    unique_ptr<int> slvItrMaxIt;
+    unique_ptr<double> slvDrtPvtThd;
+    bool check;
+    int verbose;
+    int debug;
+    bool restart;
+};
+
+ostream & operator<< (ostream & ostr, options const & opt) {
+  ostr << "OPT: A " << opt.fileA << ", B " << opt.fileB;
+  ostr << ", nbEV " << opt.nbEV << ", nbCV " << opt.nbCV << ", stdPb " << (opt.stdPb ? "yes" : "no");
+  ostr << ", symPb " << (opt.symPb ? "yes" : "no") << ", mag " << opt.mag << endl;
+  ostr << "OPT: shiftReal " << (opt.shiftReal ? "yes" : "no") << ", sigmaReal " << opt.sigmaReal;
+  ostr << ", shiftImag " << (opt.shiftImag ? "yes" : "no") << ", sigmaImag " << opt.sigmaImag;
+  ostr << ", invert " << (opt.invert ? "yes" : "no") << ", tol " << opt.tol << ", maxIt " << opt.maxIt;
+  ostr << ", " << (opt.schur ? "Schur" : "Ritz") << " vectors" << endl;
+  ostr << "OPT: slv " << opt.slv;
+  if (opt.slvItrTol)    ostr << ", slvItrTol "    << *opt.slvItrTol;
+  if (opt.slvItrMaxIt)  ostr << ", slvItrMaxIt "  << *opt.slvItrMaxIt;
+  if (opt.slvDrtPvtThd) ostr << ", slvDrtPvtThd " << *opt.slvDrtPvtThd;
+  ostr << ", check " << (opt.check ? "yes" : "no") << ", verbose " << opt.verbose << ", debug " << opt.debug;
+  ostr << ", restart " << (opt.restart ? "yes" : "no") << endl;
+  return ostr;
+}
+
+void makeZero(        double  & zero) {zero = 0.;}
+
+void makeZero(complex<double> & zero) {zero = complex<double>(0., 0.);}
+
+template<typename RC, typename EM, typename EC>
+int readMatrixMarket(string const & fileName, EM & M, int const & verbose, string const & msg) {
+  ifstream inp(fileName);
+  if (!inp) {cerr << "Error: can not open " << fileName << endl; return 1;}
+
+  a_uint l = 0, n = 0, m = 0, nnz = 0;
+  vector<a_uint> i, j;
+  vector<RC> Mij;
+  do {
+    // Skip comments.
+
+    string inpLine; getline(inp, inpLine); l++;
+    while (isspace(*inpLine.begin())) inpLine.erase(inpLine.begin()); // Suppress leading white spaces.
+    if (inpLine.length() == 0) continue; // Empty line.
+    if (inpLine[0] == '%') continue; // Comments skipped, begin reading.
+
+    // Read matrix market file.
+
+    stringstream inpSS(inpLine);
+    if (n == 0 && m == 0) { // Header.
+      inpSS >> n >> m;
+      if (!inpSS) {cerr << "Error: bad header (n, m)" << endl; return 1;}
+      if (nnz == 0) {
+        inpSS >> nnz;
+        if (inpSS) { // OK, (optional) nnz has been provided.
+          i.reserve(nnz);
+          j.reserve(nnz);
+          Mij.reserve(nnz);
+        }
+      }
+    }
+    else { // Body.
+      a_uint k = 0, l = 0;
+      RC zero; makeZero(zero);
+      RC Mkl = zero;
+      inpSS >> k >> l >> Mkl;
+      if (!inpSS) {cerr << "Error: bad line (" << fileName << ", line " << l << ")" << endl; return 1;}
+      i.push_back(k);
+      j.push_back(l);
+      Mij.push_back(Mkl);
+    }
+  }
+  while (inp);
+
+  // Handle 1-based -> 0-based.
+
+  nnz = i.size(); // In case nnz was not provided.
+  if (*max_element(begin(i), end(i)) == n || *max_element(begin(j), end(j)) == m) {
+    for (size_t k = 0; k < nnz; k++) i[k] -= 1;
+    for (size_t k = 0; k < nnz; k++) j[k] -= 1;
+  }
+
+  // Create matrix from file.
+
+  M = EM(n, m); // Set matrice dimensions.
+  vector<EC> triplets;
+  triplets.reserve(nnz);
+  for (size_t k = 0; k < nnz; k++) triplets.emplace_back(i[k], j[k], Mij[k]);
+  M.setFromTriplets(triplets.begin(), triplets.end()); // Set all (i, j, Mij).
+
+  if (verbose == 3) {
+    cout << endl << msg << endl;
+    cout << endl << M << endl;
+  }
+
+  return 0;
+}
+
+class arpackEV { // Arpack eigen values / vectors.
+  public:
+    vector<complex<double>> val; // Eigen values.
+    vector<EigVecC> vec; // Eigen vectors.
+    int nbIt;
+    double rciTime;
+};
+
+void arpackAUPD(options const & opt,
+                a_int * ido, char const * bMat, a_int nbDim, char const * which, double * resid, double * v,
+                a_int ldv, a_int * iparam, a_int * ipntr, double * workd, double * workl, a_int lworkl, double * & rwork,
+                a_int * info) {
+  assert(rwork == NULL);
+  if (opt.symPb) {
+    dsaupd_c(ido, bMat, nbDim, which, opt.nbEV, opt.tol, resid, opt.nbCV, v, ldv, iparam, ipntr, workd, workl, lworkl, info);
+  }
+  else {
+    dnaupd_c(ido, bMat, nbDim, which, opt.nbEV, opt.tol, resid, opt.nbCV, v, ldv, iparam, ipntr, workd, workl, lworkl, info);
+  }
+}
+
+void arpackAUPD(options const & opt,
+                a_int * ido, char const * bMat, a_int nbDim, char const * which, complex<double> * resid, complex<double> * v,
+                a_int ldv, a_int * iparam, a_int * ipntr, complex<double> * workd, complex<double> * workl, a_int lworkl, double * & rwork,
+                a_int * info) {
+  if (!rwork) rwork = new double[opt.nbCV];
+  znaupd_c(ido, bMat, nbDim, which, opt.nbEV, opt.tol, reinterpret_cast<_Complex double*>(resid), opt.nbCV,
+           reinterpret_cast<_Complex double*>(v), ldv, iparam, ipntr, reinterpret_cast<_Complex double*>(workd),
+           reinterpret_cast<_Complex double*>(workl), lworkl, rwork, info);
+}
+
+int arpackEUPD(options const & opt, arpackEV & out,
+               bool rvec, char const * howmny, a_int const * select, double * z,
+               a_int ldz, char const * bMat, a_int nbDim, char const * which, double * resid, double * v,
+               a_int ldv, a_int * iparam, a_int * ipntr, double * workd, double * workl, a_int lworkl, double * rwork,
+               a_int & info) {
+  assert(rwork == NULL);
+  if (opt.symPb) {
+    double * d = new double[opt.nbEV]; for (a_int k = 0; k < opt.nbEV; k++) d[k] = 0.;
+
+    dseupd_c(rvec, howmny, select, d, z, ldz, opt.sigmaReal,
+             bMat, nbDim, which, opt.nbEV, opt.tol, resid, opt.nbCV, v, ldv, iparam, ipntr, workd, workl, lworkl, &info);
+    if (info == -14) cerr << "Error: dseupd - KO: dsaupd did not find any eigenvalues to sufficient accuracy" << endl;
+    if (info < 0 && info != -14 /*-14: don't break*/) {cerr << "Error: dseupd - KO with info " << info << endl; return 1;}
+
+    // Arpack compute the whole spectrum.
+
+    a_int nbConv = iparam[4];
+    out.val.reserve(nbConv);
+    for (a_int i = 0; d && i < nbConv; i++) {
+      complex<double> lambda(d[i], 0.);
+      out.val.push_back(lambda);
+      if (out.val.size() == (size_t) opt.nbEV) break; // If more converged than requested, likely not accurate (check KO).
+    }
+
+    out.vec.reserve(nbConv);
+    for (a_int i = 0; z && i < nbConv; i++) {
+      EigVecR V = EigMpVR(z + i*nbDim, nbDim);
+      out.vec.push_back(V.cast<complex<double>>());
+      if (out.vec.size() == (size_t) opt.nbEV) break; // If more converged than requested, likely not accurate (check KO).
+    }
+
+    if (d) {delete [] d; d = NULL;}
+  }
+  else {
+    double * dr = new double[opt.nbEV+1]; for (a_int k = 0; k < opt.nbEV+1; k++) dr[k] = 0.;
+    double * di = new double[opt.nbEV+1]; for (a_int k = 0; k < opt.nbEV+1; k++) di[k] = 0.;
+    double * workev = new double[3*opt.nbCV];
+
+    dneupd_c(rvec, howmny, select, dr, di, z, ldz, opt.sigmaReal, opt.sigmaImag, workev,
+             bMat, nbDim, which, opt.nbEV, opt.tol, resid, opt.nbCV, v, ldv, iparam, ipntr, workd, workl, lworkl, &info);
+    if (info == -14) cerr << "Error: dneupd - KO: [dz]naupd did not find any eigenvalues to sufficient accuracy" << endl;
+    if (info < 0 && info != -14 /*-14: don't break*/) {cerr << "Error: dneupd - KO with info " << info << endl; return 1;}
+
+    // Arpack compute only half of the spectrum.
+
+    a_int nbConv = iparam[4];
+    out.val.reserve(nbConv);
+    for (a_int i = 0; dr && di && i <= nbConv/2; i++) { // Scan first half of the spectrum.
+      // Get first half of the spectrum.
+
+      complex<double> lambda(dr[i], di[i]);
+      out.val.push_back(lambda);
+      if (out.val.size() == (size_t) opt.nbEV) break; // If more converged than requested, likely not accurate (check KO).
+
+      // Deduce second half of the spectrum.
+
+      out.val.push_back(complex<double>(lambda.real(), -1.*lambda.imag()));
+      if (out.val.size() == (size_t) opt.nbEV) break; // If more converged than requested, likely not accurate (check KO).
+    }
+
+    out.vec.reserve(nbConv);
+    for (a_int i = 0; z && i <= nbConv/2; i++) { // Scan half spectrum.
+      // Get first half of the spectrum.
+
+      EigVecR Vr = EigMpVR(z + (2*i+0)*nbDim, nbDim); // Real part.
+      EigVecR Vi = EigMpVR(z + (2*i+1)*nbDim, nbDim); // Imaginary part.
+      complex<double> imag(0., 1.);
+      EigVecC V = Vr.cast<complex<double>>() + imag * Vi.cast<complex<double>>();
+      out.vec.push_back(V);
+      if (out.vec.size() == (size_t) opt.nbEV) break; // If more converged than requested, likely not accurate (check KO).
+
+      // Deduce second half of the spectrum.
+
+      V = Vr.cast<complex<double>>() - imag * Vi.cast<complex<double>>();
+      out.vec.push_back(V);
+      if (out.vec.size() == (size_t) opt.nbEV) break; // If more converged than requested, likely not accurate (check KO).
+    }
+
+    if (workev) {delete [] workev; workev = NULL;}
+    if (dr)     {delete [] dr;     dr     = NULL;}
+    if (di)     {delete [] di;     di     = NULL;}
+  }
+
+  return 0;
+}
+
+int arpackEUPD(options const & opt, arpackEV & out,
+               bool rvec, char const * howmny, a_int const * select, complex<double> * z,
+               a_int ldz, char const * bMat, a_int nbDim, char const * which, complex<double> * resid, complex<double> * v,
+               a_int ldv, a_int * iparam, a_int * ipntr, complex<double> * workd, complex<double> * workl, a_int lworkl, double * rwork,
+               a_int & info) {
+  complex<double> * d = new complex<double>[opt.nbEV+1]; for (a_int k = 0; k < opt.nbEV+1; k++) d[k] = complex<double>(0., 0.);
+  complex<double> * workev = new complex<double>[2*opt.nbCV];
+
+  complex<double> sigma = complex<double>(opt.sigmaReal, opt.sigmaImag);
+  zneupd_c(rvec, howmny, select, reinterpret_cast<_Complex double*>(d), reinterpret_cast<_Complex double*>(z), ldz,
+           reinterpret_cast<_Complex double &>(sigma), reinterpret_cast<_Complex double*>(workev),
+           bMat, nbDim, which, opt.nbEV, opt.tol, reinterpret_cast<_Complex double*>(resid), opt.nbCV,
+           reinterpret_cast<_Complex double*>(v), ldv, iparam, ipntr,
+           reinterpret_cast<_Complex double*>(workd), reinterpret_cast<_Complex double*>(workl), lworkl, rwork, &info);
+  if (info == -14) cerr << "Error: zneupd - KO: dsaupd did not find any eigenvalues to sufficient accuracy" << endl;
+  if (info < 0 && info != -14 /*-14: don't break*/) {cerr << "Error: zneupd - KO with info " << info << endl; return 1;}
+
+  // Arpack compute the whole spectrum.
+
+  a_int nbConv = iparam[4];
+  out.val.reserve(nbConv);
+  for (a_int i = 0; d && i < nbConv; i++) {
+    complex<double> lambda = d[i];
+    out.val.push_back(lambda);
+    if (out.val.size() == (size_t) opt.nbEV) break; // If more converged than requested, likely not accurate (check KO).
+  }
+
+  out.vec.reserve(nbConv);
+  for (a_int i = 0; z && i < nbConv; i++) {
+    EigVecC V = EigMpVC(z + i*nbDim, nbDim);
+    out.vec.push_back(V);
+    if (out.vec.size() == (size_t) opt.nbEV) break; // If more converged than requested, likely not accurate (check KO).
+  }
+
+  if (workev) {delete [] workev; workev = NULL;}
+  if (d)      {delete [] d;           d = NULL;}
+
+  return 0;
+}
+
+template<typename SLV> int arpackMode(options const & opt, int const mode,
+                                      EigMatR const & A, EigMatR const & B, SLV & solver) {
+  int rc = 1;
+
+  if (mode == 1) {
+    rc = 0;
+  }
+  else if (mode == 2 || mode == 3) {
+    if (mode == 2) { // Regular mode.
+      solver.compute(B);
+    }
+    else { // Shift invert mode.
+      if (!opt.shiftImag) { // Real shift only.
+        double sigma = opt.sigmaReal;
+        auto S = A - sigma * B;
+        solver.compute(S);
+      }
+      else { // Complex (real/imaginary) shift.
+        complex<double> sigma(opt.sigmaReal, opt.sigmaImag);
+        auto S = A.cast<complex<double>>() - sigma * B.cast<complex<double>>();
+        solver.compute(S.real()); // Real part of shifted matrix.
+      }
+    }
+
+    if (solver.info() != Eigen::Success) {cerr << "Error: decomposition KO - check A and/or B are invertible" << endl; return 1;}
+    rc = 0;
+  }
+  else {cerr << "Error: arpack mode must be 1, 2 or 3 - KO" << endl; rc = 1;}
+
+  return rc;
+}
+
+template<typename SLV> int arpackMode(options const & opt, int const mode,
+                                      EigMatC const & A, EigMatC const & B, SLV & solver) {
+  int rc = 1;
+
+  if (mode == 1) {
+    rc = 0;
+  }
+  else if (mode == 2 || mode == 3) {
+    if (mode == 2) { // Regular mode.
+      solver.compute(B);
+    }
+    else { // Shift invert mode.
+      complex<double> sigma(opt.sigmaReal, opt.sigmaImag);
+      auto S = A - sigma * B;
+      solver.compute(S);
+    }
+
+    if (solver.info() != Eigen::Success) {cerr << "Error: decomposition KO - check A and/or B are invertible" << endl; return 1;}
+    rc = 0;
+  }
+  else {cerr << "Error: arpack mode must be 1, 2 or 3 - KO" << endl; rc = 1;}
+
+  return rc;
+}
+
+template<typename RC, typename EM, typename EV, typename SLV>
+int arpackSolve(options const & opt, int const & mode,
+                EM const & A, EM const & B, SLV & solver, arpackEV & out) {
+  // Arpack set up.
+
+  // Note: all in/out parameters (all but work*) passed to d[sn][ae]upd are set to 0. before use.
+  // d[sn][ae]upd uses dgetv0 to generate a random starting vector (when info is initialized to 0).
+  // dgetv0 rely on resid/v: resid/v should be initialized to 0.0 to avoid "bad" starting random vectors.
+
+  char const * which = opt.mag.c_str();
+  a_int ido = 0; // First call to arpack.
+  char const * iMat = "I";
+  char const * gMat = "G";
+  char const * bMat = (mode == 1) ? iMat : gMat;
+  a_int nbDim = A.rows();
+  RC zero; makeZero(zero);
+  RC * resid = new RC[nbDim]; for (a_int n = 0; n < nbDim; n++) resid[n] = zero; // Avoid "bad" starting vector.
+  if (opt.restart) {
+    ifstream rfs("resid.out");
+    if (rfs.is_open()) {
+      for (a_int n = 0; n < nbDim; n++) rfs >> resid[n];
+      if (opt.verbose >= 2) {
+        cout << endl;
+        cout << "resid:" << endl;
+        for (a_int n = 0; n < nbDim; n++) cout << resid[n] << endl;
+        cout << endl;
+      }
+    }
+  }
+  a_int ldv = nbDim;
+  RC * v = new RC[ldv*opt.nbCV]; for (a_int n = 0; n < ldv*opt.nbCV; n++) v[n] = zero; // Avoid "bad" starting vector.
+  if (opt.restart) {
+    ifstream vfs("v.out");
+    if (vfs.is_open()) {
+      a_int nbCV = 0; vfs >> nbCV; if (opt.nbCV < nbCV) nbCV = opt.nbCV;
+      for (a_int n = 0; n < ldv*nbCV; n++) vfs >> v[n];
+      if (opt.verbose >= 2) {
+        cout << endl;
+        cout << "v:" << endl;
+        for (a_int n = 0; n < ldv*nbCV; n++) cout << v[n] << endl;
+        cout << endl;
+      }
+    }
+  }
+  a_int iparam[11];
+  iparam[0] = 1; // Use exact shifts (=> we'll never have ido == 3).
+  iparam[2] = opt.maxIt; // Maximum number of iterations.
+  iparam[3] = 1; // Block size.
+  iparam[4] = 0; // Number of ev found by arpack.
+  iparam[6] = mode;
+  int rc = arpackMode<SLV>(opt, mode, A, B, solver);
+  if (rc != 0) {cerr << "Error: bad arpack mode" << endl; return rc;}
+  a_int ipntr[14];
+  RC * workd = new RC[3*nbDim];
+  a_int lworkl = opt.symPb ? opt.nbCV*opt.nbCV + 8*opt.nbCV : 3*opt.nbCV*opt.nbCV + 6*opt.nbCV;
+  lworkl++; // The documentation says "LWORKL must be at least ..."
+  RC * workl = new RC[lworkl];
+  a_int info = 0; // Use random initial residual vector.
+  if (opt.restart) info = 1;
+
+  // Arpack solve.
+
+  double * rwork = NULL;
+  do {
+    // Call arpack.
+
+    arpackAUPD(opt, &ido, bMat, nbDim, which, resid, v, ldv, iparam, ipntr, workd, workl, lworkl, rwork, &info);
+    if (info ==  1) cerr << "Error: [dz][sn]aupd - KO: maximum number of iterations taken. Increase --maxIt..." << endl;
+    if (info ==  3) cerr << "Error: [dz][sn]aupd - KO: no shifts could be applied. Increase --nbCV..." << endl;
+    if (info == -9) cerr << "Error: [dz][sn]aupd - KO: starting vector is zero. Retry: play with shift..." << endl;
+    if (info < 0) {cerr << "Error: [dz][sn]aupd - KO with info " << info << ", nbIt " << iparam[2] << endl; return 1;}
+
+    // Reverse Communication Interface: perform actions according to arpack.
+
+    auto start = chrono::high_resolution_clock::now();
+
+    a_int xIdx = ipntr[0] - 1; // 0-based (Fortran is 1-based).
+    a_int yIdx = ipntr[1] - 1; // 0-based (Fortran is 1-based).
+
+    EV X(workd + xIdx, nbDim); // Arpack provides X.
+    EV Y(workd + yIdx, nbDim); // Arpack provides Y.
+
+    if (ido == -1) {
+      if (iparam[6] == 1) {
+        Y = A * X;
+      }
+      else if (iparam[6] == 2) {
+        Y = A * X;
+        auto YY = Y;          // Use copy of Y (not Y) for solve (avoid potential memory overwrite as Y is both in/out).
+        Y = solver.solve(YY); // Y = B^-1 * A * X.
+        if(solver.info() != Eigen::Success) {
+          cerr << "Error: solve KO - play with solver parameters (tol, max it, ...), or, change --slv" << endl;
+          return 1;
+        }
+      }
+      else if (iparam[6] == 3) {
+        auto Z = B * X;      // Z = B * X.
+        Y = solver.solve(Z); // Y = (A - sigma * B)^-1 * B * X.
+        if(solver.info() != Eigen::Success) {
+          cerr << "Error: solve KO - play with solver parameters (tol, max it, ...), or, change --slv" << endl;
+          return 1;
+        }
+      }
+    }
+    else if (ido == 1) {
+      if (iparam[6] == 1) {
+        Y = A * X;
+      }
+      else if (iparam[6] == 2) {
+        Y = A * X;
+        if (opt.symPb) X = Y; // Remark 5 in dsaupd documentation.
+        auto YY = Y;          // Use copy of Y (not Y) for solve (avoid potential memory overwrite as Y is both in/out).
+        Y = solver.solve(YY); // Y = B^-1 * A * X.
+        if(solver.info() != Eigen::Success) {
+          cerr << "Error: solve KO - play with solver parameters (tol, max it, ...), or, change --slv" << endl;
+          return 1;
+        }
+      }
+      else if (iparam[6] == 3) {
+        a_int zIdx = ipntr[2] - 1; // 0-based (Fortran is 1-based).
+        EV Z(workd + zIdx, nbDim); // Arpack provides Z.
+        Y = solver.solve(Z);       // Y = (A - sigma * B)^-1 * B * X.
+        if(solver.info() != Eigen::Success) {
+          cerr << "Error: solve KO - play with solver parameters (tol, max it, ...), or, change --slv" << endl;
+          return 1;
+        }
+      }
+    }
+    else if (ido == 2) {
+      if      (iparam[6] == 1) Y =     X; // Y = I * X.
+      else if (iparam[6] == 2) Y = B * X; // Y = B * X.
+      else if (iparam[6] == 3) Y = B * X; // Y = B * X.
+    }
+    else if (ido != 99) {cerr << "Error: unexpected ido " << ido << " - KO" << endl; return 1;}
+
+    auto stop = chrono::high_resolution_clock::now();
+    out.rciTime += chrono::duration_cast<chrono::milliseconds>(stop - start).count()/1000.;
+
+  } while (ido != 99);
+
+  // Get arpack results (computed eigen values and vectors).
+
+  out.nbIt = iparam[2]; // Actual number of iterations.
+  bool rvec = true;
+  char const * howmnyA = "A"; // Ritz vectors.
+  char const * howmnyP = "P"; // Schur vectors.
+  char const * howmny = opt.schur ? howmnyP : howmnyA;
+  a_int * select = new a_int[opt.nbCV]; for (a_int n = 0; n < opt.nbCV; n++) select[n] = 1;
+  a_int const nbZ = nbDim*(opt.nbEV+1); // Caution: opt.nbEV+1 for dneupd.
+  RC * z = new RC[nbZ]; for (a_int n = 0; n < nbZ; n++) z[n] = zero;
+  a_int ldz = nbDim;
+  rc = arpackEUPD(opt, out, rvec, howmny, select, z, ldz, bMat, nbDim, which, resid, v, ldv, iparam, ipntr, workd, workl, lworkl, rwork, info);
+  if (rc != 0) {cerr << "Error: bad arpack eupd" << endl; return rc;}
+
+  ofstream rfs("resid.out"); for (a_int n = 0; n < nbDim; n++) rfs << resid[n] << endl;
+  ofstream vfs("v.out"); vfs << opt.nbCV << endl; for (a_int n = 0; n < ldv*opt.nbCV; n++) vfs << v[n] << endl;
+
+  // Clean.
+
+  if (rwork)  {delete [] rwork;  rwork  = NULL;}
+  if (z)      {delete [] z;      z      = NULL;}
+  if (select) {delete [] select; select = NULL;}
+  if (workl)  {delete [] workl;  workl  = NULL;}
+  if (workd)  {delete [] workd;  workd  = NULL;}
+  if (v)      {delete [] v;      v      = NULL;}
+  if (resid)  {delete [] resid;  resid  = NULL;}
+
+  return 0;
+}
+
+template<typename EM>
+int checkArpackEigVec(options const & opt, EM const & A, EM const & B, arpackEV const & out) {
+  // Check eigen vectors.
+
+  string rs = opt.schur ? "Schur" : "Ritz";
+
+  if (opt.check && out.vec.size() == 0) {
+    cerr << "Error: no " << rs << " value / vector found" << endl;
+    return 1;
+  }
+
+  for (size_t i = 0; i < out.vec.size(); i++) {
+    EigVecC V = out.vec[i];
+    complex<double> lambda = out.val[i];
+    if (opt.verbose >= 1) {
+      cout << endl;
+      cout << rs << " value " << setw(3) << i << ": " << lambda << endl;
+      if (opt.verbose >= 2) {
+        cout << endl;
+        cout << rs << " vector " << setw(3) << i << " (norm " << V.norm() << "): " << endl;
+        cout << endl << V << endl;
+      }
+    }
+
+    if (opt.check) {
+      EigVecC left = A.template cast<complex<double>>() * V;
+      EigVecC right = opt.stdPb ? V : B.template cast<complex<double>>() * V;
+      right *= lambda;
+      EigVecC diff = left - right;
+      if (diff.norm() > sqrt(opt.tol)) {
+        cerr << endl << "Error: bad vector " << setw(3) << i << " (norm " << V.norm() << "):" << endl;
+        cerr << endl << V << endl;
+        cerr << endl << "Error: left side (A*V - norm " << left.norm() << "):" << endl;
+        cerr << endl << left << endl;
+        cerr << endl << "Error: right side (lambda*" << (opt.stdPb ? "" : "B*") << "V - norm " << right.norm() << "):" << endl;
+        cerr << endl << right << endl;
+        cerr << endl << "Error: diff (norm " << diff.norm() << ", sqrt(tol) " << sqrt(opt.tol) << "):" << endl;
+        cerr << endl << diff << endl;
+        return 1;
+      }
+      else {
+        if (opt.verbose >= 1) {
+          cout << endl << rs << " value/vector " << setw(3) << i << ": check OK";
+          cout << ", diff (norm " << diff.norm() << ", sqrt(tol) " << sqrt(opt.tol) << ")" << endl;
+        }
+      }
+    }
+  }
+
+  return 0;
+}
+
+void makeSigma(options const & opt,         double  & sigma) {sigma = opt.sigmaReal;}
+
+void makeSigma(options const & opt, complex<double> & sigma) {sigma = complex<double>(opt.sigmaReal, opt.sigmaImag);}
+
+template<typename RC, typename EM, typename EV, typename SLV>
+int arpackSolve(options const & opt, EM & A, EM const & B,
+                SLV & solver, arpackEV & out) {
+  // If needed, transform the initial problem into a new one that arpack can handle.
+
+  auto eps = numeric_limits<double>::epsilon();
+  bool shiftReal = (opt.shiftReal && fabs(opt.sigmaReal) > eps) ? true : false;
+  bool shiftImag = (opt.shiftImag && fabs(opt.sigmaImag) > eps) ? true : false;
+
+  bool backTransform = false;
+  int mode = 0;
+  if (opt.stdPb) {
+    mode = 1;
+    if (shiftReal && !shiftImag) {
+      EM I(A.rows(), A.cols());
+      I.setIdentity();
+      RC sigma; makeSigma(opt, sigma);
+      A -= sigma*I;
+      backTransform = true;
+    }
+  }
+  else {
+    mode = 2;
+    if (shiftReal || shiftImag) mode = 3;
+  }
+
+  // Solve the problem.
+
+  if (opt.verbose >= 1) {
+    cout << endl;
+    cout << "ARP: mode " << mode;
+    cout << ", nbDim " << A.rows();
+    cout << ", backTransform " << (backTransform ? "yes" : "no") << endl;
+  }
+
+  int rc = arpackSolve<RC, EM, EV, SLV>(opt, mode, A, B, solver, out);
+  if (rc != 0) {cerr << "Error: arpack solve KO" << endl; return rc;}
+
+  if (opt.verbose >= 1) {
+    cout << endl;
+    cout << "ARP: nbEV found " << out.val.size();
+    cout << ", nbIt " << out.nbIt << endl;
+  }
+
+  // If needed, transform back the arpack problem into the initial problem.
+
+  if (backTransform) {
+    for (size_t i = 0; i < out.val.size(); i++) out.val[i] += opt.sigmaReal;
+    EM I(A.rows(), A.cols());
+    I.setIdentity();
+    RC sigma; makeSigma(opt, sigma);
+    A += sigma*I; // For later checks.
+  }
+
+  // Check.
+
+  return checkArpackEigVec<EM>(opt, A, B, out);
+}
+
+template<typename RC, typename EM, typename EC, typename EV, typename SLV>
+int arpackSolve(options & opt, SLV & solver) {
+  // Read A.
+
+  EM A;
+  int rc = readMatrixMarket<RC, EM, EC>(opt.fileA, A, opt.verbose, "A:");
+  if (rc != 0) {cerr << "Error: read A KO" << endl; return rc;}
+
+  // Read B.
+
+  EM B;
+  if (!opt.stdPb) {
+    rc = readMatrixMarket<RC, EM, EC>(opt.fileB, B, opt.verbose, "B:");
+    if (rc != 0) {cerr << "Error: read B KO" << endl; return rc;}
+  }
+
+  // Check A-B compatibility.
+
+  if (!opt.stdPb) {
+    if (A.rows() != B.rows()) {cerr << "Error: A.rows() != B.rows()" << endl; return rc;}
+    if (A.cols() != B.cols()) {cerr << "Error: A.cols() != B.cols()" << endl; return rc;}
+  }
+  if (opt.nbCV > A.cols()) opt.nbCV = A.cols(); // Cut-off.
+
+  // Arpack solve.
+
+  arpackEV out;
+  out.rciTime = 0.;
+  auto start = chrono::high_resolution_clock::now();
+  rc = arpackSolve<RC, EM, EV, SLV>(opt, A, B, solver, out);
+  if (rc != 0) {cerr << "Error: arpack solve KO" << endl; return rc;}
+  auto stop = chrono::high_resolution_clock::now();
+  double fullTime = chrono::duration_cast<chrono::milliseconds>(stop - start).count()/1000.;
+  cout << endl;
+  cout << "OUT: nb EV found " << out.val.size() << ", nb iterations " << out.nbIt << endl;
+  cout << "OUT: full time " << fullTime << " s, RCI time " << out.rciTime << " s" << endl;
+
+  return 0;
+}
+
+template<typename RC, typename EM, typename EC, typename EV,
+         typename SLVBCG, typename SLVCG, typename SLVSLU, typename SLVSQR>
+int arpackSolve(options & opt) {
+  // Solve with arpack.
+
+  int rc = 0;
+  if (opt.slv == "BiCG") {
+    SLVBCG solver;
+    if (opt.slvItrTol)   solver.setTolerance(*opt.slvItrTol);
+    if (opt.slvItrMaxIt) solver.setMaxIterations(*opt.slvItrMaxIt);
+    rc = arpackSolve<RC, EM, EC, EV, SLVBCG>(opt, solver);
+  }
+  else if (opt.slv == "CG") {
+    SLVCG solver;
+    if (opt.slvItrTol)   solver.setTolerance(*opt.slvItrTol);
+    if (opt.slvItrMaxIt) solver.setMaxIterations(*opt.slvItrMaxIt);
+    rc = arpackSolve<RC, EM, EC, EV, SLVCG>(opt, solver);
+  }
+  else if (opt.slv == "LU") {
+    SLVSLU solver;
+    if (opt.slvDrtPvtThd) solver.setPivotThreshold(*opt.slvDrtPvtThd);
+    rc = arpackSolve<RC, EM, EC, EV, SLVSLU>(opt, solver);
+  }
+  else if (opt.slv == "QR") {
+    SLVSQR solver;
+    if (opt.slvDrtPvtThd) solver.setPivotThreshold(*opt.slvDrtPvtThd);
+    rc = arpackSolve<RC, EM, EC, EV, SLVSQR>(opt, solver);
+  }
+  else {cerr << "Error: unknown solver - KO" << endl; return 1;}
+  if (rc != 0) {cerr << "Error: arpack solve KO" << endl; return rc;}
+
+  return 0;
+}
+
+int main(int argc, char ** argv) {
+  // Check for options.
+
+  options opt;
+  int rc = opt.readCmdLine(argc, argv);
+  if (rc != 0) {cerr << "Error: read cmd line KO" << endl; return rc;}
+  cout << opt; // Print options.
+
+  if (opt.cpxPb) rc = arpackSolve<complex<double>, EigMatC, EigCooC, EigMpVC, EigBiCGC, EigCGC, EigSLUC, EigSQRC>(opt);
+  else           rc = arpackSolve<        double , EigMatR, EigCooR, EigMpVR, EigBiCGR, EigCGR, EigSLUR, EigSQRR>(opt);
+  if (rc != 0) {cerr << "Error: arpack solve KO" << endl; return rc;}
+
+  return 0;
+}
+
+// Local Variables:
+// mode: c++
+// c-file-style:"stroustrup"
+// show-trailing-whitespace: t
+// End:
+/* vim: set sw=2 ts=2 et smartindent :*/

EXAMPLES/MATRIX_MARKET/arpackmm.sh

@@ -0,0 +1,69 @@
+#!/bin/bash -eu
+
+for eigPb in "--A As.mtx" "--nonSymPb --A An.mtx" "--nonSymPb --cpxPb --A Az.mtx --B Bz.mtx"
+do
+  for genPb in "" "--genPb"
+  do
+    for smallMag in "" "--mag SM --noCheck" # SM is known to be difficult to converge.
+    do
+      export shiftOpt=""
+      if [[ "$eigPb" == *nonSymPb* ]]; then
+        if [[ "$genPb" == *genPb* ]]; then
+          export shiftOpt="--shiftReal 2.5 --shiftImag 2.5 --tol 0.5" # Relax tolerance, tricky to converge.
+        else
+          export shiftOpt="--shiftReal 100.0 --shiftImag 100.0"
+        fi
+      else
+        if [[ "$genPb" == *genPb* ]]; then
+          export shiftOpt="--shiftReal 50.0"
+        else
+          export shiftOpt="--shiftReal 100.0"
+        fi
+      fi
+
+      for shiftRI in "" "$shiftOpt"
+      do
+        for invert in "" "--invert"
+        do
+          for tol in "" "--tol 1.e-5"
+          do
+            for slv in "" "--slv CG --slvItrTol 1.e-06 --slvItrMaxIt 100" "--slv LU --slvDrtPvtThd 1.e-06" "--slv QR"
+            do
+              for rs in "" "--schur"
+              do
+                export extraGenPb=""
+                if [[ "$genPb" == *genPb* ]]; then
+                  export extraGenPb="$shiftOpt" # Force shift if genPb.
+                fi
+
+                if [[ "$slv" == *CG* ]]; then
+                  if [[ "$eigPb" == *nonSymPb* ]]; then
+                    continue # Skip CG that could fail (CG is meant to deal with sym matrices).
+                  fi
+                fi
+
+                # Run arpackmm: use --nbCV 6 to ease convergence, and, --verbose 3 for debug.
+                export CMD="./arpackmm $eigPb $genPb $smallMag $shiftRI $invert $tol $slv $extraGenPb $rs --nbCV 6 --verbose 3"
+                echo "$CMD"
+                eval "$CMD"
+                echo ""
+                echo "========================================================================================"
+                echo ""
+
+                # Run arpackmm: re-run with restart.
+                export CMD="$CMD --restart"
+                echo "$CMD"
+                eval "$CMD"
+                echo ""
+                echo "========================================================================================"
+                echo ""
+              done
+            done
+          done
+        done
+      done
+    done
+  done
+done
+
+echo "OK"

EXAMPLES/Makefile.am

@@ -1,2 +1,6 @@
 SUBDIRS = BAND COMPLEX NONSYM SIMPLE SVD SYM
-EXTRA_DIST = README
+if ICBEXMM
+SUBDIRS += MATRIX_MARKET
+endif
+
+EXTRA_DIST = README README.CALLING-ARPACK-FROM-C-OR-CPP

ICB/Makefile.am

@@ -0,0 +1,14 @@
+AM_CPPFLAGS = -I$(top_builddir) -I$(top_srcdir)
+
+pkgincludedir = $(includedir)/arpack
+
+pkginclude_HEADERS = debug_c.h debug_c.hpp
+pkginclude_HEADERS += stat_c.h stat_c.hpp
+pkginclude_HEADERS += arpack.h arpack.hpp
+if MPI
+pkginclude_HEADERS += parpack.h parpack.hpp
+endif
+
+noinst_LTLIBRARIES = libdbgicb.la libstaicb.la
+libdbgicb_la_SOURCES = debug_icb.F90
+libstaicb_la_SOURCES = stat_icb.F90

ICB/arpack.h

@@ -0,0 +1,27 @@
+#ifndef __ARPACK_H__
+#define __ARPACK_H__
+
+#include "arpackdef.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+void cnaupd_c(a_int* ido, char const* bmat, a_int n, char const* which, a_int nev, float tol, float _Complex* resid, a_int ncv, float _Complex* v, a_int ldv, a_int* iparam, a_int* ipntr, float _Complex* workd, float _Complex* workl, a_int lworkl, float* rwork, a_int* info);
+void cneupd_c(bool rvec, char const* howmny, a_int const* select, float _Complex* d, float _Complex* z, a_int ldz, float _Complex sigma, float _Complex* workev, char const* bmat, a_int n, char const* which, a_int nev, float tol, float _Complex* resid, a_int ncv, float _Complex* v, a_int ldv, a_int* iparam, a_int* ipntr, float _Complex* workd, float _Complex* workl, a_int lworkl, float* rwork, a_int* info);
+void dnaupd_c(a_int* ido, char const* bmat, a_int n, char const* which, a_int nev, double tol, double* resid, a_int ncv, double* v, a_int ldv, a_int* iparam, a_int* ipntr, double* workd, double* workl, a_int lworkl, a_int* info);
+void dneupd_c(bool rvec, char const* howmny, a_int const* select, double* dr, double* di, double* z, a_int ldz, double sigmar, double sigmai, double * workev, char const* bmat, a_int n, char const* which, a_int nev, double tol, double* resid, a_int ncv, double* v, a_int ldv, a_int* iparam, a_int* ipntr, double* workd, double* workl, a_int lworkl, a_int* info);
+void dsaupd_c(a_int* ido, char const* bmat, a_int n, char const* which, a_int nev, double tol, double* resid, a_int ncv, double* v, a_int ldv, a_int* iparam, a_int* ipntr, double* workd, double* workl, a_int lworkl, a_int* info);
+void dseupd_c(bool rvec, char const* howmny, a_int const* select, double* d, double* z, a_int ldz, double sigma, char const* bmat, a_int n, char const* which, a_int nev, double tol, double* resid, a_int ncv, double* v, a_int ldv, a_int* iparam, a_int* ipntr, double* workd, double* workl, a_int lworkl, a_int* info);
+void snaupd_c(a_int* ido, char const* bmat, a_int n, char const* which, a_int nev, float tol, float* resid, a_int ncv, float* v, a_int ldv, a_int* iparam, a_int* ipntr, float* workd, float* workl, a_int lworkl, a_int* info);
+void sneupd_c(bool rvec, char const* howmny, a_int const* select, float* dr, float* di, float* z, a_int ldz, float sigmar, float sigmai, float * workev, char const* bmat, a_int n, char const* which, a_int nev, float tol, float* resid, a_int ncv, float* v, a_int ldv, a_int* iparam, a_int* ipntr, float* workd, float* workl, a_int lworkl, a_int* info);
+void ssaupd_c(a_int* ido, char const* bmat, a_int n, char const* which, a_int nev, float tol, float* resid, a_int ncv, float* v, a_int ldv, a_int* iparam, a_int* ipntr, float* workd, float* workl, a_int lworkl, a_int* info);
+void sseupd_c(bool rvec, char const* howmny, a_int const* select, float* d, float* z, a_int ldz, float sigma, char const* bmat, a_int n, char const* which, a_int nev, float tol, float* resid, a_int ncv, float* v, a_int ldv, a_int* iparam, a_int* ipntr, float* workd, float* workl, a_int lworkl, a_int* info);
+void znaupd_c(a_int* ido, char const* bmat, a_int n, char const* which, a_int nev, double tol, double _Complex* resid, a_int ncv, double _Complex* v, a_int ldv, a_int* iparam, a_int* ipntr, double _Complex* workd, double _Complex* workl, a_int lworkl, double* rwork, a_int* info);
+void zneupd_c(bool rvec, char const* howmny, a_int const* select, double _Complex* d, double _Complex* z, a_int ldz, double _Complex sigma, double _Complex* workev, char const* bmat, a_int n, char const* which, a_int nev, double tol, double _Complex* resid, a_int ncv, double _Complex* v, a_int ldv, a_int* iparam, a_int* ipntr, double _Complex* workd, double _Complex* workl, a_int lworkl, double* rwork, a_int* info);
+
+#ifdef  __cplusplus
+}
+#endif
+
+#endif

ICB/arpack.hpp

@@ -1,6 +1,8 @@
 #ifndef __ARPACK_HPP__
 #define __ARPACK_HPP__
 
+#include "arpackdef.h"
+
 #include <complex.h>
 #include <complex>
 
@@ -89,97 +91,99 @@ inline char const* convert_to_char(howmny const option) {
 }
 }  // namespace internal
 
-inline void saupd(int& ido, bmat const bmat_option, int n,
-                  which const ritz_option, int nev, float tol, float* resid,
-                  int ncv, float* v, int ldv, int* iparam, int* ipntr,
-                  float* workd, float* workl, int lworkl, int& info) {
+inline void saupd(a_int& ido, bmat const bmat_option, a_int n,
+                  which const ritz_option, a_int nev, float tol, float* resid,
+                  a_int ncv, float* v, a_int ldv, a_int* iparam, a_int* ipntr,
+                  float* workd, float* workl, a_int lworkl, a_int& info) {
   internal::ssaupd_c(&ido, internal::convert_to_char(bmat_option), n,
                      internal::convert_to_char(ritz_option), nev, tol, resid,
                      ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, &info);
 }
 
-inline void seupd(bool rvec, howmny const howmny_option, int* select, float* d,
-                  float* z, int ldz, float sigma, bmat const bmat_option, int n,
-                  which const ritz_option, int nev, float tol, float* resid,
-                  int ncv, float* v, int ldv, int* iparam, int* ipntr,
-                  float* workd, float* workl, int lworkl, int& info) {
+inline void seupd(bool rvec, howmny const howmny_option, a_int* select, float* d,
+                  float* z, a_int ldz, float sigma, bmat const bmat_option, a_int n,
+                  which const ritz_option, a_int nev, float tol, float* resid,
+                  a_int ncv, float* v, a_int ldv, a_int* iparam, a_int* ipntr,
+                  float* workd, float* workl, a_int lworkl, a_int& info) {
   internal::sseupd_c(rvec, internal::convert_to_char(howmny_option), select, d,
                      z, ldz, sigma, internal::convert_to_char(bmat_option), n,
                      internal::convert_to_char(ritz_option), nev, tol, resid,
                      ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, &info);
 }
 
-inline void saupd(int& ido, bmat const bmat_option, int n,
-                  which const ritz_option, int nev, double tol, double* resid,
-                  int ncv, double* v, int ldv, int* iparam, int* ipntr,
-                  double* workd, double* workl, int lworkl, int& info) {
+inline void saupd(a_int& ido, bmat const bmat_option, a_int n,
+                  which const ritz_option, a_int nev, double tol, double* resid,
+                  a_int ncv, double* v, a_int ldv, a_int* iparam, a_int* ipntr,
+                  double* workd, double* workl, a_int lworkl, a_int& info) {
   internal::dsaupd_c(&ido, internal::convert_to_char(bmat_option), n,
                      internal::convert_to_char(ritz_option), nev, tol, resid,
                      ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, &info);
 }
 
-inline void seupd(bool rvec, howmny const howmny_option, int* select, double* d,
-                  double* z, int ldz, double sigma, bmat const bmat_option,
-                  int n, which const ritz_option, int nev, double tol,
-                  double* resid, int ncv, double* v, int ldv, int* iparam,
-                  int* ipntr, double* workd, double* workl, int lworkl,
-                  int& info) {
+inline void seupd(bool rvec, howmny const howmny_option, a_int* select, double* d,
+                  double* z, a_int ldz, double sigma, bmat const bmat_option,
+                  a_int n, which const ritz_option, a_int nev, double tol,
+                  double* resid, a_int ncv, double* v, a_int ldv, a_int* iparam,
+                  a_int* ipntr, double* workd, double* workl, a_int lworkl,
+                  a_int& info) {
   internal::dseupd_c(rvec, internal::convert_to_char(howmny_option), select, d,
                      z, ldz, sigma, internal::convert_to_char(bmat_option), n,
                      internal::convert_to_char(ritz_option), nev, tol, resid,
                      ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, &info);
 }
 
-inline void naupd(int& ido, bmat const bmat_option, int n,
-                  which const ritz_option, int nev, float tol, float* resid,
-                  int ncv, float* v, int ldv, int* iparam, int* ipntr,
-                  float* workd, float* workl, int lworkl, int& info) {
+inline void naupd(a_int& ido, bmat const bmat_option, a_int n,
+                  which const ritz_option, a_int nev, float tol, float* resid,
+                  a_int ncv, float* v, a_int ldv, a_int* iparam, a_int* ipntr,
+                  float* workd, float* workl, a_int lworkl, a_int& info) {
   internal::snaupd_c(&ido, internal::convert_to_char(bmat_option), n,
                      internal::convert_to_char(ritz_option), nev, tol, resid,
                      ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, &info);
 }
 
-inline void neupd(bool rvec, howmny const howmny_option, int* select, float* dr,
-                  float* di, float* z, int ldz, float sigmar, float sigmai,
-                  bmat const bmat_option, int n, which const ritz_option,
-                  int nev, float tol, float* resid, int ncv, float* v, int ldv,
-                  int* iparam, int* ipntr, float* workd, float* workl,
-                  int lworkl, int& info) {
+inline void neupd(bool rvec, howmny const howmny_option, a_int* select, float* dr,
+                  float* di, float* z, a_int ldz,
+                  float sigmar, float sigmai, float * workev,
+                  bmat const bmat_option, a_int n, which const ritz_option,
+                  a_int nev, float tol, float* resid, a_int ncv, float* v, a_int ldv,
+                  a_int* iparam, a_int* ipntr, float* workd, float* workl,
+                  a_int lworkl, a_int& info) {
   internal::sneupd_c(rvec, internal::convert_to_char(howmny_option), select, dr,
-                     di, z, ldz, sigmar, sigmai,
+                     di, z, ldz, sigmar, sigmai, workev,
                      internal::convert_to_char(bmat_option), n,
                      internal::convert_to_char(ritz_option), nev, tol, resid,
                      ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, &info);
 }
 
-inline void naupd(int& ido, bmat const bmat_option, int n,
-                  which const ritz_option, int nev, double tol, double* resid,
-                  int ncv, double* v, int ldv, int* iparam, int* ipntr,
-                  double* workd, double* workl, int lworkl, int& info) {
+inline void naupd(a_int& ido, bmat const bmat_option, a_int n,
+                  which const ritz_option, a_int nev, double tol, double* resid,
+                  a_int ncv, double* v, a_int ldv, a_int* iparam, a_int* ipntr,
+                  double* workd, double* workl, a_int lworkl, a_int& info) {
   internal::dnaupd_c(&ido, internal::convert_to_char(bmat_option), n,
                      internal::convert_to_char(ritz_option), nev, tol, resid,
                      ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, &info);
 }
 
-inline void neupd(bool rvec, howmny const howmny_option, int* select,
-                  double* dr, double* di, double* z, int ldz, double sigmar,
-                  double sigmai, bmat const bmat_option, int n,
-                  which const ritz_option, int nev, double tol, double* resid,
-                  int ncv, double* v, int ldv, int* iparam, int* ipntr,
-                  double* workd, double* workl, int lworkl, int& info) {
+inline void neupd(bool rvec, howmny const howmny_option, a_int* select,
+                  double* dr, double* di, double* z, a_int ldz,
+                  double sigmar, double sigmai, double * workev,
+                  bmat const bmat_option, a_int n,
+                  which const ritz_option, a_int nev, double tol, double* resid,
+                  a_int ncv, double* v, a_int ldv, a_int* iparam, a_int* ipntr,
+                  double* workd, double* workl, a_int lworkl, a_int& info) {
   internal::dneupd_c(rvec, internal::convert_to_char(howmny_option), select, dr,
-                     di, z, ldz, sigmar, sigmai,
+                     di, z, ldz, sigmar, sigmai, workev,
                      internal::convert_to_char(bmat_option), n,
                      internal::convert_to_char(ritz_option), nev, tol, resid,
                      ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, &info);
 }
 
-inline void naupd(int& ido, bmat const bmat_option, int n,
-                  which const ritz_option, int nev, float tol,
-                  std::complex<float>* resid, int ncv, std::complex<float>* v,
-                  int ldv, int* iparam, int* ipntr, std::complex<float>* workd,
-                  std::complex<float>* workl, int lworkl,
-                  float* rwork, int& info) {
+inline void naupd(a_int& ido, bmat const bmat_option, a_int n,
+                  which const ritz_option, a_int nev, float tol,
+                  std::complex<float>* resid, a_int ncv, std::complex<float>* v,
+                  a_int ldv, a_int* iparam, a_int* ipntr, std::complex<float>* workd,
+                  std::complex<float>* workl, a_int lworkl,
+                  float* rwork, a_int& info) {
   internal::cnaupd_c(&ido, internal::convert_to_char(bmat_option), n,
                      internal::convert_to_char(ritz_option), nev, tol,
                      reinterpret_cast<_Complex float*>(resid), ncv,
@@ -189,14 +193,14 @@ inline void naupd(int& ido, bmat const bmat_option, int n,
                      rwork, &info);
 }
 
-inline void neupd(bool rvec, howmny const howmny_option, int* select,
-                  std::complex<float>* d, std::complex<float>* z, int ldz,
+inline void neupd(bool rvec, howmny const howmny_option, a_int* select,
+                  std::complex<float>* d, std::complex<float>* z, a_int ldz,
                   std::complex<float> sigma, std::complex<float>* workev,
-                  bmat const bmat_option, int n, which const ritz_option,
-                  int nev, float tol, std::complex<float>* resid, int ncv,
-                  std::complex<float>* v, int ldv, int* iparam, int* ipntr,
+                  bmat const bmat_option, a_int n, which const ritz_option,
+                  a_int nev, float tol, std::complex<float>* resid, a_int ncv,
+                  std::complex<float>* v, a_int ldv, a_int* iparam, a_int* ipntr,
                   std::complex<float>* workd, std::complex<float>* workl,
-                  int lworkl, float* rwork, int& info) {
+                  a_int lworkl, float* rwork, a_int& info) {
   internal::cneupd_c(rvec, internal::convert_to_char(howmny_option), select,
                      reinterpret_cast<_Complex float*>(d),
                      reinterpret_cast<_Complex float*>(z), ldz,
@@ -211,12 +215,12 @@ inline void neupd(bool rvec, howmny const howmny_option, int* select,
                      rwork, &info);
 }
 
-inline void naupd(int& ido, bmat const bmat_option, int n,
-                  which const ritz_option, int nev, double tol,
-                  std::complex<double>* resid, int ncv, std::complex<double>* v,
-                  int ldv, int* iparam, int* ipntr, std::complex<double>* workd,
-                  std::complex<double>* workl, int lworkl,
-                  double* rwork, int& info) {
+inline void naupd(a_int& ido, bmat const bmat_option, a_int n,
+                  which const ritz_option, a_int nev, double tol,
+                  std::complex<double>* resid, a_int ncv, std::complex<double>* v,
+                  a_int ldv, a_int* iparam, a_int* ipntr, std::complex<double>* workd,
+                  std::complex<double>* workl, a_int lworkl,
+                  double* rwork, a_int& info) {
   internal::znaupd_c(&ido, internal::convert_to_char(bmat_option), n,
                      internal::convert_to_char(ritz_option), nev, tol,
                      reinterpret_cast<_Complex double*>(resid), ncv,
@@ -226,14 +230,14 @@ inline void naupd(int& ido, bmat const bmat_option, int n,
                      rwork, &info);
 }
 
-inline void neupd(bool rvec, howmny const howmny_option, int* select,
-                  std::complex<double>* d, std::complex<double>* z, int ldz,
+inline void neupd(bool rvec, howmny const howmny_option, a_int* select,
+                  std::complex<double>* d, std::complex<double>* z, a_int ldz,
                   std::complex<double> sigma, std::complex<double>* workev,
-                  bmat const bmat_option, int n, which const ritz_option,
-                  int nev, double tol, std::complex<double>* resid, int ncv,
-                  std::complex<double>* v, int ldv, int* iparam, int* ipntr,
+                  bmat const bmat_option, a_int n, which const ritz_option,
+                  a_int nev, double tol, std::complex<double>* resid, a_int ncv,
+                  std::complex<double>* v, a_int ldv, a_int* iparam, a_int* ipntr,
                   std::complex<double>* workd, std::complex<double>* workl,
-                  int lworkl, double* rwork, int& info) {
+                  a_int lworkl, double* rwork, a_int& info) {
   internal::zneupd_c(rvec, internal::convert_to_char(howmny_option), select,
                      reinterpret_cast<_Complex double*>(d),
                      reinterpret_cast<_Complex double*>(z), ldz,

ICB/debug_icb.F90

@@ -7,6 +7,7 @@ subroutine debug_c(logfil_c, ndigit_c, mgetv0_c,
                    bind(c, name="debug_c")
   use :: iso_c_binding
   implicit none
+#include "arpackdef.h"
   integer(kind=c_int), value, intent(in) :: logfil_c, ndigit_c, mgetv0_c
   integer(kind=c_int), value, intent(in) :: msaupd_c, msaup2_c, msaitr_c, mseigt_c, msapps_c, msgets_c, mseupd_c
   integer(kind=c_int), value, intent(in) :: mnaupd_c, mnaup2_c, mnaitr_c, mneigh_c, mnapps_c, mngets_c, mneupd_c

ICB/parpack.h

@@ -0,0 +1,34 @@
+#ifndef __PARPACK_H__
+#define __PARPACK_H__
+
+#include "arpackdef.h"
+
+/*
+ * IMPORTANT: MPI communicators MUST be passed from C to Fortran using MPI_Comm_c2f.
+ *            MPI_Fint MCW = MPI_Comm_c2f(MPI_COMM_WORLD);
+ */
+
+#include "mpi.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+void pcnaupd_c(MPI_Fint comm, a_int* ido, char const* bmat, a_int n, char const* which, a_int nev, float tol, float _Complex* resid, a_int ncv, float _Complex* v, a_int ldv, a_int* iparam, a_int* ipntr, float _Complex* workd, float _Complex* workl, a_int lworkl, float _Complex* rwork, a_int* info);
+void pcneupd_c(MPI_Fint comm, bool rvec, char const* howmny, a_int const* select, float _Complex* d, float _Complex* z, a_int ldz, float _Complex sigma, float _Complex* workev, char const* bmat, a_int n, char const* which, a_int nev, float tol, float _Complex* resid, a_int ncv, float _Complex* v, a_int ldv, a_int* iparam, a_int* ipntr, float _Complex* workd, float _Complex* workl, a_int lworkl, float _Complex* rwork, a_int* info);
+void pdnaupd_c(MPI_Fint comm, a_int* ido, char const* bmat, a_int n, char const* which, a_int nev, double tol, double* resid, a_int ncv, double* v, a_int ldv, a_int* iparam, a_int* ipntr, double* workd, double* workl, a_int lworkl, a_int* info);
+void pdneupd_c(MPI_Fint comm, bool rvec, char const* howmny, a_int const* select, double* dr, double* di, double* z, a_int ldz, double sigmar, double sigmai, double * workev, char const* bmat, a_int n, char const* which, a_int nev, double tol, double* resid, a_int ncv, double* v, a_int ldv, a_int* iparam, a_int* ipntr, double* workd, double* workl, a_int lworkl, a_int* info);
+void pdsaupd_c(MPI_Fint comm, a_int* ido, char const* bmat, a_int n, char const* which, a_int nev, double tol, double* resid, a_int ncv, double* v, a_int ldv, a_int* iparam, a_int* ipntr, double* workd, double* workl, a_int lworkl, a_int* info);
+void pdseupd_c(MPI_Fint comm, bool rvec, char const* howmny, a_int const* select, double* d, double* z, a_int ldz, double sigma, char const* bmat, a_int n, char const* which, a_int nev, double tol, double* resid, a_int ncv, double* v, a_int ldv, a_int* iparam, a_int* ipntr, double* workd, double* workl, a_int lworkl, a_int* info);
+void psnaupd_c(MPI_Fint comm, a_int* ido, char const* bmat, a_int n, char const* which, a_int nev, float tol, float* resid, a_int ncv, float* v, a_int ldv, a_int* iparam, a_int* ipntr, float* workd, float* workl, a_int lworkl, a_int* info);
+void psneupd_c(MPI_Fint comm, bool rvec, char const* howmny, a_int const* select, float* dr, float* di, float* z, a_int ldz, float sigmar, float sigmai, float * workev, char const* bmat, a_int n, char const* which, a_int nev, float tol, float* resid, a_int ncv, float* v, a_int ldv, a_int* iparam, a_int* ipntr, float* workd, float* workl, a_int lworkl, a_int* info);
+void pssaupd_c(MPI_Fint comm, a_int* ido, char const* bmat, a_int n, char const* which, a_int nev, float tol, float* resid, a_int ncv, float* v, a_int ldv, a_int* iparam, a_int* ipntr, float* workd, float* workl, a_int lworkl, a_int* info);
+void psseupd_c(MPI_Fint comm, bool rvec, char const* howmny, a_int const* select, float* d, float* z, a_int ldz, float sigma, char const* bmat, a_int n, char const* which, a_int nev, float tol, float* resid, a_int ncv, float* v, a_int ldv, a_int* iparam, a_int* ipntr, float* workd, float* workl, a_int lworkl, a_int* info);
+void pznaupd_c(MPI_Fint comm, a_int* ido, char const* bmat, a_int n, char const* which, a_int nev, double tol, double _Complex* resid, a_int ncv, double _Complex* v, a_int ldv, a_int* iparam, a_int* ipntr, double _Complex* workd, double _Complex* workl, a_int lworkl, double _Complex* rwork, a_int* info);
+void pzneupd_c(MPI_Fint comm, bool rvec, char const* howmny, a_int const* select, double _Complex* d, double _Complex* z, a_int ldz, double _Complex sigma, double _Complex* workev, char const* bmat, a_int n, char const* which, a_int nev, double tol, double _Complex* resid, a_int ncv, double _Complex* v, a_int ldv, a_int* iparam, a_int* ipntr, double _Complex* workd, double _Complex* workl, a_int lworkl, double _Complex* rwork, a_int* info);
+
+#ifdef  __cplusplus
+}
+#endif
+
+#endif

ICB/parpack.hpp

@@ -1,7 +1,8 @@
-
 #ifndef __PARPACK_HPP__
 #define __PARPACK_HPP__
 
+#include "arpackdef.h"
+
 #include "arpack.hpp"
 
 #include <mpi.h>
@@ -15,21 +16,21 @@ namespace internal {
 #include "parpack.h"
 }  // namespace internal
 
-inline void saupd(MPI_Fint comm, int& ido, bmat const bmat_option, int n,
-                  which const which_option, int nev, float tol, float* resid,
-                  int ncv, float* v, int ldv, int* iparam, int* ipntr,
-                  float* workd, float* workl, int lworkl, int& info) {
+inline void saupd(MPI_Fint comm, a_int& ido, bmat const bmat_option, a_int n,
+                  which const which_option, a_int nev, float tol, float* resid,
+                  a_int ncv, float* v, a_int ldv, a_int* iparam, a_int* ipntr,
+                  float* workd, float* workl, a_int lworkl, a_int& info) {
   internal::pssaupd_c(comm, &ido, internal::convert_to_char(bmat_option), n,
                       internal::convert_to_char(which_option), nev, tol, resid,
                       ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, &info);
 }
 
 inline void seupd(MPI_Fint comm, bool rvec, howmny const howmny_option,
-                  int* select, float* d, float* z, int ldz, float sigma,
-                  bmat const bmat_option, int n, which const which_option,
-                  int nev, float tol, float* resid, int ncv, float* v, int ldv,
-                  int* iparam, int* ipntr, float* workd, float* workl,
-                  int lworkl, int& info) {
+                  a_int* select, float* d, float* z, a_int ldz, float sigma,
+                  bmat const bmat_option, a_int n, which const which_option,
+                  a_int nev, float tol, float* resid, a_int ncv, float* v, a_int ldv,
+                  a_int* iparam, a_int* ipntr, float* workd, float* workl,
+                  a_int lworkl, a_int& info) {
   internal::psseupd_c(comm, rvec, internal::convert_to_char(howmny_option),
                       select, d, z, ldz, sigma,
                       internal::convert_to_char(bmat_option), n,
@@ -37,21 +38,21 @@ inline void seupd(MPI_Fint comm, bool rvec, howmny const howmny_option,
                       ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, &info);
 }
 
-inline void saupd(MPI_Fint comm, int& ido, bmat const bmat_option, int n,
-                  which const which_option, int nev, double tol, double* resid,
-                  int ncv, double* v, int ldv, int* iparam, int* ipntr,
-                  double* workd, double* workl, int lworkl, int& info) {
+inline void saupd(MPI_Fint comm, a_int& ido, bmat const bmat_option, a_int n,
+                  which const which_option, a_int nev, double tol, double* resid,
+                  a_int ncv, double* v, a_int ldv, a_int* iparam, a_int* ipntr,
+                  double* workd, double* workl, a_int lworkl, a_int& info) {
   internal::pdsaupd_c(comm, &ido, internal::convert_to_char(bmat_option), n,
                       internal::convert_to_char(which_option), nev, tol, resid,
                       ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, &info);
 }
 
 inline void seupd(MPI_Fint comm, bool rvec, howmny const howmny_option,
-                  int* select, double* d, double* z, int ldz, double sigma,
-                  bmat const bmat_option, int n, which const which_option,
-                  int nev, double tol, double* resid, int ncv, double* v,
-                  int ldv, int* iparam, int* ipntr, double* workd,
-                  double* workl, int lworkl, int& info) {
+                  a_int* select, double* d, double* z, a_int ldz, double sigma,
+                  bmat const bmat_option, a_int n, which const which_option,
+                  a_int nev, double tol, double* resid, a_int ncv, double* v,
+                  a_int ldv, a_int* iparam, a_int* ipntr, double* workd,
+                  double* workl, a_int lworkl, a_int& info) {
   internal::pdseupd_c(comm, rvec, internal::convert_to_char(howmny_option),
                       select, d, z, ldz, sigma,
                       internal::convert_to_char(bmat_option), n,
@@ -59,56 +60,58 @@ inline void seupd(MPI_Fint comm, bool rvec, howmny const howmny_option,
                       ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, &info);
 }
 
-inline void naupd(MPI_Fint comm, int& ido, bmat const bmat_option, int n,
-                  which const which_option, int nev, float tol, float* resid,
-                  int ncv, float* v, int ldv, int* iparam, int* ipntr,
-                  float* workd, float* workl, int lworkl, int& info) {
+inline void naupd(MPI_Fint comm, a_int& ido, bmat const bmat_option, a_int n,
+                  which const which_option, a_int nev, float tol, float* resid,
+                  a_int ncv, float* v, a_int ldv, a_int* iparam, a_int* ipntr,
+                  float* workd, float* workl, a_int lworkl, a_int& info) {
   internal::psnaupd_c(comm, &ido, internal::convert_to_char(bmat_option), n,
                       internal::convert_to_char(which_option), nev, tol, resid,
                       ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, &info);
 }
 
 inline void neupd(MPI_Fint comm, bool rvec, howmny const howmny_option,
-                  int* select, float* dr, float* di, float* z, int ldz,
-                  float sigmar, float sigmai, bmat const bmat_option, int n,
-                  which const which_option, int nev, float tol, float* resid,
-                  int ncv, float* v, int ldv, int* iparam, int* ipntr,
-                  float* workd, float* workl, int lworkl, int& info) {
+                  a_int* select, float* dr, float* di, float* z, a_int ldz,
+                  float sigmar, float sigmai, float * workev,
+                  bmat const bmat_option, a_int n,
+                  which const which_option, a_int nev, float tol, float* resid,
+                  a_int ncv, float* v, a_int ldv, a_int* iparam, a_int* ipntr,
+                  float* workd, float* workl, a_int lworkl, a_int& info) {
   internal::psneupd_c(comm, rvec, internal::convert_to_char(howmny_option),
-                      select, dr, di, z, ldz, sigmar, sigmai,
+                      select, dr, di, z, ldz, sigmar, sigmai, workev,
                       internal::convert_to_char(bmat_option), n,
                       internal::convert_to_char(which_option), nev, tol, resid,
                       ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, &info);
 }
 
-inline void naupd(MPI_Fint comm, int& ido, bmat const bmat_option, int n,
-                  which const which_option, int nev, double tol, double* resid,
-                  int ncv, double* v, int ldv, int* iparam, int* ipntr,
-                  double* workd, double* workl, int lworkl, int& info) {
+inline void naupd(MPI_Fint comm, a_int& ido, bmat const bmat_option, a_int n,
+                  which const which_option, a_int nev, double tol, double* resid,
+                  a_int ncv, double* v, a_int ldv, a_int* iparam, a_int* ipntr,
+                  double* workd, double* workl, a_int lworkl, a_int& info) {
   internal::pdnaupd_c(comm, &ido, internal::convert_to_char(bmat_option), n,
                       internal::convert_to_char(which_option), nev, tol, resid,
                       ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, &info);
 }
 
 inline void neupd(MPI_Fint comm, bool rvec, howmny const howmny_option,
-                  int* select, double* dr, double* di, double* z, int ldz,
-                  double sigmar, double sigmai, bmat const bmat_option, int n,
-                  which const which_option, int nev, double tol, double* resid,
-                  int ncv, double* v, int ldv, int* iparam, int* ipntr,
-                  double* workd, double* workl, int lworkl, int& info) {
+                  a_int* select, double* dr, double* di, double* z, a_int ldz,
+                  double sigmar, double sigmai, double * workev,
+                  bmat const bmat_option, a_int n,
+                  which const which_option, a_int nev, double tol, double* resid,
+                  a_int ncv, double* v, a_int ldv, a_int* iparam, a_int* ipntr,
+                  double* workd, double* workl, a_int lworkl, a_int& info) {
   internal::pdneupd_c(comm, rvec, internal::convert_to_char(howmny_option),
-                      select, dr, di, z, ldz, sigmar, sigmai,
+                      select, dr, di, z, ldz, sigmar, sigmai, workev,
                       internal::convert_to_char(bmat_option), n,
                       internal::convert_to_char(which_option), nev, tol, resid,
                       ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, &info);
 }
 
-inline void naupd(MPI_Fint comm, int& ido, bmat const bmat_option, int n,
-                  which const which_option, int nev, float tol,
-                  std::complex<float>* resid, int ncv, std::complex<float>* v,
-                  int ldv, int* iparam, int* ipntr, std::complex<float>* workd,
-                  std::complex<float>* workl, int lworkl,
-                  std::complex<float>* rwork, int& info) {
+inline void naupd(MPI_Fint comm, a_int& ido, bmat const bmat_option, a_int n,
+                  which const which_option, a_int nev, float tol,
+                  std::complex<float>* resid, a_int ncv, std::complex<float>* v,
+                  a_int ldv, a_int* iparam, a_int* ipntr, std::complex<float>* workd,
+                  std::complex<float>* workl, a_int lworkl,
+                  std::complex<float>* rwork, a_int& info) {
   internal::pcnaupd_c(comm, &ido, internal::convert_to_char(bmat_option), n,
                       internal::convert_to_char(which_option), nev, tol,
                       reinterpret_cast<_Complex float*>(resid), ncv,
@@ -119,14 +122,14 @@ inline void naupd(MPI_Fint comm, int& ido, bmat const bmat_option, int n,
 }
 
 inline void neupd(MPI_Fint comm, bool rvec, howmny const howmny_option,
-                  int* select, std::complex<float>* d, std::complex<float>* z,
-                  int ldz, std::complex<float> sigma,
-                  std::complex<float>* workev, bmat const bmat_option, int n,
-                  which const which_option, int nev, float tol,
-                  std::complex<float>* resid, int ncv, std::complex<float>* v,
-                  int ldv, int* iparam, int* ipntr, std::complex<float>* workd,
-                  std::complex<float>* workl, int lworkl,
-                  std::complex<float>* rwork, int& info)
+                  a_int* select, std::complex<float>* d, std::complex<float>* z,
+                  a_int ldz, std::complex<float> sigma,
+                  std::complex<float>* workev, bmat const bmat_option, a_int n,
+                  which const which_option, a_int nev, float tol,
+                  std::complex<float>* resid, a_int ncv, std::complex<float>* v,
+                  a_int ldv, a_int* iparam, a_int* ipntr, std::complex<float>* workd,
+                  std::complex<float>* workl, a_int lworkl,
+                  std::complex<float>* rwork, a_int& info)
 
 {
   internal::pcneupd_c(comm, rvec, internal::convert_to_char(howmny_option),
@@ -143,12 +146,12 @@ inline void neupd(MPI_Fint comm, bool rvec, howmny const howmny_option,
                       reinterpret_cast<_Complex float*>(rwork), &info);
 }
 
-inline void naupd(MPI_Fint comm, int& ido, bmat const bmat_option, int n,
-                  which const which_option, int nev, double tol,
-                  std::complex<double>* resid, int ncv, std::complex<double>* v,
-                  int ldv, int* iparam, int* ipntr, std::complex<double>* workd,
-                  std::complex<double>* workl, int lworkl,
-                  std::complex<double>* rwork, int& info) {
+inline void naupd(MPI_Fint comm, a_int& ido, bmat const bmat_option, a_int n,
+                  which const which_option, a_int nev, double tol,
+                  std::complex<double>* resid, a_int ncv, std::complex<double>* v,
+                  a_int ldv, a_int* iparam, a_int* ipntr, std::complex<double>* workd,
+                  std::complex<double>* workl, a_int lworkl,
+                  std::complex<double>* rwork, a_int& info) {
   internal::pznaupd_c(comm, &ido, internal::convert_to_char(bmat_option), n,
                       internal::convert_to_char(which_option), nev, tol,
                       reinterpret_cast<_Complex double*>(resid), ncv,
@@ -159,14 +162,14 @@ inline void naupd(MPI_Fint comm, int& ido, bmat const bmat_option, int n,
 }
 
 inline void neupd(MPI_Fint comm, bool rvec, howmny const howmny_option,
-                  int* select, std::complex<double>* d, std::complex<double>* z,
-                  int ldz, std::complex<double> sigma,
-                  std::complex<double>* workev, bmat const bmat_option, int n,
-                  which const which_option, int nev, double tol,
-                  std::complex<double>* resid, int ncv, std::complex<double>* v,
-                  int ldv, int* iparam, int* ipntr, std::complex<double>* workd,
-                  std::complex<double>* workl, int lworkl,
-                  std::complex<double>* rwork, int& info) {
+                  a_int* select, std::complex<double>* d, std::complex<double>* z,
+                  a_int ldz, std::complex<double> sigma,
+                  std::complex<double>* workev, bmat const bmat_option, a_int n,
+                  which const which_option, a_int nev, double tol,
+                  std::complex<double>* resid, a_int ncv, std::complex<double>* v,
+                  a_int ldv, a_int* iparam, a_int* ipntr, std::complex<double>* workd,
+                  std::complex<double>* workl, a_int lworkl,
+                  std::complex<double>* rwork, a_int& info) {
   internal::pzneupd_c(comm, rvec, internal::convert_to_char(howmny_option),
                       select, reinterpret_cast<_Complex double*>(d),
                       reinterpret_cast<_Complex double*>(z), ldz,

ICB/stat_icb.F90

@@ -26,6 +26,7 @@ subroutine stat_c(  nopx_c,    nbx_c, nrorth_c, nitref_c, nrstrt_c,
                   bind(c, name="stat_c")
   use :: iso_c_binding
   implicit none
+#include "arpackdef.h"
   integer(kind=c_int), intent(out) ::   nopx_c,    nbx_c, nrorth_c, nitref_c, nrstrt_c
   real(kind=c_float),  intent(out) :: tsaupd_c, tsaup2_c, tsaitr_c, tseigt_c, tsgets_c, tsapps_c, tsconv_c,&
                                       tnaupd_c, tnaup2_c, tnaitr_c, tneigh_c, tngets_c, tnapps_c, tnconv_c,&

Makefile.am

@@ -1,36 +1,28 @@
 ACLOCAL_AMFLAGS = -I m4
 
-arpack_includedir = $(includedir)/arpack
+pkgincludedir = $(includedir)/arpack
+pkginclude_HEADERS = debug.h stat.h arpackdef.h
 
 AM_DISTCHECK_CONFIGURE_FLAGS =
 if MPI
 AM_DISTCHECK_CONFIGURE_FLAGS += --enable-mpi
 endif
-arpack_include_HEADERS = debug.h stat.h
 if ICB
 AM_DISTCHECK_CONFIGURE_FLAGS += --enable-icb
-arpack_include_HEADERS += debug_c.h debug_c.hpp
-arpack_include_HEADERS += stat_c.h stat_c.hpp
-arpack_include_HEADERS += arpack.h arpack.hpp
-if MPI
-arpack_include_HEADERS += parpack.h parpack.hpp
-endif
 endif
 
-SUBDIRS = UTIL SRC . TESTS EXAMPLES PARPACK
+SUBDIRS = .
+if ICB
+SUBDIRS += ICB
+endif
+SUBDIRS += UTIL SRC TESTS EXAMPLES
+if MPI
+SUBDIRS += PARPACK
+endif
 
 EXTRA_DIST = README.md PARPACK_CHANGES CHANGES DOCUMENTS VISUAL_STUDIO \
-detect_arpack_bug.m4 CMakeLists.txt arpack-ng-config.cmake.in arpack-ng-config-version.cmake.in \
-debug.h stat.h
-if ICB
-EXTRA_DIST += debug_c.h debug_c.hpp
-EXTRA_DIST += stat_c.h stat_c.hpp
-EXTRA_DIST += arpack.h arpack.hpp
-if MPI
-EXTRA_DIST += parpack.h parpack.hpp
-endif
-endif
+detect_arpack_bug.m4 CMakeLists.txt arpack-ng-config.cmake.in arpack-ng-config-version.cmake.in
 
 pkgconfig_DATA = arpack@LIBSUFFIX@.pc
 
-DISTCLEANFILES = $(pkgconfig_DATA)
+DISTCLEANFILES = $(pkgconfig_DATA) arpackdef.h

PARPACK/EXAMPLES/MPI/Makefile.am

@@ -9,10 +9,14 @@ CNDRV = pcndrv1
 ZNDRV = pzndrv1
 
 check_PROGRAMS = $(SNDRV) $(DNDRV) $(SSDRV) $(DSDRV) $(CNDRV) $(ZNDRV)
-# Disable tests for now. Needs to be run using mpirun
-#TESTS = $(check_PROGRAMS)
 
-EXTRA_DIST = debug.h stat.h
+# Run MPI tests with "mpirun -n 2"
+LOG_COMPILER = mpirun
+LOG_FLAGS = -n 2
+
+TESTS = $(check_PROGRAMS)
+
+EXTRA_DIST = debug.h stat.h README.CALLING-PARPACK-FROM-C-OR-CPP
 
 # Simple nonsymmetric problem using single precision
 psndrv1_SOURCES = psndrv1.f

PARPACK/EXAMPLES/MPI/pcndrv1.f

@@ -60,11 +60,17 @@ c
       include 'debug.h'
       include 'stat.h'
 c
-c     %---------------%
-c     | MPI INTERFACE |
-c     %---------------%
+c     %-------------------------------%
+c     | MPI INTERFACE                 |
+c     | ILP64 is not supported by MPI |
+c     | integer*4 must be imposed in  |
+c     | all calls involving MPI.      |
+c     |                               |
+c     | Use ierr for MPI calls.       |
+c     %-------------------------------%
+
+      integer*4         comm, myid, nprocs, rc, ierr
 
-      integer           comm, myid, nprocs, rc, nloc
 c     %-----------------------------%
 c     | Define maximum dimensions   |
 c     | for all arrays.             |
@@ -91,13 +97,16 @@ c
       Real
      &                  rwork(maxncv), rd(maxncv,3)
 c
-c     %---------------%
-c     | Local Scalars |
-c     %---------------%
+c     %------------------------------------%
+c     | Local Scalars                      |
+c     |                                    |
+c     | Use info if ILP64 can be supported |
+c     | (call to BLAS, LAPACK, ARPACK).    |
+c     %------------------------------------%
 c
       character         bmat*1, which*2
       integer           ido, n, nx, nev, ncv, lworkl, info, j,
-     &                  ierr, nconv, maxitr, ishfts, mode
+     &                  nloc, nconv, maxitr, ishfts, mode
       Complex
      &                  sigma
       Real
@@ -284,7 +293,7 @@ c
          call pcneupd (comm, rvec, 'A', select, d, v, ldv, sigma,
      &        workev, bmat, nloc, which, nev, tol, resid, ncv,
      &        v, ldv, iparam, ipntr, workd, workl, lworkl,
-     &        rwork, ierr)
+     &        rwork, info)
 c
 c        %----------------------------------------------%
 c        | Eigenvalues are returned in the one          |
@@ -297,7 +306,7 @@ c        | corresponding to the eigenvalues in D is     |
 c        | returned in V.                               |
 c        %----------------------------------------------%
 c
-         if ( ierr .ne. 0) then
+         if ( info .ne. 0) then
 c
 c           %------------------------------------%
 c           | Error condition:                   |
@@ -306,7 +315,7 @@ c           %------------------------------------%
 c
             if ( myid .eq. 0 ) then
                 print *, ' '
-                print *, ' Error with _neupd, info = ', ierr
+                print *, ' Error with _neupd, info = ', info
                 print *, ' Check the documentation of _neupd. '
                 print *, ' '
             endif
@@ -424,7 +433,7 @@ c----------------------------------------------------------------------------
 c
 c     .. MPI Declarations ...
       include           'mpif.h'
-      integer           comm, nprocs, myid, ierr,
+      integer*4         comm, nprocs, myid, ierr,
      &                  status(MPI_STATUS_SIZE)
 c
       integer           nloc, nx, np, j, lo, next, prev

PARPACK/EXAMPLES/MPI/pdndrv1.f

@@ -60,10 +60,17 @@ c
       include 'debug.h'
       include 'stat.h'
 
-c     %---------------%
-c     | MPI INTERFACE |
-c     %---------------%
-      integer           comm, myid, nprocs, rc, nloc
+c     %-------------------------------%
+c     | MPI INTERFACE                 |
+c     | ILP64 is not supported by MPI |
+c     | integer*4 must be imposed in  |
+c     | all calls involving MPI.      |
+c     |                               |
+c     | Use ierr for MPI calls.       |
+c     %-------------------------------%
+
+      integer*4         comm, myid, nprocs, rc, ierr
+
 c
 c     %-----------------------------%
 c     | Define maximum dimensions   |
@@ -90,13 +97,16 @@ c
      &                  workev(3*maxncv),
      &                  workl(3*maxncv*maxncv+6*maxncv)
 c
-c     %---------------%
-c     | Local Scalars |
-c     %---------------%
+c     %------------------------------------%
+c     | Local Scalars                      |
+c     |                                    |
+c     | Use info if ILP64 can be supported |
+c     | (call to BLAS, LAPACK, ARPACK).    |
+c     %------------------------------------%
 c
       character         bmat*1, which*2
       integer           ido, n, nx, nev, ncv, lworkl, info, j,
-     &                  ierr, nconv, maxitr, ishfts, mode
+     &                  nloc, nconv, maxitr, ishfts, mode
       Double precision
      &                  tol, sigmar, sigmai
       logical           first, rvec
@@ -295,7 +305,7 @@ c
          call pdneupd ( comm, rvec, 'A', select, d, d(1,2), v, ldv,
      &        sigmar, sigmai, workev, bmat, nloc, which, nev, tol,
      &        resid, ncv, v, ldv, iparam, ipntr, workd, workl,
-     &        lworkl, ierr )
+     &        lworkl, info )
 c
 c        %-----------------------------------------------%
 c        | The real part of the eigenvalue is returned   |
@@ -309,7 +319,7 @@ c        | for the invariant subspace corresponding to   |
 c        | the eigenvalues in D is returned in V.        |
 c        %-----------------------------------------------%
 c
-         if ( ierr .ne. 0) then
+         if ( info .ne. 0) then
 c
 c           %------------------------------------%
 c           | Error condition:                   |
@@ -318,7 +328,7 @@ c           %------------------------------------%
 c
          	if ( myid .eq. 0 ) then
              	print *, ' '
-             	print *, ' Error with _neupd, info = ', ierr
+             	print *, ' Error with _neupd, info = ', info
              	print *, ' Check the documentation of _neupd. '
              	print *, ' '
             endif
@@ -468,7 +478,7 @@ c----------------------------------------------------------------------------
 c
 c     .. MPI Declarations ...
       include           'mpif.h'
-      integer           comm, nprocs, myid, ierr,
+      integer*4         comm, nprocs, myid, ierr,
      &                  status(MPI_STATUS_SIZE)
 c
       integer           nloc, nx, np, j, lo, next, prev

PARPACK/EXAMPLES/MPI/pdndrv3.f

@@ -59,11 +59,17 @@ c
       include 'debug.h'
       include 'stat.h'
 
-c     %---------------%
-c     | MPI INTERFACE |
-c     %---------------%
+c     %-------------------------------%
+c     | MPI INTERFACE                 |
+c     | ILP64 is not supported by MPI |
+c     | integer*4 must be imposed in  |
+c     | all calls involving MPI.      |
+c     |                               |
+c     | Use ierr for MPI calls.       |
+c     %-------------------------------%
+
+      integer*4         comm, myid, nprocs, rc, ierr
 
-      integer           comm, myid, nprocs, rc, nloc
 c
 c     %-----------------------------%
 c     | Define leading dimensions   |
@@ -91,12 +97,15 @@ c
      &                  workl(3*maxncv*maxncv+6*maxncv),
      &                  md(maxn), me(maxn-1), temp(maxn), temp_buf(maxn)
 c
-c     %---------------%
-c     | Local Scalars |
-c     %---------------%
+c     %------------------------------------%
+c     | Local Scalars                      |
+c     |                                    |
+c     | Use info if ILP64 can be supported |
+c     | (call to BLAS, LAPACK, ARPACK).    |
+c     %------------------------------------%
 c
       character         bmat*1, which*2
-      integer           ido, n, nev, ncv, lworkl, info, ierr, j,
+      integer           ido, n, nev, ncv, lworkl, info, nloc, j,
      &                  nconv, maxitr, ishfts, mode, blk
       Double precision
      &                  tol, sigmar, sigmai
@@ -184,8 +193,8 @@ c
   20  continue
       md(n) = 4.0*one
 c
-      call dpttrf(n, md, me, ierr)
-      if ( ierr .ne. 0 ) then
+      call dpttrf(n, md, me, info)
+      if ( info .ne. 0 ) then
          print*, ' '
          print*, ' ERROR with _pttrf. '
          print*, ' '
@@ -265,8 +274,8 @@ c======== Hack for Linear system ======= ccc
             call MPI_ALLREDUCE( temp_buf, temp, n,
      &            MPI_DOUBLE_PRECISION, MPI_SUM, comm, ierr )
             call dpttrs(n, 1, md, me, temp, n,
-     &                  ierr)
-            if ( ierr .ne. 0 ) then
+     &                  info)
+            if ( info .ne. 0 ) then
                print*, ' '
                print*, ' ERROR with _pttrs. '
                print*, ' '
@@ -338,7 +347,7 @@ c
          call pdneupd ( comm, rvec, 'A', select, d, d(1,2), v, ldv,
      &        sigmar, sigmai, workev, bmat, nloc, which, nev, tol,
      &        resid, ncv, v, ldv, iparam, ipntr, workd,
-     &        workl, lworkl, ierr )
+     &        workl, lworkl, info )
 c
 c        %-----------------------------------------------%
 c        | The real part of the eigenvalue is returned   |
@@ -352,7 +361,7 @@ c        | for the invariant subspace corresponding to   |
 c        | the eigenvalues in D is returned in V.        |
 c        %-----------------------------------------------%
 c
-         if ( ierr .ne. 0 ) then
+         if ( info .ne. 0 ) then
 c
 c           %------------------------------------%
 c           | Error condition:                   |
@@ -361,7 +370,7 @@ c           %------------------------------------%
 c
             if ( myid .eq. 0 ) then
                print *, ' '
-               print *, ' Error with _neupd, info = ', ierr
+               print *, ' Error with _neupd, info = ', info
                print *, ' Check the documentation of _neupd'
                print *, ' '
             endif
@@ -499,7 +508,7 @@ c
 c
 c     .. MPI Declarations ...
       include           'mpif.h'
-      integer           comm, nprocs, myid, ierr,
+      integer*4         comm, nprocs, myid, ierr,
      &                  status(MPI_STATUS_SIZE)
 c
       integer           nloc, n, j, next, prev
@@ -557,7 +566,7 @@ c
 c
 c     .. MPI Declarations ...
       include           'mpif.h'
-      integer           comm, nprocs, myid, ierr,
+      integer*4         comm, nprocs, myid, ierr,
      &                  status(MPI_STATUS_SIZE)
 c
       integer           nloc, j, next, prev

PARPACK/EXAMPLES/MPI/pdsdrv1.f

@@ -57,11 +57,17 @@ c
       include 'debug.h'
       include 'stat.h'
 
-c     %---------------%
-c     | MPI INTERFACE |
-c     %---------------%
+c     %-------------------------------%
+c     | MPI INTERFACE                 |
+c     | ILP64 is not supported by MPI |
+c     | integer*4 must be imposed in  |
+c     | all calls involving MPI.      |
+c     |                               |
+c     | Use ierr for MPI calls.       |
+c     %-------------------------------%
+
+      integer*4         comm, myid, nprocs, rc, ierr
 
-      integer           comm, myid, nprocs, rc, nloc
 c
 c     %-----------------------------%
 c     | Define leading dimensions   |
@@ -87,12 +93,15 @@ c
       logical          select(maxncv)
       integer          iparam(11), ipntr(11)
 c
-c     %---------------%
-c     | Local Scalars |
-c     %---------------%
+c     %------------------------------------%
+c     | Local Scalars                      |
+c     |                                    |
+c     | Use info if ILP64 can be supported |
+c     | (call to BLAS, LAPACK, ARPACK).    |
+c     %------------------------------------%
 c
       character        bmat*1, which*2
-      integer          ido, n, nev, ncv, lworkl, info, ierr, j,
+      integer          ido, n, nev, ncv, lworkl, info, nloc, j,
      &                 nx, nconv, maxitr, mode, ishfts
       logical          rvec
       Double precision
@@ -293,7 +302,7 @@ c
          call pdseupd ( comm, rvec, 'All', select,
      &        d, v, ldv, sigma,
      &        bmat, nloc, which, nev, tol, resid, ncv, v, ldv,
-     &        iparam, ipntr, workd, workl, lworkl, ierr )
+     &        iparam, ipntr, workd, workl, lworkl, info )
 c        %----------------------------------------------%
 c        | Eigenvalues are returned in the first column |
 c        | of the two dimensional array D and the       |
@@ -305,7 +314,7 @@ c        | corresponding to the eigenvalues in D is     |
 c        | returned in V.                               |
 c        %----------------------------------------------%
 c
-         if ( ierr .ne. 0) then
+         if ( info .ne. 0) then
 c
 c            %------------------------------------%
 c            | Error condition:                   |
@@ -315,7 +324,7 @@ c
 c
             if ( myid .eq. 0 ) then
              	print *, ' '
-             	print *, ' Error with _seupd, info = ', ierr
+             	print *, ' Error with _seupd, info = ', info
              	print *, ' Check the documentation of _seupd. '
              	print *, ' '
             endif
@@ -421,7 +430,7 @@ c
 c
 c     .. MPI Declarations ...
       include           'mpif.h'
-      integer           comm, nprocs, myid, ierr,
+      integer*4         comm, nprocs, myid, ierr,
      &                  status(MPI_STATUS_SIZE)
       integer           nloc, nx, np, j, lo, next, prev
       Double precision

PARPACK/EXAMPLES/MPI/psndrv1.f

@@ -60,10 +60,17 @@ c
       include 'debug.h'
       include 'stat.h'
 
-c     %---------------%
-c     | MPI INTERFACE |
-c     %---------------%
-      integer           comm, myid, nprocs, rc, nloc
+c     %-------------------------------%
+c     | MPI INTERFACE                 |
+c     | ILP64 is not supported by MPI |
+c     | integer*4 must be imposed in  |
+c     | all calls involving MPI.      |
+c     |                               |
+c     | Use ierr for MPI calls.       |
+c     %-------------------------------%
+
+      integer*4         comm, myid, nprocs, rc, ierr
+
 c
 c     %-----------------------------%
 c     | Define maximum dimensions   |
@@ -90,13 +97,16 @@ c
      &                  workev(3*maxncv),
      &                  workl(3*maxncv*maxncv+6*maxncv)
 c
-c     %---------------%
-c     | Local Scalars |
-c     %---------------%
+c     %------------------------------------%
+c     | Local Scalars                      |
+c     |                                    |
+c     | Use info if ILP64 can be supported |
+c     | (call to BLAS, LAPACK, ARPACK).    |
+c     %------------------------------------%
 c
       character         bmat*1, which*2
       integer           ido, n, nx, nev, ncv, lworkl, info, j,
-     &                  ierr, nconv, maxitr, ishfts, mode
+     &                  nloc, nconv, maxitr, ishfts, mode
       Real
      &                  tol, sigmar, sigmai
       logical           first, rvec
@@ -295,7 +305,7 @@ c
          call psneupd ( comm, rvec, 'A', select, d, d(1,2), v, ldv,
      &        sigmar, sigmai, workev, bmat, nloc, which, nev, tol,
      &        resid, ncv, v, ldv, iparam, ipntr, workd, workl,
-     &        lworkl, ierr )
+     &        lworkl, info )
 c
 c        %-----------------------------------------------%
 c        | The real part of the eigenvalue is returned   |
@@ -309,7 +319,7 @@ c        | for the invariant subspace corresponding to   |
 c        | the eigenvalues in D is returned in V.        |
 c        %-----------------------------------------------%
 c
-         if ( ierr .ne. 0) then
+         if ( info .ne. 0) then
 c
 c           %------------------------------------%
 c           | Error condition:                   |
@@ -318,7 +328,7 @@ c           %------------------------------------%
 c
          	if ( myid .eq. 0 ) then
              	print *, ' '
-             	print *, ' Error with _neupd, info = ', ierr
+             	print *, ' Error with _neupd, info = ', info
              	print *, ' Check the documentation of _neupd. '
              	print *, ' '
             endif
@@ -468,7 +478,7 @@ c----------------------------------------------------------------------------
 c
 c     .. MPI Declarations ...
       include           'mpif.h'
-      integer           comm, nprocs, myid, ierr,
+      integer*4         comm, nprocs, myid, ierr,
      &                  status(MPI_STATUS_SIZE)
 c
       integer           nloc, nx, np, j, lo, next, prev

PARPACK/EXAMPLES/MPI/psndrv3.f

@@ -59,11 +59,17 @@ c
       include 'debug.h'
       include 'stat.h'
 
-c     %---------------%
-c     | MPI INTERFACE |
-c     %---------------%
+c     %-------------------------------%
+c     | MPI INTERFACE                 |
+c     | ILP64 is not supported by MPI |
+c     | integer*4 must be imposed in  |
+c     | all calls involving MPI.      |
+c     |                               |
+c     | Use ierr for MPI calls.       |
+c     %-------------------------------%
+
+      integer*4         comm, myid, nprocs, rc, ierr
 
-      integer           comm, myid, nprocs, rc, nloc
 c
 c     %-----------------------------%
 c     | Define leading dimensions   |
@@ -91,12 +97,15 @@ c
      &                  workl(3*maxncv*maxncv+6*maxncv),
      &                  md(maxn), me(maxn-1), temp(maxn), temp_buf(maxn)
 c
-c     %---------------%
-c     | Local Scalars |
-c     %---------------%
+c     %------------------------------------%
+c     | Local Scalars                      |
+c     |                                    |
+c     | Use info if ILP64 can be supported |
+c     | (call to BLAS, LAPACK, ARPACK).    |
+c     %------------------------------------%
 c
       character         bmat*1, which*2
-      integer           ido, n, nev, ncv, lworkl, info, ierr, j,
+      integer           ido, n, nev, ncv, lworkl, info, nloc, j,
      &                  nconv, maxitr, ishfts, mode, blk
       Real
      &                  tol, sigmar, sigmai
@@ -184,8 +193,8 @@ c
   20  continue
       md(n) = 4.0*one
 c
-      call spttrf(n, md, me, ierr)
-      if ( ierr .ne. 0 ) then
+      call spttrf(n, md, me, info)
+      if ( info .ne. 0 ) then
          print*, ' '
          print*, ' ERROR with _pttrf. '
          print*, ' '
@@ -265,8 +274,8 @@ c======== Hack for Linear system ======= ccc
             call MPI_ALLREDUCE( temp_buf, temp, n,
      &            MPI_REAL, MPI_SUM, comm, ierr )
             call spttrs(n, 1, md, me, temp, n,
-     &                  ierr)
-            if ( ierr .ne. 0 ) then
+     &                  info)
+            if ( info .ne. 0 ) then
                print*, ' '
                print*, ' ERROR with _pttrs. '
                print*, ' '
@@ -338,7 +347,7 @@ c
          call psneupd ( comm, rvec, 'A', select, d, d(1,2), v, ldv,
      &        sigmar, sigmai, workev, bmat, nloc, which, nev, tol,
      &        resid, ncv, v, ldv, iparam, ipntr, workd,
-     &        workl, lworkl, ierr )
+     &        workl, lworkl, info )
 c
 c        %-----------------------------------------------%
 c        | The real part of the eigenvalue is returned   |
@@ -352,7 +361,7 @@ c        | for the invariant subspace corresponding to   |
 c        | the eigenvalues in D is returned in V.        |
 c        %-----------------------------------------------%
 c
-         if ( ierr .ne. 0 ) then
+         if ( info .ne. 0 ) then
 c
 c           %------------------------------------%
 c           | Error condition:                   |
@@ -361,7 +370,7 @@ c           %------------------------------------%
 c
             if ( myid .eq. 0 ) then
                print *, ' '
-               print *, ' Error with _neupd, info = ', ierr
+               print *, ' Error with _neupd, info = ', info
                print *, ' Check the documentation of _neupd'
                print *, ' '
             endif
@@ -499,7 +508,7 @@ c
 c
 c     .. MPI Declarations ...
       include           'mpif.h'
-      integer           comm, nprocs, myid, ierr,
+      integer*4         comm, nprocs, myid, ierr,
      &                  status(MPI_STATUS_SIZE)
 c
       integer           nloc, n, j, next, prev
@@ -557,7 +566,7 @@ c
 c
 c     .. MPI Declarations ...
       include           'mpif.h'
-      integer           comm, nprocs, myid, ierr,
+      integer*4         comm, nprocs, myid, ierr,
      &                  status(MPI_STATUS_SIZE)
 c
       integer           nloc, j, next, prev

PARPACK/EXAMPLES/MPI/pssdrv1.f

@@ -57,11 +57,17 @@ c
       include 'debug.h'
       include 'stat.h'
 
-c     %---------------%
-c     | MPI INTERFACE |
-c     %---------------%
+c     %-------------------------------%
+c     | MPI INTERFACE                 |
+c     | ILP64 is not supported by MPI |
+c     | integer*4 must be imposed in  |
+c     | all calls involving MPI.      |
+c     |                               |
+c     | Use ierr for MPI calls.       |
+c     %-------------------------------%
+
+      integer*4         comm, myid, nprocs, rc, ierr
 
-      integer           comm, myid, nprocs, rc, nloc
 c
 c     %-----------------------------%
 c     | Define leading dimensions   |
@@ -87,12 +93,15 @@ c
       logical          select(maxncv)
       integer          iparam(11), ipntr(11)
 c
-c     %---------------%
-c     | Local Scalars |
-c     %---------------%
+c     %------------------------------------%
+c     | Local Scalars                      |
+c     |                                    |
+c     | Use info if ILP64 can be supported |
+c     | (call to BLAS, LAPACK, ARPACK).    |
+c     %------------------------------------%
 c
       character        bmat*1, which*2
-      integer          ido, n, nev, ncv, lworkl, info, ierr, j,
+      integer          ido, n, nev, ncv, lworkl, info, nloc, j,
      &                 nx, nconv, maxitr, mode, ishfts
       logical          rvec
       Real
@@ -293,7 +302,7 @@ c
          call psseupd ( comm, rvec, 'All', select,
      &        d, v, ldv, sigma,
      &        bmat, nloc, which, nev, tol, resid, ncv, v, ldv,
-     &        iparam, ipntr, workd, workl, lworkl, ierr )
+     &        iparam, ipntr, workd, workl, lworkl, info )
 c        %----------------------------------------------%
 c        | Eigenvalues are returned in the first column |
 c        | of the two dimensional array D and the       |
@@ -305,7 +314,7 @@ c        | corresponding to the eigenvalues in D is     |
 c        | returned in V.                               |
 c        %----------------------------------------------%
 c
-         if ( ierr .ne. 0) then
+         if ( info .ne. 0) then
 c
 c            %------------------------------------%
 c            | Error condition:                   |
@@ -315,7 +324,7 @@ c
 c
             if ( myid .eq. 0 ) then
              	print *, ' '
-             	print *, ' Error with _seupd, info = ', ierr
+             	print *, ' Error with _seupd, info = ', info
              	print *, ' Check the documentation of _seupd. '
              	print *, ' '
             endif
@@ -421,7 +430,7 @@ c
 c
 c     .. MPI Declarations ...
       include           'mpif.h'
-      integer           comm, nprocs, myid, ierr,
+      integer*4         comm, nprocs, myid, ierr,
      &                  status(MPI_STATUS_SIZE)
       integer           nloc, nx, np, j, lo, next, prev
       Real

PARPACK/EXAMPLES/MPI/pzndrv1.f

@@ -60,11 +60,17 @@ c
       include 'debug.h'
       include 'stat.h'
 c
-c     %---------------%
-c     | MPI INTERFACE |
-c     %---------------%
+c     %-------------------------------%
+c     | MPI INTERFACE                 |
+c     | ILP64 is not supported by MPI |
+c     | integer*4 must be imposed in  |
+c     | all calls involving MPI.      |
+c     |                               |
+c     | Use ierr for MPI calls.       |
+c     %-------------------------------%
+
+      integer*4         comm, myid, nprocs, rc, ierr
 
-      integer           comm, myid, nprocs, rc, nloc
 c     %-----------------------------%
 c     | Define maximum dimensions   |
 c     | for all arrays.             |
@@ -91,13 +97,16 @@ c
       Double precision
      &                  rwork(maxncv), rd(maxncv,3)
 c
-c     %---------------%
-c     | Local Scalars |
-c     %---------------%
+c     %------------------------------------%
+c     | Local Scalars                      |
+c     |                                    |
+c     | Use info if ILP64 can be supported |
+c     | (call to BLAS, LAPACK, ARPACK).    |
+c     %------------------------------------%
 c
       character         bmat*1, which*2
       integer           ido, n, nx, nev, ncv, lworkl, info, j,
-     &                  ierr, nconv, maxitr, ishfts, mode
+     &                  nloc, nconv, maxitr, ishfts, mode
       Complex*16
      &                  sigma
       Double precision
@@ -284,7 +293,7 @@ c
          call pzneupd (comm, rvec, 'A', select, d, v, ldv, sigma,
      &        workev, bmat, nloc, which, nev, tol, resid, ncv,
      &        v, ldv, iparam, ipntr, workd, workl, lworkl,
-     &        rwork, ierr)
+     &        rwork, info)
 c
 c        %----------------------------------------------%
 c        | Eigenvalues are returned in the one          |
@@ -297,7 +306,7 @@ c        | corresponding to the eigenvalues in D is     |
 c        | returned in V.                               |
 c        %----------------------------------------------%
 c
-         if ( ierr .ne. 0) then
+         if ( info .ne. 0) then
 c
 c           %------------------------------------%
 c           | Error condition:                   |
@@ -306,7 +315,7 @@ c           %------------------------------------%
 c
             if ( myid .eq. 0 ) then
                 print *, ' '
-                print *, ' Error with _neupd, info = ', ierr
+                print *, ' Error with _neupd, info = ', info
                 print *, ' Check the documentation of _neupd. '
                 print *, ' '
             endif
@@ -424,7 +433,7 @@ c----------------------------------------------------------------------------
 c
 c     .. MPI Declarations ...
       include           'mpif.h'
-      integer           comm, nprocs, myid, ierr,
+      integer*4         comm, nprocs, myid, ierr,
      &                  status(MPI_STATUS_SIZE)
 c
       integer           nloc, nx, np, j, lo, next, prev

PARPACK/SRC/BLACS/Makefile.am

@@ -20,4 +20,4 @@ libparpack@LIBSUFFIX@_la_LIBADD = \
 	$(top_builddir)/PARPACK/UTIL/BLACS/libparpackutil.la \
 	$(top_builddir)/SRC/libarpack@LIBSUFFIX@.la \
 	$(LAPACK_LIBS) $(BLAS_LIBS)
-libparpack@LIBSUFFIX@_la_LDFLAGS = -no-undefined -version-info 2:0
+libparpack@LIBSUFFIX@_la_LDFLAGS = -no-undefined -version-info 3:0:1

PARPACK/SRC/MPI/Makefile.am

@@ -1,4 +1,4 @@
-AUTOMAKE_OPTIONS = subdir-objects # Needed as debug_init/icb.f90 are not in current directory.
+AUTOMAKE_OPTIONS = subdir-objects # Needed as debug/stat* are not in current directory.
 
 F77 = $(MPIF77)
 FFLAGS_SAV = @FFLAGS@
@@ -20,10 +20,10 @@ ZSRC = pznaitr.f pznapps.f pznaup2.f pznaupd.f pzneigh.f pzneupd.f pzngets.f \
        pzgetv0.f pzlarnv.f pdznorm2.f
 
 if ICB
-SSRC += icbpss.f90 icbpsn.f90
-DSRC += icbpds.f90 icbpdn.f90
-CSRC += icbpcn.f90
-ZSRC += icbpzn.f90
+SSRC += icbpss.F90 icbpsn.F90
+DSRC += icbpds.F90 icbpdn.F90
+CSRC += icbpcn.F90
+ZSRC += icbpzn.F90
 endif
 
 EXTRA_DIST = debug.h stat.h pcontext.h
@@ -34,14 +34,15 @@ libparpack@LIBSUFFIX@_noopt_la_FFLAGS = -O0
 
 lib_LTLIBRARIES = libparpack@LIBSUFFIX@.la
 libparpack@LIBSUFFIX@_la_SOURCES = $(PSRC) $(SSRC) $(DSRC) $(CSRC) $(ZSRC)
-libparpack@LIBSUFFIX@_la_SOURCES += $(top_builddir)/debug_init.f90
-if ICB
-libparpack@LIBSUFFIX@_la_SOURCES += $(top_builddir)/debug_icb.f90
-libparpack@LIBSUFFIX@_la_SOURCES += $(top_builddir)/stat_icb.f90
-endif
+libparpack@LIBSUFFIX@_la_SOURCES += $(top_builddir)/dbgini.f
+libparpack@LIBSUFFIX@_la_SOURCES += $(top_builddir)/staini.f
 libparpack@LIBSUFFIX@_la_FFLAGS = $(FFLAGS_SAV)
 libparpack@LIBSUFFIX@_la_LIBADD = libparpack@LIBSUFFIX@_noopt.la \
 	$(top_builddir)/PARPACK/UTIL/MPI/libparpackutil.la \
 	$(top_builddir)/SRC/libarpack@LIBSUFFIX@.la \
 	$(LAPACK_LIBS) $(BLAS_LIBS) $(MPI_Fortran_LIBS)
-libparpack@LIBSUFFIX@_la_LDFLAGS = -no-undefined -version-info 2:0
+libparpack@LIBSUFFIX@_la_LDFLAGS = -no-undefined -version-info 3:0:1
+if ICB
+libparpack@LIBSUFFIX@_la_LIBADD += $(top_builddir)/ICB/libdbgicb.la $(top_builddir)/ICB/libstaicb.la
+libparpack@LIBSUFFIX@_la_CPPFLAGS = $(AM_CPPFLAGS) -I$(top_builddir)
+endif

PARPACK/SRC/MPI/icbpcn.F90

@@ -5,6 +5,7 @@ subroutine pcnaupd_c(comm, ido, bmat, n, which, nev, tol, resid, ncv, v, ldv,&
                      bind(c, name="pcnaupd_c")
   use :: iso_c_binding
   implicit none
+#include "arpackdef.h"
   integer(kind=c_int),          value,              intent(in)    :: comm
   integer(kind=c_int),                              intent(inout) :: ido
   character(kind=c_char),       dimension(1),       intent(in)    :: bmat
@@ -33,6 +34,7 @@ subroutine pcneupd_c(comm, rvec, howmny, select, d, z, ldz, sigma, workev,&
                      bind(c, name="pcneupd_c")
   use :: iso_c_binding
   implicit none
+#include "arpackdef.h"
   integer(kind=c_int),          value,              intent(in)    :: comm
   logical(kind=c_bool),         value,              intent(in)    :: rvec
   character(kind=c_char),       dimension(1),       intent(in)    :: howmny

PARPACK/SRC/MPI/icbpdn.F90

@@ -5,6 +5,7 @@ subroutine pdnaupd_c(comm, ido, bmat, n, which, nev, tol, resid, ncv, v, ldv,&
                      bind(c, name="pdnaupd_c")
   use :: iso_c_binding
   implicit none
+#include "arpackdef.h"
   integer(kind=c_int),    value,               intent(in)    :: comm
   integer(kind=c_int),                         intent(inout) :: ido
   character(kind=c_char), dimension(1),        intent(in)    :: bmat
@@ -27,12 +28,13 @@ subroutine pdnaupd_c(comm, ido, bmat, n, which, nev, tol, resid, ncv, v, ldv,&
 end subroutine pdnaupd_c
 
 subroutine pdneupd_c(comm, rvec, howmny, select,                  &
-                     dr, di, z, ldz, sigmar, sigmai,              &
+                     dr, di, z, ldz, sigmar, sigmai, workev,      &
                      bmat, n, which, nev, tol, resid, ncv, v, ldv,&
                      iparam, ipntr, workd, workl, lworkl, info)   &
                      bind(c, name="pdneupd_c")
   use :: iso_c_binding
   implicit none
+#include "arpackdef.h"
   integer(kind=c_int),    value,               intent(in)    :: comm
   logical(kind=c_bool),   value,               intent(in)    :: rvec
   character(kind=c_char), dimension(1),        intent(in)    :: howmny
@@ -43,6 +45,7 @@ subroutine pdneupd_c(comm, rvec, howmny, select,                  &
   integer(kind=c_int),    value,               intent(in)    :: ldz
   real(kind=c_double),    value,               intent(in)    :: sigmar
   real(kind=c_double),    value,               intent(in)    :: sigmai
+  real(kind=c_double),    dimension(3*ncv),    intent(out)   :: workev
   character(kind=c_char), dimension(1),        intent(in)    :: bmat
   integer(kind=c_int),    value,               intent(in)    :: n
   character(kind=c_char), dimension(2),        intent(in)    :: which
@@ -59,7 +62,7 @@ subroutine pdneupd_c(comm, rvec, howmny, select,                  &
   integer(kind=c_int),    value,               intent(in)    :: lworkl
   integer(kind=c_int),                         intent(inout) :: info
   call pdneupd(comm, rvec, howmny, select,                  &
-               dr, di, z, ldz, sigmar, sigmai,              &
+               dr, di, z, ldz, sigmar, sigmai, workev,      &
                bmat, n, which, nev, tol, resid, ncv, v, ldv,&
                iparam, ipntr, workd, workl, lworkl, info)
 end subroutine pdneupd_c

PARPACK/SRC/MPI/icbpds.F90

@@ -5,6 +5,7 @@ subroutine pdsaupd_c(comm, ido, bmat, n, which, nev, tol, resid, ncv, v, ldv,&
                      bind(c, name="pdsaupd_c")
   use :: iso_c_binding
   implicit none
+#include "arpackdef.h"
   integer(kind=c_int),    value,               intent(in)    :: comm
   integer(kind=c_int),                         intent(inout) :: ido
   character(kind=c_char), dimension(1),        intent(in)    :: bmat
@@ -32,6 +33,7 @@ subroutine pdseupd_c(comm, rvec, howmny, select, d, z, ldz, sigma,&
                      bind(c, name="pdseupd_c")
   use :: iso_c_binding
   implicit none
+#include "arpackdef.h"
   integer(kind=c_int),    value,               intent(in)    :: comm
   logical(kind=c_bool),   value,               intent(in)    :: rvec
   character(kind=c_char), dimension(1),        intent(in)    :: howmny

PARPACK/SRC/MPI/icbpsn.F90

@@ -5,6 +5,7 @@ subroutine psnaupd_c(comm, ido, bmat, n, which, nev, tol, resid, ncv, v, ldv,&
                      bind(c, name="psnaupd_c")
   use :: iso_c_binding
   implicit none
+#include "arpackdef.h"
   integer(kind=c_int),    value,               intent(in)    :: comm
   integer(kind=c_int),                         intent(inout) :: ido
   character(kind=c_char), dimension(1),        intent(in)    :: bmat
@@ -27,12 +28,13 @@ subroutine psnaupd_c(comm, ido, bmat, n, which, nev, tol, resid, ncv, v, ldv,&
 end subroutine psnaupd_c
 
 subroutine psneupd_c(comm, rvec, howmny, select,                  &
-                     dr, di, z, ldz, sigmar, sigmai,              &
+                     dr, di, z, ldz, sigmar, sigmai, workev,      &
                      bmat, n, which, nev, tol, resid, ncv, v, ldv,&
                      iparam, ipntr, workd, workl, lworkl, info)   &
                      bind(c, name="psneupd_c")
   use :: iso_c_binding
   implicit none
+#include "arpackdef.h"
   integer(kind=c_int),    value,               intent(in)    :: comm
   logical(kind=c_bool),   value,               intent(in)    :: rvec
   character(kind=c_char), dimension(1),        intent(in)    :: howmny
@@ -43,6 +45,7 @@ subroutine psneupd_c(comm, rvec, howmny, select,                  &
   integer(kind=c_int),    value,               intent(in)    :: ldz
   real(kind=c_float),     value,               intent(in)    :: sigmar
   real(kind=c_float),     value,               intent(in)    :: sigmai
+  real(kind=c_float),     dimension(3*ncv),    intent(out)   :: workev
   character(kind=c_char), dimension(1),        intent(in)    :: bmat
   integer(kind=c_int),    value,               intent(in)    :: n
   character(kind=c_char), dimension(2),        intent(in)    :: which
@@ -59,7 +62,7 @@ subroutine psneupd_c(comm, rvec, howmny, select,                  &
   integer(kind=c_int),    value,               intent(in)    :: lworkl
   integer(kind=c_int),                         intent(inout) :: info
   call psneupd(comm, rvec, howmny, select,                  &
-               dr, di, z, ldz, sigmar, sigmai,              &
+               dr, di, z, ldz, sigmar, sigmai, workev,      &
                bmat, n, which, nev, tol, resid, ncv, v, ldv,&
                iparam, ipntr, workd, workl, lworkl, info)
 end subroutine psneupd_c

PARPACK/SRC/MPI/icbpss.F90

@@ -5,6 +5,7 @@ subroutine pssaupd_c(comm, ido, bmat, n, which, nev, tol, resid, ncv, v, ldv,&
                      bind(c, name="pssaupd_c")
   use :: iso_c_binding
   implicit none
+#include "arpackdef.h"
   integer(kind=c_int),    value,               intent(in)    :: comm
   integer(kind=c_int),                         intent(inout) :: ido
   character(kind=c_char), dimension(1),        intent(in)    :: bmat
@@ -32,6 +33,7 @@ subroutine psseupd_c(comm, rvec, howmny, select, d, z, ldz, sigma,&
                      bind(c, name="psseupd_c")
   use :: iso_c_binding
   implicit none
+#include "arpackdef.h"
   integer(kind=c_int),    value,               intent(in)    :: comm
   logical(kind=c_bool),   value,               intent(in)    :: rvec
   character(kind=c_char), dimension(1),        intent(in)    :: howmny

PARPACK/SRC/MPI/icbpzn.F90

@@ -5,6 +5,7 @@ subroutine pznaupd_c(comm, ido, bmat, n, which, nev, tol, resid, ncv, v, ldv,&
                      bind(c, name="pznaupd_c")
   use :: iso_c_binding
   implicit none
+#include "arpackdef.h"
   integer(kind=c_int),              value,              intent(in)    :: comm
   integer(kind=c_int),                                  intent(inout) :: ido
   character(kind=c_char),           dimension(1),       intent(in)    :: bmat
@@ -33,6 +34,7 @@ subroutine pzneupd_c(comm, rvec, howmny, select, d, z, ldz, sigma, workev,&
                      bind(c, name="pzneupd_c")
   use :: iso_c_binding
   implicit none
+#include "arpackdef.h"
   integer(kind=c_int),              value,              intent(in)    :: comm
   logical(kind=c_bool),             value,              intent(in)    :: rvec
   character(kind=c_char),           dimension(1),       intent(in)    :: howmny

PARPACK/TESTS/MPI/Makefile.am

@@ -20,11 +20,11 @@ issue46_SOURCES= issue46.f
 if ICB
 icb_parpack_c_SOURCES = icb_parpack_c.c
 icb_parpack_c_LDADD = $(LDADD) $(MPI_C_LIBS)
-icb_parpack_c_CFLAGS = -I$(top_builddir)
+icb_parpack_c_CPPFLAGS = $(AM_CPPFLAGS) -I$(top_builddir) -I$(top_srcdir)/ICB
 
 icb_parpack_cpp_SOURCES = icb_parpack_cpp.cpp
 icb_parpack_cpp_LDADD = $(LDADD) $(MPI_CXX_LIBS)
-icb_parpack_cpp_CPPFLAGS = -I$(top_builddir)
+icb_parpack_cpp_CPPFLAGS = $(AM_CPPFLAGS) -I$(top_builddir) -I$(top_srcdir)/ICB
 endif
 
 TESTS = $(check_PROGRAMS)

PARPACK/TESTS/MPI/icb_parpack_c.c

@@ -21,9 +21,6 @@
  * A*x = lambda*x where A is the diagonal matrix
  * with entries 1000, 999, ... , 2, 1 on the diagonal.
  * */
-#ifndef BLASINT
-#define BLASINT int
-#endif
 
 void dMatVec(double * x, double * y) {
   int i;
@@ -32,25 +29,25 @@ void dMatVec(double * x, double * y) {
 };
 
 int ds() {
-  BLASINT ido = 0;
+  a_int ido = 0;
   char bmat[] = "I";
-  BLASINT N = 1000;
+  a_int N = 1000;
   char which[] = "LM";
-  BLASINT nev = 3;
+  a_int nev = 3;
   double tol = 0;
   double resid[N];
-  BLASINT ncv = 2*nev+1;
+  a_int ncv = 2*nev+1;
   double V[ncv*N];
-  BLASINT ldv = N;
-  BLASINT iparam[11];
-  BLASINT ipntr[14];
+  a_int ldv = N;
+  a_int iparam[11];
+  a_int ipntr[14];
   double workd[3*N];
   bool rvec = true;
   char howmny[] = "A";
   double* d = (double*) malloc((nev+1)*sizeof(double));
-  int select[ncv];
+  a_int select[ncv];
   double z[(N+1)*(nev+1)];
-  BLASINT ldz = N+1;
+  a_int ldz = N+1;
   double sigma=0;
   int k;
   for (k=0; k < 3*N; ++k )
@@ -58,8 +55,8 @@ int ds() {
   double workl[3*(ncv*ncv) + 6*ncv];
   for (k=0; k < 3*(ncv*ncv) + 6*ncv; ++k )
     workl[k] = 0;
-  BLASINT lworkl = 3*(ncv*ncv) + 6*ncv;
-  BLASINT info = 0;
+  a_int lworkl = 3*(ncv*ncv) + 6*ncv;
+  a_int info = 0;
   int rank; MPI_Comm_rank(MPI_COMM_WORLD, &rank);
 
   iparam[0] = 1;
@@ -101,36 +98,36 @@ void zMatVec(double _Complex * x, double _Complex * y) {
 };
 
 int zn() {
-  BLASINT ido = 0;
+  a_int ido = 0;
   char bmat[] = "I";
-  BLASINT N = 1000;
+  a_int N = 1000;
   char which[] = "LM";
-  BLASINT nev = 1;
+  a_int nev = 1;
   double tol = 0;
   double _Complex resid[N];
-  BLASINT ncv = 2*nev+1;
+  a_int ncv = 2*nev+1;
   double _Complex V[ncv*N];
-  BLASINT ldv = N;
-  BLASINT iparam[11];
-  BLASINT ipntr[14];
+  a_int ldv = N;
+  a_int iparam[11];
+  a_int ipntr[14];
   double _Complex workd[3*N];
   bool rvec = true;
   char howmny[] = "A";
   double _Complex* d = (double _Complex*) malloc((nev+1)*sizeof(double _Complex));
-  int select[ncv];
+  a_int select[ncv];
   double _Complex z[(N+1)*(nev+1)];
-  BLASINT ldz = N+1;
-  double sigma=0;
+  a_int ldz = N+1;
+  double _Complex sigma=0. + I*0.;
   int k;
   for (k=0; k < 3*N; ++k )
     workd[k] = 0. + I * 0.;
   double _Complex workl[3*(ncv*ncv) + 6*ncv];
   for (k=0; k < 3*(ncv*ncv) + 6*ncv; ++k )
     workl[k] = 0. + I * 0.;
-  BLASINT lworkl = 3*(ncv*ncv) + 6*ncv;
+  a_int lworkl = 3*(ncv*ncv) + 6*ncv;
   double _Complex rwork[ncv];
   double _Complex workev[2*ncv];
-  BLASINT info = 0;
+  a_int info = 0;
   int rank; MPI_Comm_rank(MPI_COMM_WORLD, &rank);
 
   iparam[0] = 1;

PARPACK/TESTS/MPI/icb_parpack_cpp.cpp

@@ -21,10 +21,6 @@
 #include "debug_c.hpp"  // debug parpack.
 #include "stat_c.hpp"   // arpack statistics.
 
-#ifndef BLASINT
-#define BLASINT int
-#endif
-
 void diagonal_matrix_vector_product(float const* const x, float* const y) {
   for (int i = 0; i < 1000; ++i) {
     y[i] = static_cast<float>(i + 1) * x[i];
@@ -32,18 +28,18 @@ void diagonal_matrix_vector_product(float const* const x, float* const y) {
 }
 
 void real_symmetric_runner() {
-  BLASINT N = 1000;
-  BLASINT nev = 3;
-  BLASINT ncv = 2 * nev + 1;
-  BLASINT ldz = N + 1;
-  BLASINT lworkl = 3 * (ncv * ncv) + 6 * ncv;
-  BLASINT ldv = N;
+  a_int N = 1000;
+  a_int nev = 3;
+  a_int ncv = 2 * nev + 1;
+  a_int ldz = N + 1;
+  a_int lworkl = 3 * (ncv * ncv) + 6 * ncv;
+  a_int ldv = N;
 
   bool rvec = true;
   float tol = 0.0f;
   float sigma = 0.0f;
 
-  std::array<BLASINT, 14> ipntr;
+  std::array<a_int, 14> ipntr;
 
   std::vector<float> workd(3 * N, 0.0f);
   std::vector<float> workl(3 * (ncv * ncv) + 6 * ncv, 0.0f);
@@ -51,14 +47,14 @@ void real_symmetric_runner() {
   std::vector<float> d(nev + 1);
   std::vector<float> z((N + 1) * (nev + 1));
   std::vector<float> resid(N);
-  std::vector<int> select(ncv);
+  std::vector<a_int> select(ncv);
 
-  BLASINT info = 0;
+  a_int info = 0;
 
   int rank;
   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
 
-  std::array<BLASINT, 11> iparam;
+  std::array<a_int, 11> iparam;
   iparam[0] = 1;
   iparam[2] = 10 * N;
   iparam[3] = 1;
@@ -67,7 +63,7 @@ void real_symmetric_runner() {
 
   MPI_Fint MCW = MPI_Comm_c2f(MPI_COMM_WORLD);
 
-  BLASINT ido = 0;
+  a_int ido = 0;
 
   while (ido != 99) {
     arpack::saupd(MCW, ido, arpack::bmat::identity, N,
@@ -105,46 +101,46 @@ void diagonal_matrix_vector_product(std::complex<float> const* const x,
 }
 
 void complex_symmetric_runner() {
-  BLASINT N = 1000;
-  BLASINT nev = 1;
-  BLASINT ncv = 2 * nev + 1;
-  BLASINT ldv = N;
-  BLASINT ldz = N + 1;
+  a_int N = 1000;
+  a_int nev = 1;
+  a_int ncv = 2 * nev + 1;
+  a_int ldv = N;
+  a_int ldz = N + 1;
 
   float tol = 0.0f;
   bool rvec = true;
-  float sigma = 0.0f;
+  std::complex<float> sigma(0.0f, 0.0f);
 
   std::vector<std::complex<float>> resid(N);
   std::vector<std::complex<float>> V(ncv * N);
   std::vector<std::complex<float>> workd(3 * N);
   std::vector<std::complex<float>> d(nev + 1);
   std::vector<std::complex<float>> z((N + 1) * (nev + 1));
-  std::vector<int> select(ncv);
+  std::vector<a_int> select(ncv);
 
-  BLASINT lworkl = 3 * (ncv * ncv) + 6 * ncv;
+  a_int lworkl = 3 * (ncv * ncv) + 6 * ncv;
   std::vector<std::complex<float>> workl(lworkl);
 
   std::vector<std::complex<float>> rwork(ncv);
   std::vector<std::complex<float>> workev(2 * ncv);
 
-  BLASINT info = 0;
+  a_int info = 0;
 
   int rank;
   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
 
-  std::array<BLASINT, 11> iparam;
+  std::array<a_int, 11> iparam;
   iparam[0] = 1;
   iparam[2] = 10 * N;
   iparam[3] = 1;
   iparam[4] = 0;  // number of ev found by arpack.
   iparam[6] = 1;
 
-  std::array<BLASINT, 14> ipntr;
+  std::array<a_int, 14> ipntr;
 
   MPI_Fint MCW = MPI_Comm_c2f(MPI_COMM_WORLD);
 
-  BLASINT ido = 0;
+  a_int ido = 0;
 
   while (ido != 99) {
     arpack::naupd(MCW, ido, arpack::bmat::identity, N,

PARPACK/TESTS/MPI/issue46.f

@@ -3,7 +3,15 @@ c----------- Example to reproduce issue#46 -------------------------------------
 c
       program issue46
       include 'mpif.h'
-      integer ierr, comm, color, key, myid, nprocs, cnprocs
+c     %-------------------------------%
+c     | MPI INTERFACE                 |
+c     | ILP64 is not supported by MPI |
+c     | integer*4 must be imposed in  |
+c     | all calls involving MPI.      |
+c     |                               |
+c     | Use ierr for MPI calls.       |
+c     %-------------------------------%
+      integer*4 ierr, comm, myid, nprocs, cnprocs, color
       call MPI_INIT( ierr )
       call MPI_COMM_RANK( MPI_COMM_WORLD, myid, ierr )
       call MPI_COMM_SIZE( MPI_COMM_WORLD, nprocs, ierr )
@@ -18,7 +26,8 @@ c     Create commuticator to run arnoldi only on the first CPU
       call MPI_COMM_SIZE( comm, cnprocs, ierr )
       if(color .eq. 1) call parnoldi(comm)
       call MPI_Barrier(MPI_COMM_WORLD, ierr)
-      call parnoldi(MPI_COMM_WORLD)
+      comm = MPI_COMM_WORLD
+      call parnoldi(comm)
       call MPI_FINALIZE(ierr)
       end
 
@@ -28,7 +37,7 @@ c     Create commuticator to run arnoldi only on the first CPU
       include 'debug.h'
       include 'stat.h'
 
-      integer           comm, myid, nprocs, rc, nloc
+      integer*4         comm, myid, nprocs, rc, ierr
 c
 c     %-----------------------------%
 c     | Define leading dimensions   |
@@ -54,12 +63,15 @@ c
       logical          select(maxncv)
       integer          iparam(11), ipntr(11)
 c
-c     %---------------%
-c     | Local Scalars |
-c     %---------------%
+c     %------------------------------------%
+c     | Local Scalars                      |
+c     |                                    |
+c     | Use info if ILP64 can be supported |
+c     | (call to BLAS, LAPACK, ARPACK).    |
+c     %------------------------------------%
 c
       character        bmat*1, which*2
-      integer          ido, n, nev, ncv, lworkl, info, ierr, j,
+      integer          ido, n, nev, ncv, lworkl, info, nloc, j,
      &                 nx, nconv, maxitr, mode, ishfts
       logical          rvec
       Double precision
@@ -255,10 +267,10 @@ c        %-------------------------------------------%
 c
          rvec = .true.
 c
-         call pdseupd ( comm, rvec, 'All', select,
+         call pdseupd ( comm, rvec, 'A', select,
      &        d, v, ldv, sigma,
      &        bmat, nloc, which, nev, tol, resid, ncv, v, ldv,
-     &        iparam, ipntr, workd, workl, lworkl, ierr )
+     &        iparam, ipntr, workd, workl, lworkl, info )
 c        %----------------------------------------------%
 c        | Eigenvalues are returned in the first column |
 c        | of the two dimensional array D and the       |
@@ -270,7 +282,7 @@ c        | corresponding to the eigenvalues in D is     |
 c        | returned in V.                               |
 c        %----------------------------------------------%
 c
-         if ( ierr .ne. 0) then
+         if ( info .ne. 0) then
 c
 c            %------------------------------------%
 c            | Error condition:                   |
@@ -280,7 +292,7 @@ c
 c
             if ( myid .eq. 0 ) then
              	print *, ' '
-             	print *, ' Error with _seupd, info = ', ierr
+             	print *, ' Error with _seupd, info = ', info
              	print *, ' Check the documentation of _seupd. '
              	print *, ' '
             endif
@@ -385,7 +397,7 @@ c
 c
 c     .. MPI Declarations ...
       include           'mpif.h'
-      integer           comm, nprocs, myid, ierr,
+      integer*4         comm, nprocs, myid, ierr,
      &                  status(MPI_STATUS_SIZE)
       integer           nloc, nx, np, j, lo, next, prev
       Double precision

README.md

@@ -15,23 +15,39 @@ Important Features:
 * Routines for The Singular Value Decomposition.
 * Example driver routines that may be used as templates to implement numerous
   Shift-Invert strategies for all problem types, data types and precision.
+* arpackmm: utility to test arpack with matrix market files.
+  Note: to run this utility, you need the eigen library (to handle RCI).
+* ILP64 support:
+  * users: set INTERFACE64 at configure time.
+  * developers:
+    * all files which needs ILP64 support must include "arpackdef.h".
+    * when coding, use a_int (defined in arpackdef.h) instead of int.
+      a_int stands for "architecture int": it's #defined to int or int64_t according
+      to the architecture.
+  * example: to test arpack with sequential ILP64 MKL assuming you use gnu compilers
+    ```$ ./bootstrap
+    $ export FFLAGS='-I/usr/include/mkl'
+    $ export FCFLAGS='-I/usr/include/mkl'
+    $ export LIBS='-Wl,--no-as-needed -lmkl_sequential -lmkl_core -lpthread -lm -ldl'
+    $ export INTERFACE64=1
+    $ ./configure --with-blas=mkl_gf_ilp64 --with-lapack=mkl_gf_ilp64
+    $ make all check```
 
-This project is a joint project between Debian, Octave and Scilab in order to
+This project started as a joint project between Debian, Octave and Scilab in order to
 provide a common and maintained version of arpack.
+This is now a community project maintained by a few volunteers.
 
 Indeed, no single release has been published by Rice university for the last
 few years and since many software (Octave, Scilab, R, Matlab...) forked it and
 implemented their own modifications, arpack-ng aims to tackle this by providing
-a common repository and maintained versions.
+a common repository, maintained versions with a testsuite.
 
 arpack-ng is replacing arpack almost everywhere.
 
-
 1. You have successfully unbundled ARPACK-NG and are now in the ARPACK-NG
    directory that was created for you.
 
-2.
-   The directory SRC contains the top level routines including
+2. The directory SRC contains the top level routines including
    the highest level reverse communication interface routines
 
 * ssaupd, dsaupd - symmetric single and double precision
@@ -70,9 +86,9 @@ arpack-ng is replacing arpack almost everywhere.
 
    The following instructions explain how to make the ARPACK library.
 
-4. Unlike ARPACK, ARPACK-NG is providing autotools based build system
-   and iso_c_binding support (which enables to call fortran subroutines
-   natively from C or C++).
+4. Unlike ARPACK, ARPACK-NG is providing autotools and cmake based build
+   system and iso_c_binding support (which enables to call fortran 
+   subroutines natively from C or C++).
 
 Therefore, the classical commands should work as expected:
    
@@ -120,13 +136,5 @@ To get iso_c_binding support:
    for templates on how to invoke the computational modes of ARPACK.
    Also look in the README.MD file for explanations concerning the
    other documents.
-
-Authors:
-* Danny Sorensen   at  sorensen@caam.rice.edu
-* Richard Lehoucq  at  rblehou@sandia.gov
-* Chao Yang        at  cyang@lbl.gov
-* Kristi Maschhoff at  kristyn@tera.com
-* Sylvestre Ledru  at  sylvestre@debian.org (current maintainer)
-* Allan Cornet
-
+ 
  Good luck and enjoy.

SRC/Makefile.am

@@ -1,4 +1,4 @@
-AUTOMAKE_OPTIONS = subdir-objects # Needed as debug_init/icb.f90 are not in current directory.
+AUTOMAKE_OPTIONS = subdir-objects # Needed as debug/stat* are not in current directory.
 
 SSRC = snaitr.f snapps.f snaup2.f snaupd.f snconv.f sneigh.f sneupd.f sngets.f sstatn.f \
        ssaitr.f ssapps.f ssaup2.f ssaupd.f ssconv.f sseigt.f sseupd.f ssgets.f sstats.f \
@@ -15,20 +15,21 @@ ZSRC = znaitr.f znapps.f znaup2.f znaupd.f zneigh.f zneupd.f zngets.f zstatn.f \
        zgetv0.f zsortc.f
 
 if ICB
-SSRC += icbass.f90 icbasn.f90
-DSRC += icbads.f90 icbadn.f90
-CSRC += icbacn.f90
-ZSRC += icbazn.f90
+SSRC += icbass.F90 icbasn.F90
+DSRC += icbads.F90 icbadn.F90
+CSRC += icbacn.F90
+ZSRC += icbazn.F90
 endif
 
 EXTRA_DIST = debug.h stat.h version.h
 
 lib_LTLIBRARIES = libarpack@LIBSUFFIX@.la
 libarpack@LIBSUFFIX@_la_SOURCES = $(SSRC) $(DSRC) $(CSRC) $(ZSRC)
-libarpack@LIBSUFFIX@_la_SOURCES += $(top_builddir)/debug_init.f90
-if ICB
-libarpack@LIBSUFFIX@_la_SOURCES += $(top_builddir)/debug_icb.f90
-libarpack@LIBSUFFIX@_la_SOURCES += $(top_builddir)/stat_icb.f90
-endif
+libarpack@LIBSUFFIX@_la_SOURCES += $(top_builddir)/dbgini.f
+libarpack@LIBSUFFIX@_la_SOURCES += $(top_builddir)/staini.f
 libarpack@LIBSUFFIX@_la_LIBADD = $(top_builddir)/UTIL/libarpackutil.la $(LAPACK_LIBS) $(BLAS_LIBS)
-libarpack@LIBSUFFIX@_la_LDFLAGS = -no-undefined -version-info 2:0
+libarpack@LIBSUFFIX@_la_LDFLAGS = -no-undefined -version-info 3:0:1
+if ICB
+libarpack@LIBSUFFIX@_la_LIBADD += $(top_builddir)/ICB/libdbgicb.la $(top_builddir)/ICB/libstaicb.la
+libarpack@LIBSUFFIX@_la_CPPFLAGS = $(AM_CPPFLAGS) -I$(top_builddir)
+endif

SRC/cgetv0.f

@@ -239,12 +239,16 @@ c        | the generalized problem when B is possibly (singular).   |
 c        %----------------------------------------------------------%
 c
          call arscnd (t2)
-         nopx = nopx + 1
-         ipntr(1) = 1
-         ipntr(2) = n + 1
-         call ccopy (n, resid, 1, workd, 1)
-         ido = -1
-         go to 9000
+         if (itry .eq. 1) then
+            nopx = nopx + 1
+            ipntr(1) = 1
+            ipntr(2) = n + 1
+            call ccopy (n, resid, 1, workd, 1)
+            ido = -1
+            go to 9000
+         else if (itry .gt. 1 .and. bmat .eq. 'G') then
+            call ccopy (n, resid, 1, workd(n + 1), 1)
+         end if
       end if
 c
 c     %----------------------------------------%
@@ -269,7 +273,7 @@ c     %------------------------------------------------------%
 c
       call arscnd (t2)
       first = .TRUE.
-      call ccopy (n, workd(n+1), 1, resid, 1)
+      if (itry .eq. 1) call ccopy (n, workd(n + 1), 1, resid, 1)
       if (bmat .eq. 'G') then
          nbx = nbx + 1
          ipntr(1) = n + 1

SRC/dgetv0.f

@@ -242,12 +242,16 @@ c        | the generalized problem when B is possibly (singular).   |
 c        %----------------------------------------------------------%
 c
          call arscnd (t2)
-         nopx = nopx + 1
-         ipntr(1) = 1
-         ipntr(2) = n + 1
-         call dcopy (n, resid, 1, workd, 1)
-         ido = -1
-         go to 9000
+         if (itry .eq. 1) then
+            nopx = nopx + 1
+            ipntr(1) = 1
+            ipntr(2) = n + 1
+            call dcopy (n, resid, 1, workd, 1)
+            ido = -1
+            go to 9000
+         else if (itry .gt. 1 .and. bmat .eq. 'G') then
+            call dcopy (n, resid, 1, workd(n + 1), 1)
+         end if
       end if
 c
 c     %-----------------------------------------%
@@ -274,7 +278,7 @@ c     %------------------------------------------------------%
 c
       call arscnd (t2)
       first = .TRUE.
-      call dcopy (n, workd(n+1), 1, resid, 1)
+      if (itry .eq. 1) call dcopy (n, workd(n + 1), 1, resid, 1)
       if (bmat .eq. 'G') then
          nbx = nbx + 1
          ipntr(1) = n + 1

SRC/icbacn.F90

@@ -5,6 +5,7 @@ subroutine cnaupd_c(ido, bmat, n, which, nev, tol, resid, ncv, v, ldv,&
                     bind(c, name="cnaupd_c")
   use :: iso_c_binding
   implicit none
+#include "arpackdef.h"
   integer(kind=c_int),                              intent(inout) :: ido
   character(kind=c_char),       dimension(1),       intent(in)    :: bmat
   integer(kind=c_int),          value,              intent(in)    :: n
@@ -32,6 +33,7 @@ subroutine cneupd_c(rvec, howmny, select, d, z, ldz, sigma, workev,  &
                     bind(c, name="cneupd_c")
   use :: iso_c_binding
   implicit none
+#include "arpackdef.h"
   logical(kind=c_bool),         value,              intent(in)    :: rvec
   character(kind=c_char),       dimension(1),       intent(in)    :: howmny
   logical(kind=c_bool),         dimension(ncv),     intent(in)    :: select

SRC/icbadn.F90

@@ -5,6 +5,7 @@ subroutine dnaupd_c(ido, bmat, n, which, nev, tol, resid, ncv, v, ldv,&
                     bind(c, name="dnaupd_c")
   use :: iso_c_binding
   implicit none
+#include "arpackdef.h"
   integer(kind=c_int),                         intent(inout) :: ido
   character(kind=c_char), dimension(1),        intent(in)    :: bmat
   integer(kind=c_int),    value,               intent(in)    :: n
@@ -26,12 +27,13 @@ subroutine dnaupd_c(ido, bmat, n, which, nev, tol, resid, ncv, v, ldv,&
 end subroutine dnaupd_c
 
 subroutine dneupd_c(rvec, howmny, select,                        &
-                    dr, di, z, ldz, sigmar, sigmai,              &
+                    dr, di, z, ldz, sigmar, sigmai, workev,      &
                     bmat, n, which, nev, tol, resid, ncv, v, ldv,&
                     iparam, ipntr, workd, workl, lworkl, info)   &
                     bind(c, name="dneupd_c")
   use :: iso_c_binding
   implicit none
+#include "arpackdef.h"
   logical(kind=c_bool),   value,               intent(in)    :: rvec
   character(kind=c_char), dimension(1),        intent(in)    :: howmny
   logical(kind=c_bool),   dimension(ncv),      intent(in)    :: select
@@ -41,6 +43,7 @@ subroutine dneupd_c(rvec, howmny, select,                        &
   integer(kind=c_int),    value,               intent(in)    :: ldz
   real(kind=c_double),    value,               intent(in)    :: sigmar
   real(kind=c_double),    value,               intent(in)    :: sigmai
+  real(kind=c_double),    dimension(3*ncv),    intent(out)   :: workev
   character(kind=c_char), dimension(1),        intent(in)    :: bmat
   integer(kind=c_int),    value,               intent(in)    :: n
   character(kind=c_char), dimension(2),        intent(in)    :: which
@@ -57,7 +60,7 @@ subroutine dneupd_c(rvec, howmny, select,                        &
   integer(kind=c_int),    value,               intent(in)    :: lworkl
   integer(kind=c_int),                         intent(inout) :: info
   call dneupd(rvec, howmny, select,                        &
-              dr, di, z, ldz, sigmar, sigmai,              &
+              dr, di, z, ldz, sigmar, sigmai, workev,      &
               bmat, n, which, nev, tol, resid, ncv, v, ldv,&
               iparam, ipntr, workd, workl, lworkl, info)
 end subroutine dneupd_c

SRC/icbads.F90

@@ -5,6 +5,7 @@ subroutine dsaupd_c(ido, bmat, n, which, nev, tol, resid, ncv, v, ldv,&
                     bind(c, name="dsaupd_c")
   use :: iso_c_binding
   implicit none
+#include "arpackdef.h"
   integer(kind=c_int),                         intent(inout) :: ido
   character(kind=c_char), dimension(1),        intent(in)    :: bmat
   integer(kind=c_int),    value,               intent(in)    :: n
@@ -31,6 +32,7 @@ subroutine dseupd_c(rvec, howmny, select, d, z, ldz, sigma,      &
                     bind(c, name="dseupd_c")
   use :: iso_c_binding
   implicit none
+#include "arpackdef.h"
   logical(kind=c_bool),   value,               intent(in)    :: rvec
   character(kind=c_char), dimension(1),        intent(in)    :: howmny
   logical(kind=c_bool),   dimension(ncv),      intent(in)    :: select

SRC/icbasn.F90

@@ -5,6 +5,7 @@ subroutine snaupd_c(ido, bmat, n, which, nev, tol, resid, ncv, v, ldv,&
                     bind(c, name="snaupd_c")
   use :: iso_c_binding
   implicit none
+#include "arpackdef.h"
   integer(kind=c_int),                         intent(inout) :: ido
   character(kind=c_char), dimension(1),        intent(in)    :: bmat
   integer(kind=c_int),    value,               intent(in)    :: n
@@ -26,12 +27,13 @@ subroutine snaupd_c(ido, bmat, n, which, nev, tol, resid, ncv, v, ldv,&
 end subroutine snaupd_c
 
 subroutine sneupd_c(rvec, howmny, select,                        &
-                    dr, di, z, ldz, sigmar, sigmai,              &
+                    dr, di, z, ldz, sigmar, sigmai, workev,      &
                     bmat, n, which, nev, tol, resid, ncv, v, ldv,&
                     iparam, ipntr, workd, workl, lworkl, info)   &
                     bind(c, name="sneupd_c")
   use :: iso_c_binding
   implicit none
+#include "arpackdef.h"
   logical(kind=c_bool),   value,               intent(in)    :: rvec
   character(kind=c_char), dimension(1),        intent(in)    :: howmny
   logical(kind=c_bool),   dimension(ncv),      intent(in)    :: select
@@ -41,6 +43,7 @@ subroutine sneupd_c(rvec, howmny, select,                        &
   integer(kind=c_int),    value,               intent(in)    :: ldz
   real(kind=c_float),     value,               intent(in)    :: sigmar
   real(kind=c_float),     value,               intent(in)    :: sigmai
+  real(kind=c_float),     dimension(3*ncv),    intent(out)   :: workev
   character(kind=c_char), dimension(1),        intent(in)    :: bmat
   integer(kind=c_int),    value,               intent(in)    :: n
   character(kind=c_char), dimension(2),        intent(in)    :: which
@@ -57,7 +60,7 @@ subroutine sneupd_c(rvec, howmny, select,                        &
   integer(kind=c_int),    value,               intent(in)    :: lworkl
   integer(kind=c_int),                         intent(inout) :: info
   call sneupd(rvec, howmny, select,                        &
-              dr, di, z, ldz, sigmar, sigmai,              &
+              dr, di, z, ldz, sigmar, sigmai, workev,      &
               bmat, n, which, nev, tol, resid, ncv, v, ldv,&
               iparam, ipntr, workd, workl, lworkl, info)
 end subroutine sneupd_c

SRC/icbass.F90

@@ -5,6 +5,7 @@ subroutine ssaupd_c(ido, bmat, n, which, nev, tol, resid, ncv, v, ldv,&
                     bind(c, name="ssaupd_c")
   use :: iso_c_binding
   implicit none
+#include "arpackdef.h"
   integer(kind=c_int),                         intent(inout) :: ido
   character(kind=c_char), dimension(1),        intent(in)    :: bmat
   integer(kind=c_int),    value,               intent(in)    :: n
@@ -31,6 +32,7 @@ subroutine sseupd_c(rvec, howmny, select, d, z, ldz, sigma,      &
                     bind(c, name="sseupd_c")
   use :: iso_c_binding
   implicit none
+#include "arpackdef.h"
   logical(kind=c_bool),   value,               intent(in)    :: rvec
   character(kind=c_char), dimension(1),        intent(in)    :: howmny
   logical(kind=c_bool),   dimension(ncv),      intent(in)    :: select

SRC/icbazn.F90

@@ -5,6 +5,7 @@ subroutine znaupd_c(ido, bmat, n, which, nev, tol, resid, ncv, v, ldv,&
                     bind(c, name="znaupd_c")
   use :: iso_c_binding
   implicit none
+#include "arpackdef.h"
   integer(kind=c_int),                                  intent(inout) :: ido
   character(kind=c_char),           dimension(1),       intent(in)    :: bmat
   integer(kind=c_int),              value,              intent(in)    :: n
@@ -32,6 +33,7 @@ subroutine zneupd_c(rvec, howmny, select, d, z, ldz, sigma, workev,  &
                     bind(c, name="zneupd_c")
   use :: iso_c_binding
   implicit none
+#include "arpackdef.h"
   logical(kind=c_bool),             value,              intent(in)    :: rvec
   character(kind=c_char),           dimension(1),       intent(in)    :: howmny
   logical(kind=c_bool),             dimension(ncv),     intent(in)    :: select

SRC/sgetv0.f

@@ -242,12 +242,16 @@ c        | the generalized problem when B is possibly (singular).   |
 c        %----------------------------------------------------------%
 c
          call arscnd (t2)
-         nopx = nopx + 1
-         ipntr(1) = 1
-         ipntr(2) = n + 1
-         call scopy (n, resid, 1, workd, 1)
-         ido = -1
-         go to 9000
+         if (itry .eq. 1) then
+            nopx = nopx + 1
+            ipntr(1) = 1
+            ipntr(2) = n + 1
+            call scopy (n, resid, 1, workd, 1)
+            ido = -1
+            go to 9000
+         else if (itry .gt. 1 .and. bmat .eq. 'G') then
+            call scopy (n, resid, 1, workd(n + 1), 1)
+         end if
       end if
 c
 c     %-----------------------------------------%
@@ -274,7 +278,7 @@ c     %------------------------------------------------------%
 c
       call arscnd (t2)
       first = .TRUE.
-      call scopy (n, workd(n+1), 1, resid, 1)
+      if (itry .eq. 1) call scopy (n, workd(n + 1), 1, resid, 1)
       if (bmat .eq. 'G') then
          nbx = nbx + 1
          ipntr(1) = n + 1

SRC/zgetv0.f

@@ -239,12 +239,16 @@ c        | the generalized problem when B is possibly (singular).   |
 c        %----------------------------------------------------------%
 c
          call arscnd (t2)
-         nopx = nopx + 1
-         ipntr(1) = 1
-         ipntr(2) = n + 1
-         call zcopy (n, resid, 1, workd, 1)
-         ido = -1
-         go to 9000
+         if (itry .eq. 1) then
+            nopx = nopx + 1
+            ipntr(1) = 1
+            ipntr(2) = n + 1
+            call zcopy (n, resid, 1, workd, 1)
+            ido = -1
+            go to 9000
+         else if (itry .gt. 1 .and. bmat .eq. 'G') then
+            call zcopy (n, resid, 1, workd(n + 1), 1)
+         end if
       end if
 c
 c     %----------------------------------------%
@@ -269,7 +273,7 @@ c     %------------------------------------------------------%
 c
       call arscnd (t2)
       first = .TRUE.
-      call zcopy (n, workd(n+1), 1, resid, 1)
+      if (itry .eq. 1) call zcopy (n, workd(n + 1), 1, resid, 1)
       if (bmat .eq. 'G') then
          nbx = nbx + 1
          ipntr(1) = n + 1

TESTS/Makefile.am

@@ -1,6 +1,6 @@
 LDADD=$(top_builddir)/SRC/libarpack$(LIBSUFFIX).la $(LAPACK_LIBS) $(BLAS_LIBS)
 
-check_PROGRAMS = dnsimp bug_1323 bug_58_double bug_79_double_complex
+check_PROGRAMS = dnsimp bug_1323 bug_58_double bug_79_double_complex bug_142 bug_142_gen
 if ICB
 check_PROGRAMS += icb_arpack_c
 check_PROGRAMS += icb_arpack_cpp
@@ -8,11 +8,9 @@ else
 check_PROGRAMS += bug_1315_single
 check_PROGRAMS += bug_1315_double
 endif
-
+dist_check_DATA = testA.mtx
 TESTS = $(check_PROGRAMS)
 
-EXTRA_DIST = testA.mtx
-
 dnsimp_SOURCES = dnsimp.f mmio.f debug.h
 
 bug_1323_SOURCES = bug_1323.f
@@ -21,17 +19,23 @@ bug_58_double_SOURCES = bug_58_double.f
 
 bug_79_double_complex_SOURCES = bug_79_double_complex.f
 
+bug_142_SOURCES = bug_142.f
+
+bug_142_gen_SOURCES = bug_142_gen.f
+
 if ICB
 icb_arpack_c_SOURCES = icb_arpack_c.c
 icb_arpack_c_LDADD = $(top_builddir)/SRC/libarpack$(LIBSUFFIX).la
-icb_arpack_c_CFLAGS = -I$(top_builddir)
+icb_arpack_c_CPPFLAGS = $(AM_CPPFLAGS) -I$(top_builddir) -I$(top_srcdir)/ICB
 
 icb_arpack_cpp_SOURCES = icb_arpack_cpp.cpp
 icb_arpack_cpp_LDADD = $(top_builddir)/SRC/libarpack$(LIBSUFFIX).la
-icb_arpack_cpp_CPPFLAGS = -I$(top_builddir)
+icb_arpack_cpp_CPPFLAGS = $(AM_CPPFLAGS) -I$(top_builddir) -I$(top_srcdir)/ICB
 else
 bug_1315_single_SOURCES = bug_1315_single.c
 bug_1315_double_SOURCES = bug_1315_double.c
 bug_1315_single_LDADD = $(top_builddir)/SRC/libarpack$(LIBSUFFIX).la
 bug_1315_double_LDADD = $(top_builddir)/SRC/libarpack$(LIBSUFFIX).la
+bug_1315_single_CPPFLAGS = $(AM_CPPFLAGS) -I$(top_builddir)
+bug_1315_double_CPPFLAGS = $(AM_CPPFLAGS) -I$(top_builddir)
 endif

TESTS/bug_1315_double.c

@@ -1,3 +1,5 @@
+#include "arpackdef.h"
+
 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>
@@ -16,18 +18,15 @@
  * This is not efficient since the problem is
  * symmetric but is done to exhibit the bug.
  * */
-#ifndef BLASINT
-#define BLASINT int
-#endif
 
-extern void dnaupd(BLASINT *, char *, BLASINT *, char *, BLASINT *,
-		   double *, double *, BLASINT *, double *,
-		   BLASINT *, BLASINT *, BLASINT *, double *,
-		   double *, BLASINT *, BLASINT *);
+extern void dnaupd(a_int *, char *, a_int *, char *, a_int *,
+		   double *, double *, a_int *, double *,
+		   a_int *, a_int *, a_int *, double *,
+		   double *, a_int *, a_int *);
 
-extern void dneupd( BLASINT*, char*, BLASINT *, double *, double *, double *, BLASINT*, double *,
-		double *, double *, char *, BLASINT *, char *, BLASINT *, double *, double *, BLASINT *,
-		double *, BLASINT *, BLASINT *, BLASINT *, double *, double *, BLASINT *, BLASINT * );
+extern void dneupd( a_int*, char*, a_int *, double *, double *, double *, a_int*, double *,
+		double *, double *, char *, a_int *, char *, a_int *, double *, double *, a_int *,
+		double *, a_int *, a_int *, a_int *, double *, double *, a_int *, a_int * );
 
 void matVec(double * x, double * y) {
   int i;
@@ -36,26 +35,26 @@ void matVec(double * x, double * y) {
 };
 
 int main() {
-  BLASINT ido = 0;
+  a_int ido = 0;
   char bmat[] = "I";
-  BLASINT N = 1000;
+  a_int N = 1000;
   char which[] = "LM";
-  BLASINT nev = 9;
+  a_int nev = 9;
   double tol = 0;
   double resid[N];
-  BLASINT ncv = 2*nev+1;
+  a_int ncv = 2*nev+1;
   double V[ncv*N];
-  BLASINT ldv = N;
-  BLASINT iparam[11];
-  BLASINT ipntr[14];
+  a_int ldv = N;
+  a_int iparam[11];
+  a_int ipntr[14];
   double workd[3*N];
-  BLASINT rvec = 1;
+  a_int rvec = 1;
   char howmny[] = "A";
   double* dr = (double*) malloc((nev+1)*sizeof(double));
   double* di = (double*) malloc((nev+1)*sizeof(double));
-  BLASINT select[3*ncv];
+  a_int select[3*ncv];
   double z[(N+1)*(nev+1)];
-  BLASINT ldz = N+1;
+  a_int ldz = N+1;
   double sigmar=0;
   double sigmai=0;
   double workev[3*ncv];
@@ -65,8 +64,8 @@ int main() {
   double workl[3*(ncv*ncv) + 6*ncv];
   for (k=0; k < 3*(ncv*ncv) + 6*ncv; ++k )
     workl[k] = 0;
-  BLASINT lworkl = 3*(ncv*ncv) + 6*ncv;
-  BLASINT info = 0;
+  a_int lworkl = 3*(ncv*ncv) + 6*ncv;
+  a_int info = 0;
 
   iparam[0] = 1;
   iparam[2] = 10*N;

TESTS/bug_1315_single.c

@@ -1,3 +1,5 @@
+#include "arpackdef.h"
+
 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>
@@ -16,18 +18,15 @@
  * This is not efficient since the problem is
  * symmetric but is done to exhibit the bug.
  */
-#ifndef BLASINT
-#define BLASINT int
-#endif
 
-extern void snaupd(BLASINT *, char *, BLASINT *, char *, BLASINT *,
-		   float *, float *, BLASINT *, float *,
-		   BLASINT *, BLASINT *, BLASINT *, float *,
-		   float *, BLASINT *, BLASINT *);
+extern void snaupd(a_int *, char *, a_int *, char *, a_int *,
+		   float *, float *, a_int *, float *,
+		   a_int *, a_int *, a_int *, float *,
+		   float *, a_int *, a_int *);
 
-extern void sneupd( BLASINT*, char*, BLASINT *, float *, float *, float *, BLASINT*, float *,
-		float *, float *, char *, BLASINT *, char *, BLASINT *, float *, float *, BLASINT *,
-		float *, BLASINT *, BLASINT *, BLASINT *, float *, float *, BLASINT *, BLASINT * );
+extern void sneupd( a_int*, char*, a_int *, float *, float *, float *, a_int*, float *,
+		float *, float *, char *, a_int *, char *, a_int *, float *, float *, a_int *,
+		float *, a_int *, a_int *, a_int *, float *, float *, a_int *, a_int * );
 
 
 void matVec(float * x, float * y) {
@@ -37,26 +36,26 @@ void matVec(float * x, float * y) {
 };
 
 int main() {
-  BLASINT ido = 0;
+  a_int ido = 0;
   char bmat[] = "I";
-  BLASINT N = 1000;
+  a_int N = 1000;
   char which[] = "LM";
-  BLASINT nev = 9;
+  a_int nev = 9;
   float tol = 0;
   float resid[N];
-  BLASINT ncv = 2*nev+1;
+  a_int ncv = 2*nev+1;
   float V[ncv*N];
-  BLASINT ldv = N;
-  BLASINT iparam[11];
-  BLASINT ipntr[14];
+  a_int ldv = N;
+  a_int iparam[11];
+  a_int ipntr[14];
   float workd[3*N];
-  BLASINT rvec = 1;
+  a_int rvec = 1;
   char howmny[] = "A";
   float* dr = (float*) malloc((nev+1)*sizeof(float));
   float* di = (float*) malloc((nev+1)*sizeof(float));
-  BLASINT select[3*ncv];
+  a_int select[3*ncv];
   float z[(N+1)*(nev+1)];
-  BLASINT ldz = N+1;
+  a_int ldz = N+1;
   float sigmar=0;
   float sigmai=0;
   float workev[3*ncv];
@@ -66,8 +65,8 @@ int main() {
   float workl[3*(ncv*ncv) + 6*ncv];
   for (k=0; k < 3*(ncv*ncv) + 6*ncv; ++k )
     workl[k] = 0;
-  BLASINT lworkl = 3*(ncv*ncv) + 6*ncv;
-  BLASINT info = 0;
+  a_int lworkl = 3*(ncv*ncv) + 6*ncv;
+  a_int info = 0;
 
   iparam[0] = 1;
   iparam[2] = 10*N;

TESTS/bug_142.f

@@ -0,0 +1,414 @@
+      program bug_142
+c
+c     Avoid taking the initial vector in the range of OP after a restart
+c     (standard case)
+c
+c     Example program to illustrate the idea of reverse communication
+c     for a standard nonsymmetric eigenvalue problem.
+c
+c     We implement example one of ex-nonsym.doc in DOCUMENTS directory
+c
+c\Example-1
+c     ... Suppose we want to solve A*x = lambda*x in regular mode,
+c         where A is obtained from the standard central difference
+c         discretization of the convection-diffusion operator 
+c                 (Laplacian u) + rho*(du / dx)
+c         on the unit square [0,1]x[0,1] with zero Dirichlet boundary 
+c         condition.
+c
+c     ... OP = A  and  B = I.
+c
+c     ... Assume "call av (nx,x,y)" computes y = A*x.c
+c
+c     ... Use mode 1 of DNAUPD.
+c
+c\BeginLib
+c
+c\Routines called:
+c     dnaupd  ARPACK reverse communication interface routine.
+c     dneupd  ARPACK routine that returns Ritz values and (optionally)
+c             Ritz vectors.
+c     dlapy2  LAPACK routine to compute sqrt(x**2+y**2) carefully.
+c     daxpy   Level 1 BLAS that computes y <- alpha*x+y.
+c     dnrm2   Level 1 BLAS that computes the norm of a vector.
+c     av      Matrix vector multiplication routine that computes A*x.
+c     tv      Matrix vector multiplication routine that computes T*x, 
+c             where T is a tridiagonal matrix.  It is used in routine
+c             av.
+c
+c\Author
+c     Richard Lehoucq
+c     Danny Sorensen
+c     Chao Yang
+c     Dept. of Computational &
+c     Applied Mathematics
+c     Rice University
+c     Houston, Texas
+c
+c\SCCS Information: @(#)
+c FILE: ndrv1.F   SID: 2.5   DATE OF SID: 10/17/00   RELEASE: 2
+c
+c\Remarks
+c     1. None
+c
+c\EndLib
+c---------------------------------------------------------------------------
+c
+c     %-----------------------------%
+c     | Define maximum dimensions   |
+c     | for all arrays.             |
+c     | MAXN:   Maximum dimension   |
+c     |         of the A allowed.   |
+c     | MAXNEV: Maximum NEV allowed |
+c     | MAXNCV: Maximum NCV allowed |
+c     %-----------------------------%
+c
+      integer           maxn, maxnev, maxncv, ldv
+      parameter         (maxn=256, maxnev=12, maxncv=30, ldv=maxn)
+c
+c     %--------------%
+c     | Local Arrays |
+c     %--------------%
+c
+      integer           iparam(11), ipntr(14)
+      logical           select(maxncv)
+      Double precision
+     &                  ax(maxn), d(maxncv,3), resid(maxn), 
+     &                  v(ldv,maxncv), workd(3*maxn), 
+     &                  workev(3*maxncv), 
+     &                  workl(3*maxncv*maxncv+6*maxncv), a(maxn, maxn)
+c
+c     %---------------%
+c     | Local Scalars |
+c     %---------------%
+c
+      character         bmat*1, which*2
+      integer           ido, n, nx, nev, ncv, lworkl, info, j,
+     &                  ierr, nconv, maxitr, ishfts, mode
+      Double precision
+     &                  tol, sigmar, sigmai
+      logical           first, rvec
+c
+c     %------------%
+c     | Parameters |
+c     %------------%
+c
+      Double precision
+     &                  zero
+      parameter         (zero = 0.0D+0)
+c
+c     %-----------------------------%
+c     | BLAS & LAPACK routines used |
+c     %-----------------------------%
+c
+      Double precision
+     &                  dlapy2, dnrm2
+      external          dlapy2, dnrm2, daxpy
+c
+c     %--------------------%
+c     | Intrinsic function |
+c     %--------------------%
+c
+      intrinsic         abs
+c
+c     %-----------------------%
+c     | Executable Statements |
+c     %-----------------------%
+c 
+c     %--------------------------------------------------%
+c     | The number NX is the number of interior points   |
+c     | in the discretization of the 2-dimensional       |
+c     | convection-diffusion operator on the unit        |
+c     | square with zero Dirichlet boundary condition.   | 
+c     | The number N(=NX*NX) is the dimension of the     |
+c     | matrix.  A standard eigenvalue problem is        |
+c     | solved (BMAT = 'I').  NEV is the number of       |
+c     | eigenvalues to be approximated.  The user can    |
+c     | modify NX, NEV, NCV, WHICH to solve problems of  |
+c     | different sizes, and to get different parts of   |
+c     | the spectrum.  However, The following            |
+c     | conditions must be satisfied:                    |
+c     |                   N <= MAXN                      |
+c     |                 NEV <= MAXNEV                    |
+c     |           NEV + 2 <= NCV <= MAXNCV               | 
+c     %--------------------------------------------------% 
+c
+      nx    = 10 
+      n     = 11
+      nev   = 1
+      ncv   = 11 
+      do i = 1,n
+         do j = 1,n
+            a(i,j) = 0.15d0/11
+         end do
+      end do
+      do j = 2,n
+         a(1,j) = a(1,j) + 0.85d0
+      end do
+      do i = 2,n
+         a(i,1) = (1-a(1,1))/10
+      end do
+
+      if ( n .gt. maxn ) then
+         print *, ' ERROR with _NDRV1: N is greater than MAXN '
+         go to 9000
+      else if ( nev .gt. maxnev ) then
+         print *, ' ERROR with _NDRV1: NEV is greater than MAXNEV '
+         go to 9000
+      else if ( ncv .gt. maxncv ) then
+         print *, ' ERROR with _NDRV1: NCV is greater than MAXNCV '
+         go to 9000
+      end if
+      bmat  = 'I'
+      which = 'LM'
+c
+c     %-----------------------------------------------------%
+c     | The work array WORKL is used in DNAUPD as           |  
+c     | workspace.  Its dimension LWORKL is set as          |
+c     | illustrated below.  The parameter TOL determines    |
+c     | the stopping criterion. If TOL<=0, machine          |
+c     | precision is used.  The variable IDO is used for    |
+c     | reverse communication, and is initially set to 0.   |
+c     | Setting INFO=0 indicates that a random vector is    |
+c     | generated in DNAUPD to start the Arnoldi iteration. | 
+c     %-----------------------------------------------------%
+c
+      lworkl  = 3*ncv**2+6*ncv 
+      tol    = zero 
+      ido    = 0
+      info   = 0
+c
+c     %---------------------------------------------------%
+c     | This program uses exact shifts with respect to    |
+c     | the current Hessenberg matrix (IPARAM(1) = 1).    |
+c     | IPARAM(3) specifies the maximum number of Arnoldi |
+c     | iterations allowed.  Mode 1 of DNAUPD is used     |
+c     | (IPARAM(7) = 1). All these options can be changed |
+c     | by the user. For details see the documentation in |
+c     | DNAUPD.                                           |
+c     %---------------------------------------------------%
+c
+      ishfts = 1
+      maxitr = 300
+      mode   = 1
+c
+      iparam(1) = ishfts
+      iparam(3) = maxitr 
+      iparam(7) = mode
+c
+c     %-------------------------------------------%
+c     | M A I N   L O O P (Reverse communication) | 
+c     %-------------------------------------------%
+c
+ 10   continue
+c
+c        %---------------------------------------------%
+c        | Repeatedly call the routine DNAUPD and take |
+c        | actions indicated by parameter IDO until    |
+c        | either convergence is indicated or maxitr   |
+c        | has been exceeded.                          |
+c        %---------------------------------------------%
+c
+         call dnaupd ( ido, bmat, n, which, nev, tol, resid, 
+     &        ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, 
+     &        info )
+c
+         if (ido .eq. -1 .or. ido .eq. 1) then
+c
+c           %-------------------------------------------%
+c           | Perform matrix vector multiplication      |
+c           |                y <--- OP*x                |
+c           | The user should supply his/her own        |
+c           | matrix vector multiplication routine here |
+c           | that takes workd(ipntr(1)) as the input   |
+c           | vector, and return the matrix vector      |
+c           | product to workd(ipntr(2)).               | 
+c           %-------------------------------------------%
+c
+            call dgemv('N',n,n,1.0d0,a,maxn,
+     &           workd(ipntr(1)),1,0.0d0,workd(ipntr(2)),1)
+c
+c           %-----------------------------------------%
+c           | L O O P   B A C K to call DNAUPD again. |
+c           %-----------------------------------------%
+c
+            go to 10
+c
+         else if (ido .eq. 2) then
+
+            call dcopy (n, workd(ipntr(1)), 1, workd(ipntr(2)), 1)
+
+            go to 10
+
+         end if 
+c 
+c     %----------------------------------------%
+c     | Either we have convergence or there is |
+c     | an error.                              |
+c     %----------------------------------------%
+c
+      if ( info .lt. 0 ) then
+c
+c        %--------------------------%
+c        | Error message, check the |
+c        | documentation in DNAUPD. |
+c        %--------------------------%
+c
+         print *, ' '
+         print *, ' Error with _naupd, info = ', info
+         print *, ' Check the documentation of _naupd'
+         print *, ' '
+         stop 1
+c
+      else 
+c
+c        %-------------------------------------------%
+c        | No fatal errors occurred.                 |
+c        | Post-Process using DNEUPD.                |
+c        |                                           |
+c        | Computed eigenvalues may be extracted.    |
+c        |                                           |
+c        | Eigenvectors may also be computed now if  |
+c        | desired.  (indicated by rvec = .true.)    |
+c        %-------------------------------------------%
+c
+         rvec = .true.
+c
+         call dneupd ( rvec, 'A', select, d, d(1,2), v, ldv, 
+     &        sigmar, sigmai, workev, bmat, n, which, nev, tol, 
+     &        resid, ncv, v, ldv, iparam, ipntr, workd, workl,
+     &        lworkl, ierr )
+c
+c        %-----------------------------------------------%
+c        | The real part of the eigenvalue is returned   |
+c        | in the first column of the two dimensional    |
+c        | array D, and the imaginary part is returned   |
+c        | in the second column of D.  The corresponding |
+c        | eigenvectors are returned in the first NEV    |
+c        | columns of the two dimensional array V if     |
+c        | requested.  Otherwise, an orthogonal basis    |
+c        | for the invariant subspace corresponding to   |
+c        | the eigenvalues in D is returned in V.        |
+c        %-----------------------------------------------%
+c
+         if ( ierr .ne. 0) then
+c
+c           %------------------------------------%
+c           | Error condition:                   |
+c           | Check the documentation of DNEUPD. |
+c           %------------------------------------%
+c
+             print *, ' '
+             print *, ' Error with _neupd, info = ', ierr
+             print *, ' Check the documentation of _neupd. '
+             print *, ' '
+c
+         else 
+c
+             first  = .true.
+             nconv  = iparam(5)
+             do 20 j=1, nconv
+c
+c               %---------------------------%
+c               | Compute the residual norm |
+c               |                           |
+c               |   ||  A*x - lambda*x ||   |
+c               |                           |
+c               | for the NCONV accurately  |
+c               | computed eigenvalues and  |
+c               | eigenvectors.  (iparam(5) |
+c               | indicates how many are    |
+c               | accurate to the requested |
+c               | tolerance)                |
+c               %---------------------------%
+c
+                if (d(j,2) .eq. zero)  then
+c
+c                  %--------------------%
+c                  | Ritz value is real |
+c                  %--------------------%
+c
+                   call dgemv('N',n,n,1.0d0,a,maxn,
+     &                  v(1,j),1,0.0d0,ax,1)
+                   call daxpy(n, -d(j,1), v(1,j), 1, ax, 1)
+                   d(j,3) = dnrm2(n, ax, 1)
+                   d(j,3) = d(j,3) / abs(d(j,1))
+c
+                else if (first) then
+c
+c                  %------------------------%
+c                  | Ritz value is complex. |
+c                  | Residual of one Ritz   |
+c                  | value of the conjugate |
+c                  | pair is computed.      | 
+c                  %------------------------%
+c        
+                   call dgemv('N',n,n,1.0d0,a,maxn,
+     &                  v(1,j),1,0.0d0,ax,1)
+                   call daxpy(n, -d(j,1), v(1,j), 1, ax, 1)
+                   call daxpy(n, d(j,2), v(1,j+1), 1, ax, 1)
+                   d(j,3) = dnrm2(n, ax, 1)
+                   call dgemv('N',n,n,1.0d0,a,maxn,
+     &                  v(1,j+1),1,0.0d0,ax,1)
+                   call daxpy(n, -d(j,2), v(1,j), 1, ax, 1)
+                   call daxpy(n, -d(j,1), v(1,j+1), 1, ax, 1)
+                   d(j,3) = dlapy2( d(j,3), dnrm2(n, ax, 1) )
+                   d(j,3) = d(j,3) / dlapy2(d(j,1),d(j,2))
+                   d(j+1,3) = d(j,3)
+                   first = .false.
+                else
+                   first = .true.
+                end if
+c
+ 20          continue
+c
+c            %-----------------------------%
+c            | Display computed residuals. |
+c            %-----------------------------%
+c
+             call dmout(6, nconv, 3, d, maxncv, -6,
+     &            'Ritz values (Real,Imag) and relative residuals')
+          end if
+c
+c        %-------------------------------------------%
+c        | Print additional convergence information. |
+c        %-------------------------------------------%
+c
+         if ( info .eq. 1) then
+             print *, ' '
+             print *, ' Maximum number of iterations reached.'
+             print *, ' '
+         else if ( info .eq. 3) then
+             print *, ' ' 
+             print *, ' No shifts could be applied during implicit',
+     &                ' Arnoldi update, try increasing NCV.'
+             print *, ' '
+         end if      
+c
+         print *, ' '
+         print *, ' _NDRV1 '
+         print *, ' ====== '
+         print *, ' ' 
+         print *, ' Size of the matrix is ', n
+         print *, ' The number of Ritz values requested is ', nev
+         print *, ' The number of Arnoldi vectors generated',
+     &            ' (NCV) is ', ncv
+         print *, ' What portion of the spectrum: ', which
+         print *, ' The number of converged Ritz values is ', 
+     &              nconv 
+         print *, ' The number of Implicit Arnoldi update',
+     &            ' iterations taken is ', iparam(3)
+         print *, ' The number of OP*x is ', iparam(9)
+         print *, ' The convergence criterion is ', tol
+         print *, ' '
+c
+      end if
+c
+c     %---------------------------%
+c     | Done with program dndrv1. |
+c     %---------------------------%
+c
+ 9000 continue
+c
+      end
+c 

TESTS/bug_142_gen.f

@@ -0,0 +1,414 @@
+      program bug_142_gen
+c
+c     Avoid taking the initial vector in the range of OP after a restart
+c     (generalized case)
+c
+c     Example program to illustrate the idea of reverse communication
+c     for a standard nonsymmetric eigenvalue problem.
+c
+c     We implement example one of ex-nonsym.doc in DOCUMENTS directory
+c
+c\Example-1
+c     ... Suppose we want to solve A*x = lambda*x in regular mode,
+c         where A is obtained from the standard central difference
+c         discretization of the convection-diffusion operator 
+c                 (Laplacian u) + rho*(du / dx)
+c         on the unit square [0,1]x[0,1] with zero Dirichlet boundary 
+c         condition.
+c
+c     ... OP = A  and  B = I.
+c
+c     ... Assume "call av (nx,x,y)" computes y = A*x.c
+c
+c     ... Use mode 1 of DNAUPD.
+c
+c\BeginLib
+c
+c\Routines called:
+c     dnaupd  ARPACK reverse communication interface routine.
+c     dneupd  ARPACK routine that returns Ritz values and (optionally)
+c             Ritz vectors.
+c     dlapy2  LAPACK routine to compute sqrt(x**2+y**2) carefully.
+c     daxpy   Level 1 BLAS that computes y <- alpha*x+y.
+c     dnrm2   Level 1 BLAS that computes the norm of a vector.
+c     av      Matrix vector multiplication routine that computes A*x.
+c     tv      Matrix vector multiplication routine that computes T*x, 
+c             where T is a tridiagonal matrix.  It is used in routine
+c             av.
+c
+c\Author
+c     Richard Lehoucq
+c     Danny Sorensen
+c     Chao Yang
+c     Dept. of Computational &
+c     Applied Mathematics
+c     Rice University
+c     Houston, Texas
+c
+c\SCCS Information: @(#)
+c FILE: ndrv1.F   SID: 2.5   DATE OF SID: 10/17/00   RELEASE: 2
+c
+c\Remarks
+c     1. None
+c
+c\EndLib
+c---------------------------------------------------------------------------
+c
+c     %-----------------------------%
+c     | Define maximum dimensions   |
+c     | for all arrays.             |
+c     | MAXN:   Maximum dimension   |
+c     |         of the A allowed.   |
+c     | MAXNEV: Maximum NEV allowed |
+c     | MAXNCV: Maximum NCV allowed |
+c     %-----------------------------%
+c
+      integer           maxn, maxnev, maxncv, ldv
+      parameter         (maxn=256, maxnev=12, maxncv=30, ldv=maxn)
+c
+c     %--------------%
+c     | Local Arrays |
+c     %--------------%
+c
+      integer           iparam(11), ipntr(14)
+      logical           select(maxncv)
+      Double precision
+     &                  ax(maxn), d(maxncv,3), resid(maxn), 
+     &                  v(ldv,maxncv), workd(3*maxn), 
+     &                  workev(3*maxncv), 
+     &                  workl(3*maxncv*maxncv+6*maxncv), a(maxn, maxn)
+c
+c     %---------------%
+c     | Local Scalars |
+c     %---------------%
+c
+      character         bmat*1, which*2
+      integer           ido, n, nx, nev, ncv, lworkl, info, j,
+     &                  ierr, nconv, maxitr, ishfts, mode
+      Double precision
+     &                  tol, sigmar, sigmai
+      logical           first, rvec
+c
+c     %------------%
+c     | Parameters |
+c     %------------%
+c
+      Double precision
+     &                  zero
+      parameter         (zero = 0.0D+0)
+c
+c     %-----------------------------%
+c     | BLAS & LAPACK routines used |
+c     %-----------------------------%
+c
+      Double precision
+     &                  dlapy2, dnrm2
+      external          dlapy2, dnrm2, daxpy
+c
+c     %--------------------%
+c     | Intrinsic function |
+c     %--------------------%
+c
+      intrinsic         abs
+c
+c     %-----------------------%
+c     | Executable Statements |
+c     %-----------------------%
+c 
+c     %--------------------------------------------------%
+c     | The number NX is the number of interior points   |
+c     | in the discretization of the 2-dimensional       |
+c     | convection-diffusion operator on the unit        |
+c     | square with zero Dirichlet boundary condition.   | 
+c     | The number N(=NX*NX) is the dimension of the     |
+c     | matrix.  A standard eigenvalue problem is        |
+c     | solved (BMAT = 'I').  NEV is the number of       |
+c     | eigenvalues to be approximated.  The user can    |
+c     | modify NX, NEV, NCV, WHICH to solve problems of  |
+c     | different sizes, and to get different parts of   |
+c     | the spectrum.  However, The following            |
+c     | conditions must be satisfied:                    |
+c     |                   N <= MAXN                      |
+c     |                 NEV <= MAXNEV                    |
+c     |           NEV + 2 <= NCV <= MAXNCV               | 
+c     %--------------------------------------------------% 
+c
+      nx    = 10 
+      n     = 11
+      nev   = 1
+      ncv   = 11 
+      do i = 1,n
+         do j = 1,n
+            a(i,j) = 0.15d0/11
+         end do
+      end do
+      do j = 2,n
+         a(1,j) = a(1,j) + 0.85d0
+      end do
+      do i = 2,n
+         a(i,1) = (1-a(1,1))/10
+      end do
+
+      if ( n .gt. maxn ) then
+         print *, ' ERROR with _NDRV1: N is greater than MAXN '
+         go to 9000
+      else if ( nev .gt. maxnev ) then
+         print *, ' ERROR with _NDRV1: NEV is greater than MAXNEV '
+         go to 9000
+      else if ( ncv .gt. maxncv ) then
+         print *, ' ERROR with _NDRV1: NCV is greater than MAXNCV '
+         go to 9000
+      end if
+      bmat  = 'G'
+      which = 'LM'
+c
+c     %-----------------------------------------------------%
+c     | The work array WORKL is used in DNAUPD as           |  
+c     | workspace.  Its dimension LWORKL is set as          |
+c     | illustrated below.  The parameter TOL determines    |
+c     | the stopping criterion. If TOL<=0, machine          |
+c     | precision is used.  The variable IDO is used for    |
+c     | reverse communication, and is initially set to 0.   |
+c     | Setting INFO=0 indicates that a random vector is    |
+c     | generated in DNAUPD to start the Arnoldi iteration. | 
+c     %-----------------------------------------------------%
+c
+      lworkl  = 3*ncv**2+6*ncv 
+      tol    = zero 
+      ido    = 0
+      info   = 0
+c
+c     %---------------------------------------------------%
+c     | This program uses exact shifts with respect to    |
+c     | the current Hessenberg matrix (IPARAM(1) = 1).    |
+c     | IPARAM(3) specifies the maximum number of Arnoldi |
+c     | iterations allowed.  Mode 1 of DNAUPD is used     |
+c     | (IPARAM(7) = 1). All these options can be changed |
+c     | by the user. For details see the documentation in |
+c     | DNAUPD.                                           |
+c     %---------------------------------------------------%
+c
+      ishfts = 1
+      maxitr = 300
+      mode   = 2
+c
+      iparam(1) = ishfts
+      iparam(3) = maxitr 
+      iparam(7) = mode
+c
+c     %-------------------------------------------%
+c     | M A I N   L O O P (Reverse communication) | 
+c     %-------------------------------------------%
+c
+ 10   continue
+c
+c        %---------------------------------------------%
+c        | Repeatedly call the routine DNAUPD and take |
+c        | actions indicated by parameter IDO until    |
+c        | either convergence is indicated or maxitr   |
+c        | has been exceeded.                          |
+c        %---------------------------------------------%
+c
+         call dnaupd ( ido, bmat, n, which, nev, tol, resid, 
+     &        ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, 
+     &        info )
+c
+         if (ido .eq. -1 .or. ido .eq. 1) then
+c
+c           %-------------------------------------------%
+c           | Perform matrix vector multiplication      |
+c           |                y <--- OP*x                |
+c           | The user should supply his/her own        |
+c           | matrix vector multiplication routine here |
+c           | that takes workd(ipntr(1)) as the input   |
+c           | vector, and return the matrix vector      |
+c           | product to workd(ipntr(2)).               | 
+c           %-------------------------------------------%
+c
+            call dgemv('N',n,n,1.0d0,a,maxn,
+     &           workd(ipntr(1)),1,0.0d0,workd(ipntr(2)),1)
+c
+c           %-----------------------------------------%
+c           | L O O P   B A C K to call DNAUPD again. |
+c           %-----------------------------------------%
+c
+            go to 10
+c
+         else if (ido .eq. 2) then
+
+            call dcopy (n, workd(ipntr(1)), 1, workd(ipntr(2)), 1)
+
+            go to 10
+
+         end if 
+c 
+c     %----------------------------------------%
+c     | Either we have convergence or there is |
+c     | an error.                              |
+c     %----------------------------------------%
+c
+      if ( info .lt. 0 ) then
+c
+c        %--------------------------%
+c        | Error message, check the |
+c        | documentation in DNAUPD. |
+c        %--------------------------%
+c
+         print *, ' '
+         print *, ' Error with _naupd, info = ', info
+         print *, ' Check the documentation of _naupd'
+         print *, ' '
+         stop 1
+c
+      else 
+c
+c        %-------------------------------------------%
+c        | No fatal errors occurred.                 |
+c        | Post-Process using DNEUPD.                |
+c        |                                           |
+c        | Computed eigenvalues may be extracted.    |
+c        |                                           |
+c        | Eigenvectors may also be computed now if  |
+c        | desired.  (indicated by rvec = .true.)    |
+c        %-------------------------------------------%
+c
+         rvec = .true.
+c
+         call dneupd ( rvec, 'A', select, d, d(1,2), v, ldv, 
+     &        sigmar, sigmai, workev, bmat, n, which, nev, tol, 
+     &        resid, ncv, v, ldv, iparam, ipntr, workd, workl,
+     &        lworkl, ierr )
+c
+c        %-----------------------------------------------%
+c        | The real part of the eigenvalue is returned   |
+c        | in the first column of the two dimensional    |
+c        | array D, and the imaginary part is returned   |
+c        | in the second column of D.  The corresponding |
+c        | eigenvectors are returned in the first NEV    |
+c        | columns of the two dimensional array V if     |
+c        | requested.  Otherwise, an orthogonal basis    |
+c        | for the invariant subspace corresponding to   |
+c        | the eigenvalues in D is returned in V.        |
+c        %-----------------------------------------------%
+c
+         if ( ierr .ne. 0) then
+c
+c           %------------------------------------%
+c           | Error condition:                   |
+c           | Check the documentation of DNEUPD. |
+c           %------------------------------------%
+c
+             print *, ' '
+             print *, ' Error with _neupd, info = ', ierr
+             print *, ' Check the documentation of _neupd. '
+             print *, ' '
+c
+         else 
+c
+             first  = .true.
+             nconv  = iparam(5)
+             do 20 j=1, nconv
+c
+c               %---------------------------%
+c               | Compute the residual norm |
+c               |                           |
+c               |   ||  A*x - lambda*x ||   |
+c               |                           |
+c               | for the NCONV accurately  |
+c               | computed eigenvalues and  |
+c               | eigenvectors.  (iparam(5) |
+c               | indicates how many are    |
+c               | accurate to the requested |
+c               | tolerance)                |
+c               %---------------------------%
+c
+                if (d(j,2) .eq. zero)  then
+c
+c                  %--------------------%
+c                  | Ritz value is real |
+c                  %--------------------%
+c
+                   call dgemv('N',n,n,1.0d0,a,maxn,
+     &                  v(1,j),1,0.0d0,ax,1)
+                   call daxpy(n, -d(j,1), v(1,j), 1, ax, 1)
+                   d(j,3) = dnrm2(n, ax, 1)
+                   d(j,3) = d(j,3) / abs(d(j,1))
+c
+                else if (first) then
+c
+c                  %------------------------%
+c                  | Ritz value is complex. |
+c                  | Residual of one Ritz   |
+c                  | value of the conjugate |
+c                  | pair is computed.      | 
+c                  %------------------------%
+c        
+                   call dgemv('N',n,n,1.0d0,a,maxn,
+     &                  v(1,j),1,0.0d0,ax,1)
+                   call daxpy(n, -d(j,1), v(1,j), 1, ax, 1)
+                   call daxpy(n, d(j,2), v(1,j+1), 1, ax, 1)
+                   d(j,3) = dnrm2(n, ax, 1)
+                   call dgemv('N',n,n,1.0d0,a,maxn,
+     &                  v(1,j+1),1,0.0d0,ax,1)
+                   call daxpy(n, -d(j,2), v(1,j), 1, ax, 1)
+                   call daxpy(n, -d(j,1), v(1,j+1), 1, ax, 1)
+                   d(j,3) = dlapy2( d(j,3), dnrm2(n, ax, 1) )
+                   d(j,3) = d(j,3) / dlapy2(d(j,1),d(j,2))
+                   d(j+1,3) = d(j,3)
+                   first = .false.
+                else
+                   first = .true.
+                end if
+c
+ 20          continue
+c
+c            %-----------------------------%
+c            | Display computed residuals. |
+c            %-----------------------------%
+c
+             call dmout(6, nconv, 3, d, maxncv, -6,
+     &            'Ritz values (Real,Imag) and relative residuals')
+          end if
+c
+c        %-------------------------------------------%
+c        | Print additional convergence information. |
+c        %-------------------------------------------%
+c
+         if ( info .eq. 1) then
+             print *, ' '
+             print *, ' Maximum number of iterations reached.'
+             print *, ' '
+         else if ( info .eq. 3) then
+             print *, ' ' 
+             print *, ' No shifts could be applied during implicit',
+     &                ' Arnoldi update, try increasing NCV.'
+             print *, ' '
+         end if      
+c
+         print *, ' '
+         print *, ' _NDRV1 '
+         print *, ' ====== '
+         print *, ' ' 
+         print *, ' Size of the matrix is ', n
+         print *, ' The number of Ritz values requested is ', nev
+         print *, ' The number of Arnoldi vectors generated',
+     &            ' (NCV) is ', ncv
+         print *, ' What portion of the spectrum: ', which
+         print *, ' The number of converged Ritz values is ', 
+     &              nconv 
+         print *, ' The number of Implicit Arnoldi update',
+     &            ' iterations taken is ', iparam(3)
+         print *, ' The number of OP*x is ', iparam(9)
+         print *, ' The convergence criterion is ', tol
+         print *, ' '
+c
+      end if
+c
+c     %---------------------------%
+c     | Done with program dndrv1. |
+c     %---------------------------%
+c
+ 9000 continue
+c
+      end
+c 

TESTS/icb_arpack_c.c

@@ -19,9 +19,6 @@
  * A*x = lambda*x where A is the diagonal matrix
  * with entries 1000, 999, ... , 2, 1 on the diagonal.
  * */
-#ifndef BLASINT
-#define BLASINT int
-#endif
 
 void dMatVec(double * x, double * y) {
   int i;
@@ -30,25 +27,25 @@ void dMatVec(double * x, double * y) {
 };
 
 int ds() {
-  BLASINT ido = 0;
+  a_int ido = 0;
   char bmat[] = "I";
-  BLASINT N = 1000;
+  a_int N = 1000;
   char which[] = "LM";
-  BLASINT nev = 9;
+  a_int nev = 9;
   double tol = 0;
   double resid[N];
-  BLASINT ncv = 2*nev+1;
+  a_int ncv = 2*nev+1;
   double V[ncv*N];
-  BLASINT ldv = N;
-  BLASINT iparam[11];
-  BLASINT ipntr[14];
+  a_int ldv = N;
+  a_int iparam[11];
+  a_int ipntr[14];
   double workd[3*N];
   bool rvec = true;
   char howmny[] = "A";
   double* d = (double*) malloc((nev+1)*sizeof(double));
-  int select[ncv];
+  a_int select[ncv];
   double z[(N+1)*(nev+1)];
-  BLASINT ldz = N+1;
+  a_int ldz = N+1;
   double sigma=0;
   int k;
   for (k=0; k < 3*N; ++k )
@@ -56,8 +53,8 @@ int ds() {
   double workl[3*(ncv*ncv) + 6*ncv];
   for (k=0; k < 3*(ncv*ncv) + 6*ncv; ++k )
     workl[k] = 0;
-  BLASINT lworkl = 3*(ncv*ncv) + 6*ncv;
-  BLASINT info = 0;
+  a_int lworkl = 3*(ncv*ncv) + 6*ncv;
+  a_int info = 0;
 
   iparam[0] = 1;
   iparam[2] = 10*N;
@@ -97,36 +94,36 @@ void zMatVec(double _Complex * x, double _Complex * y) {
 };
 
 int zn() {
-  BLASINT ido = 0;
+  a_int ido = 0;
   char bmat[] = "I";
-  BLASINT N = 1000;
+  a_int N = 1000;
   char which[] = "LM";
-  BLASINT nev = 9;
+  a_int nev = 9;
   double tol = 0;
   double _Complex resid[N];
-  BLASINT ncv = 2*nev+1;
+  a_int ncv = 2*nev+1;
   double _Complex V[ncv*N];
-  BLASINT ldv = N;
-  BLASINT iparam[11];
-  BLASINT ipntr[14];
+  a_int ldv = N;
+  a_int iparam[11];
+  a_int ipntr[14];
   double _Complex workd[3*N];
   bool rvec = true;
   char howmny[] = "A";
   double _Complex* d = (double _Complex*) malloc((nev+1)*sizeof(double _Complex));
-  int select[ncv];
+  a_int select[ncv];
   double _Complex z[(N+1)*(nev+1)];
-  BLASINT ldz = N+1;
-  double sigma=0;
+  a_int ldz = N+1;
+  double _Complex sigma=0. + I*0.;
   int k;
   for (k=0; k < 3*N; ++k )
     workd[k] = 0;
   double _Complex workl[3*(ncv*ncv) + 6*ncv];
   for (k=0; k < 3*(ncv*ncv) + 6*ncv; ++k )
     workl[k] = 0;
-  BLASINT lworkl = 3*(ncv*ncv) + 6*ncv;
+  a_int lworkl = 3*(ncv*ncv) + 6*ncv;
   double rwork[ncv];
   double _Complex workev[2*ncv];
-  BLASINT info = 0;
+  a_int info = 0;
 
   iparam[0] = 1;
   iparam[2] = 10*N;

TESTS/icb_arpack_cpp.cpp

@@ -19,10 +19,6 @@
 #include "debug_c.hpp"  // debug arpack.
 #include "stat_c.hpp"   // arpack statistics.
 
-#ifndef BLASINT
-#define BLASINT int
-#endif
-
 template<typename Real>
 void diagonal_matrix_vector_product(Real const* const x, Real* const y) {
   for (int i = 0; i < 1000; ++i) {
@@ -32,15 +28,15 @@ void diagonal_matrix_vector_product(Real const* const x, Real* const y) {
 
 template<typename Real>
 void real_symmetric_runner() {
-  BLASINT const N = 1000;
-  BLASINT const nev = 9;
+  a_int const N = 1000;
+  a_int const nev = 9;
 
-  BLASINT const ncv = 2 * nev + 1;
-  BLASINT const ldv = N;
+  a_int const ncv = 2 * nev + 1;
+  a_int const ldv = N;
 
-  BLASINT const ldz = N + 1;
+  a_int const ldz = N + 1;
 
-  BLASINT const lworkl = 3 * (ncv * ncv) + 6 * ncv;
+  a_int const lworkl = 3 * (ncv * ncv) + 6 * ncv;
 
   Real const tol = 0.0;
   Real const sigma = 0.0;
@@ -54,7 +50,7 @@ void real_symmetric_runner() {
   std::vector<Real> d((nev + 1));
   std::vector<Real> z((N + 1) * (nev + 1));
 
-  std::array<BLASINT, 11> iparam{};
+  std::array<a_int, 11> iparam{};
 
   iparam[0] = 1;
   iparam[2] = 10 * N;
@@ -62,9 +58,9 @@ void real_symmetric_runner() {
   iparam[4] = 0;  // number of ev found by arpack.
   iparam[6] = 1;
 
-  std::array<BLASINT, 14> ipntr{};
+  std::array<a_int, 14> ipntr{};
 
-  BLASINT info = 0, ido = 0;
+  a_int info = 0, ido = 0;
 
   while (ido != 99) {
     arpack::saupd(ido, arpack::bmat::identity, N,
@@ -81,7 +77,7 @@ void real_symmetric_runner() {
     throw std::domain_error("Error inside ARPACK routines");
   }
 
-  std::vector<int> select(ncv);
+  std::vector<a_int> select(ncv);
 
   arpack::seupd(rvec, arpack::howmny::ritz_vectors, select.data(), d.data(),
                 z.data(), ldz, sigma, arpack::bmat::identity, N,
@@ -109,18 +105,18 @@ void diagonal_matrix_vector_product(std::complex<Real> const* const x,
 
 template<typename Real>
 void complex_symmetric_runner() {
-  BLASINT const N = 1000;
-  BLASINT const nev = 9;
+  a_int const N = 1000;
+  a_int const nev = 9;
 
-  BLASINT const ncv = 2 * nev + 1;
-  BLASINT const ldv = N;
+  a_int const ncv = 2 * nev + 1;
+  a_int const ldv = N;
 
-  BLASINT const ldz = N + 1;
+  a_int const ldz = N + 1;
 
-  BLASINT const lworkl = 3 * (ncv * ncv) + 6 * ncv;
+  a_int const lworkl = 3 * (ncv * ncv) + 6 * ncv;
 
   Real const tol = 0.0;
-  Real const sigma = 0.0;
+  std::complex<Real> const sigma(0.0, 0.0);
 
   bool const rvec = true;
 
@@ -133,16 +129,16 @@ void complex_symmetric_runner() {
   std::vector<Real> rwork(ncv);
   std::vector<std::complex<Real>> workev(2 * ncv);
 
-  std::array<BLASINT, 11> iparam{};
+  std::array<a_int, 11> iparam{};
   iparam[0] = 1;
   iparam[2] = 10 * N;
   iparam[3] = 1;
   iparam[4] = 0;  // number of ev found by arpack.
   iparam[6] = 1;
 
-  std::array<BLASINT, 14> ipntr{};
+  std::array<a_int, 14> ipntr{};
 
-  BLASINT info = 0, ido = 0;
+  a_int info = 0, ido = 0;
 
   while (ido != 99) {
     arpack::naupd(ido, arpack::bmat::identity, N,
@@ -159,7 +155,7 @@ void complex_symmetric_runner() {
     throw std::domain_error("Error inside ARPACK routines");
   }
 
-  std::vector<int> select(ncv);
+  std::vector<a_int> select(ncv);
 
   arpack::neupd(rvec, arpack::howmny::ritz_vectors, select.data(), d.data(),
                 z.data(), ldz, sigma, workev.data(), arpack::bmat::identity, N,

arpack.h

@@ -1,25 +0,0 @@
-#ifndef __ARPACK_H__
-#define __ARPACK_H__
-
-#ifdef __cplusplus
-extern "C" {
-#endif
-
-void cnaupd_c(int* ido, char const* bmat, int n, char const* which, int nev, float tol, float _Complex* resid, int ncv, float _Complex* v, int ldv, int* iparam, int* ipntr, float _Complex* workd, float _Complex* workl, int lworkl, float* rwork, int* info);
-void cneupd_c(bool rvec, char const* howmny, int const* select, float _Complex* d, float _Complex* z, int ldz, float _Complex sigma, float _Complex* workev, char const* bmat, int n, char const* which, int nev, float tol, float _Complex* resid, int ncv, float _Complex* v, int ldv, int* iparam, int* ipntr, float _Complex* workd, float _Complex* workl, int lworkl, float* rwork, int* info);
-void dnaupd_c(int* ido, char const* bmat, int n, char const* which, int nev, double tol, double* resid, int ncv, double* v, int ldv, int* iparam, int* ipntr, double* workd, double* workl, int lworkl, int* info);
-void dneupd_c(bool rvec, char const* howmny, int const* select, double* dr, double* di, double* z, int ldz, double sigmar, double sigmai, char const* bmat, int n, char const* which, int nev, double tol, double* resid, int ncv, double* v, int ldv, int* iparam, int* ipntr, double* workd, double* workl, int lworkl, int* info);
-void dsaupd_c(int* ido, char const* bmat, int n, char const* which, int nev, double tol, double* resid, int ncv, double* v, int ldv, int* iparam, int* ipntr, double* workd, double* workl, int lworkl, int* info);
-void dseupd_c(bool rvec, char const* howmny, int const* select, double* d, double* z, int ldz, double sigma, char const* bmat, int n, char const* which, int nev, double tol, double* resid, int ncv, double* v, int ldv, int* iparam, int* ipntr, double* workd, double* workl, int lworkl, int* info);
-void snaupd_c(int* ido, char const* bmat, int n, char const* which, int nev, float tol, float* resid, int ncv, float* v, int ldv, int* iparam, int* ipntr, float* workd, float* workl, int lworkl, int* info);
-void sneupd_c(bool rvec, char const* howmny, int const* select, float* dr, float* di, float* z, int ldz, float sigmar, float sigmai, char const* bmat, int n, char const* which, int nev, float tol, float* resid, int ncv, float* v, int ldv, int* iparam, int* ipntr, float* workd, float* workl, int lworkl, int* info);
-void ssaupd_c(int* ido, char const* bmat, int n, char const* which, int nev, float tol, float* resid, int ncv, float* v, int ldv, int* iparam, int* ipntr, float* workd, float* workl, int lworkl, int* info);
-void sseupd_c(bool rvec, char const* howmny, int const* select, float* d, float* z, int ldz, float sigma, char const* bmat, int n, char const* which, int nev, float tol, float* resid, int ncv, float* v, int ldv, int* iparam, int* ipntr, float* workd, float* workl, int lworkl, int* info);
-void znaupd_c(int* ido, char const* bmat, int n, char const* which, int nev, double tol, double _Complex* resid, int ncv, double _Complex* v, int ldv, int* iparam, int* ipntr, double _Complex* workd, double _Complex* workl, int lworkl, double* rwork, int* info);
-void zneupd_c(bool rvec, char const* howmny, int const* select, double _Complex* d, double _Complex* z, int ldz, double _Complex sigma, double _Complex* workev, char const* bmat, int n, char const* which, int nev, double tol, double _Complex* resid, int ncv, double _Complex* v, int ldv, int* iparam, int* ipntr, double _Complex* workd, double _Complex* workl, int lworkl, double* rwork, int* info);
-
-#ifdef  __cplusplus
-}
-#endif
-
-#endif

arpack.pc.in

@@ -1,7 +1,7 @@
 prefix=@prefix@
 exec_prefix=@exec_prefix@
 libdir=@libdir@
-includedir=@includedir@/arpack
+includedir=@includedir@
 
 Name: @PACKAGE_NAME@
 Description: Collection of Fortran77 subroutines designed to solve large scale eigenvalue problems
@@ -9,4 +9,4 @@ Version: @PACKAGE_VERSION@
 URL: @PACKAGE_URL@
 Libs: -L${libdir} -larpack@LIBSUFFIX@
 Libs.private: @LAPACK_LIBS@ @BLAS_LIBS@
-Cflags: -I${includedir}
+Cflags: -I${includedir}/arpack

arpackdef.h.in

@@ -0,0 +1,15 @@
+#ifndef __ARPACKDEF_H__
+#define __ARPACKDEF_H__
+
+#define INTERFACE64 @INTERFACE64@
+
+#if INTERFACE64
+#define c_int  c_int64_t
+#define a_int    int64_t
+#define a_uint  uint64_t
+#else
+#define a_int            int
+#define a_uint  unsigned int
+#endif
+
+#endif

configure.ac

@@ -1,9 +1,8 @@
 AC_PREREQ(2.67)
-AC_INIT([ARPACK-NG],[3.6.0],[https://github.com/opencollab/arpack-ng/issues/],[arpack-ng],[https://github.com/opencollab/arpack-ng/])
+AC_INIT([ARPACK-NG],[3.7.0],[https://github.com/opencollab/arpack-ng/issues/],[arpack-ng],[https://github.com/opencollab/arpack-ng/])
 AC_CONFIG_SRCDIR([SRC/version.h])
 AC_CONFIG_AUX_DIR([build-aux])
 AC_CONFIG_MACRO_DIR([m4])
-AC_CONFIG_LINKS([TESTS/testA.mtx:TESTS/testA.mtx])
 AM_INIT_AUTOMAKE([foreign])
 AM_MAINTAINER_MODE
 
@@ -12,6 +11,7 @@ dnl Autoconf was designed with the fundamental fact that compiler (= language su
 AC_PROG_FC
 AC_PROG_CC
 AC_PROG_CXX
+PKG_PROG_PKG_CONFIG
 
 dnl Check for BLAS libraries
 AX_BLAS([], [AC_MSG_ERROR([cannot find BLAS libraries])])
@@ -21,10 +21,19 @@ AX_LAPACK([], [AC_MSG_ERROR([cannot find LAPACK libraries])])
 
 AC_ARG_VAR(INTERFACE64, [set to 1 to use the 64-bit integer interface (ILP64) for ARPACK, BLAS and LAPACK])
 if test x"$INTERFACE64" == x"1"; then
+  AC_LANG_PUSH([Fortran 77])
   AX_CHECK_COMPILE_FLAG(-fdefault-integer-8, FFLAGS="$FFLAGS -fdefault-integer-8",
                         AX_CHECK_COMPILE_FLAG(-i8, FFLAGS="$FFLAGS -i8",
                                               AC_MSG_WARN([configure does not know how to make your Fortran compiler use 64-bit integers: set it manually via FFLAGS.])))
-  CFLAGS="$CFLAGS -DBLASINT=int64_t"
+  AC_LANG_POP([Fortran 77])
+
+  AC_LANG_PUSH([Fortran])
+  AX_CHECK_COMPILE_FLAG(-fdefault-integer-8, FCFLAGS="$FCFLAGS -fdefault-integer-8",
+                        AX_CHECK_COMPILE_FLAG(-i8, FCFLAGS="$FCFLAGS -i8",
+                                              AC_MSG_WARN([configure does not know how to make your Fortran compiler use 64-bit integers: set it manually via FCFLAGS.])))
+  AC_LANG_POP([Fortran])
+else
+  INTERFACE64=0
 fi
 
 AC_ARG_VAR(LIBSUFFIX, [suffix to add to ARPACK libraries names])
@@ -53,17 +62,40 @@ AS_IF([test x"$enable_mpi" != x"no"], [
 	AC_LANG_POP([Fortran 77])
 ])
 
-AM_CONDITIONAL([MPI], [test x"$enable_mpi" != x"no"])
-AM_CONDITIONAL([BLACS], [false])
-
 dnl See if compiling with ISO_C_BINDING support
 AC_ARG_ENABLE([icb],
-              [AC_HELP_STRING([--enable-icb],
+              [AS_HELP_STRING([--enable-icb],
                               [build support for *[ae]upd_c with ISO_C_BINDING])],
-              [AS_VAR_SET([enable_icb], [$enableval])],
+              [],
               [AS_VAR_SET([enable_icb], [no])])
 
-AM_CONDITIONAL([ICB], [test x"$enable_icb" != x"no"])
+dnl See if compiling matrix market example based on ICB
+AC_ARG_ENABLE([icb-exmm],
+              [AS_HELP_STRING([--enable-icb-exmm],
+                              [matrix market example based on ISO_C_BINDING])],
+              [AS_VAR_SET([enable_icb], [$enableval])],
+              [AS_VAR_SET([enable_icb_exmm], [no])])
+
+if test x"$enable_icb_exmm" != x"no"; then
+	PKG_CHECK_MODULES([EIGEN3], [eigen3 >= 3.2])
+	AC_LANG_PUSH([C++])
+	CPPFLAGS_SAVE=$CPPFLAGS
+	CPPFLAGS=$EIGEN3_CFLAGS
+	AC_CHECK_HEADER([Eigen/Sparse],
+		[],
+		[AC_MSG_ERROR([Eigen/Sparse not found])])
+	AC_CHECK_HEADER([Eigen/IterativeLinearSolvers],
+		[],
+		[AC_MSG_ERROR([Eigen/IterativeLinearSolvers not found])])
+	AC_CHECK_HEADER([Eigen/SparseLU],
+		[],
+		[AC_MSG_ERROR([Eigen/SparseLU not found])])
+	AC_CHECK_HEADER([Eigen/SparseQR],
+		[],
+		[AC_MSG_ERROR([Eigen/SparseQR not found])])
+	CPPFLAGS=$CPPFLAGS_SAVE
+	AC_LANG_POP([C++])
+fi
 
 if test x"$enable_icb" != x"no"; then
   AC_LANG_PUSH([Fortran])
@@ -88,8 +120,6 @@ if test x"$enable_icb" != x"no"; then
                 )
   AC_LANG_POP([Fortran])
 
-  AC_CONFIG_FILES([arpack.h:arpack.h])
-  AC_CONFIG_FILES([arpack.hpp:arpack.hpp])
   if test x"$enable_mpi" != x"no"; then
     AC_LANG_PUSH([Fortran])
     AX_MPI([], AC_MSG_ERROR([could not compile a Fortran MPI test program]))
@@ -97,13 +127,12 @@ if test x"$enable_icb" != x"no"; then
     AC_FC_WRAPPERS dnl set FCLIBS
     AC_LANG_POP([Fortran])
 
-    AC_CONFIG_FILES([parpack.h:parpack.h])
-    AC_CONFIG_FILES([parpack.hpp:parpack.hpp])
     AC_LANG_PUSH([C])
     AX_MPI([], AC_MSG_ERROR([could not compile a C MPI test program]))
     AC_SUBST([MPI_C_LIBS], ["$MPILIBS $FCLIBS"])
     CC=$MPICC
     AC_LANG_POP([C])
+
     AC_LANG_PUSH([C++])
     AX_MPI([], AC_MSG_ERROR([could not compile a C++ MPI test program]))
     AC_SUBST([MPI_CXX_LIBS], ["$MPILIBS $FCLIBS"])
@@ -136,9 +165,11 @@ else
     AC_FC_LINE_LENGTH(unlimited)
     FFLAGS="$FFLAGS $ac_cv_fc_line_length"
 
+    AC_LANG_PUSH([Fortran 77])
     AX_CHECK_COMPILE_FLAG(-cpp, FFLAGS="$FFLAGS -cpp",
                           AX_CHECK_COMPILE_FLAG(-fpp, FFLAGS="$FFLAGS -fpp",
                                                 AC_MSG_WARN([configure does not know how to enable your Fortran compiler's preprocessor: set it manually via FFLAGS.])))
+    AC_LANG_POP([Fortran 77])
 
     SCALARFUNS="naitr napps naup2 naupd nconv neigh neupd ngets statn saitr sapps saup2 saupd sconv seigt seupd sgets stats getv0 sortc sortr sesrt stqrb"
     COMPLEXFUNS="naitr napps naup2 naupd neigh neupd ngets statn getv0 sortc"
@@ -182,30 +213,22 @@ else
   AC_F77_FUNC(dneupd$SYMBOLSUFFIX, dneupdsuff)
   AC_F77_FUNC(cheev$SYMBOLSUFFIX, cheevsuff)
   dnl Needed for BLAS check and for tests (even when suffix is empty)
-  CFLAGS="$CFLAGS -Dsgemm=$sgemmsuff -Dsnaupd=$snaupdsuff -Dsneupd=$sneupdsuff -Ddnaupd=$dnaupdsuff -Ddneupd=$dneupdsuff -Dcheev=$cheevsuff"
+  CPPFLAGS="$CPPFLAGS -Dsgemm=$sgemmsuff -Dsnaupd=$snaupdsuff -Dsneupd=$sneupdsuff -Ddnaupd=$dnaupdsuff -Ddneupd=$dneupdsuff -Dcheev=$cheevsuff"
 fi
 
-dnl Add support for debuging arpack
-AC_CONFIG_FILES([debug.h:debug.h])
-if test x"$enable_icb" != x"no"; then
-  AC_CONFIG_FILES([debug_c.h:debug_c.h])
-  AC_CONFIG_FILES([debug_c.hpp:debug_c.hpp])
-fi
-
-dnl Add support for arpack statistics
-AC_CONFIG_FILES([stat.h:stat.h])
-if test x"$enable_icb" != x"no"; then
-  AC_CONFIG_FILES([stat_c.h:stat_c.h])
-  AC_CONFIG_FILES([stat_c.hpp:stat_c.hpp])
-fi
+AM_CONDITIONAL([MPI], [test x"$enable_mpi" != x"no"])
+AM_CONDITIONAL([BLACS], [false])
+AM_CONDITIONAL([ICB], [test x"$enable_icb" != x"no"])
+AM_CONDITIONAL([ICBEXMM], [test x"$enable_icb_exmm" != x"no"])
 
 m4_ifdef([PKG_INSTALLDIR], [PKG_INSTALLDIR], [AC_SUBST([pkgconfigdir], [${libdir}/pkgconfig])])
 AC_CONFIG_FILES([arpack$LIBSUFFIX.pc:arpack.pc.in], [], [LIBSUFFIX="$LIBSUFFIX"])
-
 AC_CONFIG_FILES([
+	arpackdef.h
 	Makefile
-	SRC/Makefile
+	ICB/Makefile
 	UTIL/Makefile
+	SRC/Makefile
 	TESTS/Makefile
 	EXAMPLES/Makefile
 	EXAMPLES/BAND/Makefile
@@ -214,6 +237,7 @@ AC_CONFIG_FILES([
 	EXAMPLES/SIMPLE/Makefile
 	EXAMPLES/SVD/Makefile
 	EXAMPLES/SYM/Makefile
+	EXAMPLES/MATRIX_MARKET/Makefile
 	PARPACK/Makefile
 	PARPACK/SRC/MPI/Makefile
 	PARPACK/SRC/BLACS/Makefile
@@ -223,23 +247,38 @@ AC_CONFIG_FILES([
 	PARPACK/EXAMPLES/MPI/Makefile
 	PARPACK/EXAMPLES/BLACS/Makefile
 ])
+AC_CONFIG_LINKS([
+	TESTS/testA.mtx:TESTS/testA.mtx
+	EXAMPLES/MATRIX_MARKET/An.mtx:EXAMPLES/MATRIX_MARKET/An.mtx
+	EXAMPLES/MATRIX_MARKET/As.mtx:EXAMPLES/MATRIX_MARKET/As.mtx
+	EXAMPLES/MATRIX_MARKET/Az.mtx:EXAMPLES/MATRIX_MARKET/Az.mtx
+	EXAMPLES/MATRIX_MARKET/B.mtx:EXAMPLES/MATRIX_MARKET/B.mtx
+	EXAMPLES/MATRIX_MARKET/Bz.mtx:EXAMPLES/MATRIX_MARKET/Bz.mtx
+])
 AC_OUTPUT
 
 AC_MSG_RESULT([
 --------------------------------------------------
-Configuration summary for $PACKAGE_NAME ($VERSION)
+Configuration summary for $PACKAGE_STRING
 --------------------------------------------------
 Installation prefix : $prefix
 MPI enabled         : $enable_mpi
 ICB enabled         : $enable_icb
+INTERFACE64         : $INTERFACE64
+F77                 : $F77
+FFLAGS              : $FFLAGS
 FC                  : $FC
 FCFLAGS             : $FCFLAGS
 CC                  : $CC
 CFLAGS              : $CFLAGS
 CXX                 : $CXX
 CXXFLAGS            : $CXXFLAGS
+CPPFLAGS            : $CPPFLAGS
 BLAS                : $BLAS_LIBS
 LAPACK              : $LAPACK_LIBS
+EIGEN               : $EIGEN3_CFLAGS
+LIBS                : $LIBS
+LDADD               : $LDADD
 --------------------------------------------------
 Configuration OK
 --------------------------------------------------

dbgini.f

@@ -0,0 +1,13 @@
+c Initialisation of the debug common block to "no debug".
+      block data dbgini
+      common /debug/ logfil, ndigit, mgetv0,
+     &     msaupd, msaup2, msaitr, mseigt, msapps, msgets, mseupd,
+     &     mnaupd, mnaup2, mnaitr, mneigh, mnapps, mngets, mneupd,
+     &     mcaupd, mcaup2, mcaitr, mceigh, mcapps, mcgets, mceupd
+      data logfil, ndigit, mgetv0,
+     &     msaupd, msaup2, msaitr, mseigt, msapps, msgets, mseupd,
+     &     mnaupd, mnaup2, mnaitr, mneigh, mnapps, mngets, mneupd,
+     &     mcaupd, mcaup2, mcaitr, mceigh, mcapps, mcgets, mceupd
+     &     / 6, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+     &     0, 0, 0, 0, 0, 0 /
+      end

debug_init.f90

@@ -1,12 +0,0 @@
-! Initialisation of the debug common block to "no debug".
-block data debug_init
-  common /debug/ logfil, ndigit, mgetv0,                                &
-                 msaupd, msaup2, msaitr, mseigt, msapps, msgets, mseupd,&
-                 mnaupd, mnaup2, mnaitr, mneigh, mnapps, mngets, mneupd,&
-                 mcaupd, mcaup2, mcaitr, mceigh, mcapps, mcgets, mceupd
-  data logfil, ndigit, mgetv0,                                &
-       msaupd, msaup2, msaitr, mseigt, msapps, msgets, mseupd,&
-       mnaupd, mnaup2, mnaitr, mneigh, mnapps, mngets, mneupd,&
-       mcaupd, mcaup2, mcaitr, mceigh, mcapps, mcgets, mceupd &
-       / 6, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 /
-end block data debug_init

parpack.h

@@ -1,32 +0,0 @@
-#ifndef __PARPACK_H__
-#define __PARPACK_H__
-
-/*
- * IMPORTANT: MPI communicators MUST be passed from C to Fortran using MPI_Comm_c2f.
- *            MPI_Fint MCW = MPI_Comm_c2f(MPI_COMM_WORLD);
- */
-
-#include "mpi.h"
-
-#ifdef __cplusplus
-extern "C" {
-#endif
-
-void pcnaupd_c(MPI_Fint comm, int* ido, char const* bmat, int n, char const* which, int nev, float tol, float _Complex* resid, int ncv, float _Complex* v, int ldv, int* iparam, int* ipntr, float _Complex* workd, float _Complex* workl, int lworkl, float _Complex* rwork, int* info);
-void pcneupd_c(MPI_Fint comm, bool rvec, char const* howmny, int const* select, float _Complex* d, float _Complex* z, int ldz, float _Complex sigma, float _Complex* workev, char const* bmat, int n, char const* which, int nev, float tol, float _Complex* resid, int ncv, float _Complex* v, int ldv, int* iparam, int* ipntr, float _Complex* workd, float _Complex* workl, int lworkl, float _Complex* rwork, int* info);
-void pdnaupd_c(MPI_Fint comm, int* ido, char const* bmat, int n, char const* which, int nev, double tol, double* resid, int ncv, double* v, int ldv, int* iparam, int* ipntr, double* workd, double* workl, int lworkl, int* info);
-void pdneupd_c(MPI_Fint comm, bool rvec, char const* howmny, int const* select, double* dr, double* di, double* z, int ldz, double sigmar, double sigmai, char const* bmat, int n, char const* which, int nev, double tol, double* resid, int ncv, double* v, int ldv, int* iparam, int* ipntr, double* workd, double* workl, int lworkl, int* info);
-void pdsaupd_c(MPI_Fint comm, int* ido, char const* bmat, int n, char const* which, int nev, double tol, double* resid, int ncv, double* v, int ldv, int* iparam, int* ipntr, double* workd, double* workl, int lworkl, int* info);
-void pdseupd_c(MPI_Fint comm, bool rvec, char const* howmny, int const* select, double* d, double* z, int ldz, double sigma, char const* bmat, int n, char const* which, int nev, double tol, double* resid, int ncv, double* v, int ldv, int* iparam, int* ipntr, double* workd, double* workl, int lworkl, int* info);
-void psnaupd_c(MPI_Fint comm, int* ido, char const* bmat, int n, char const* which, int nev, float tol, float* resid, int ncv, float* v, int ldv, int* iparam, int* ipntr, float* workd, float* workl, int lworkl, int* info);
-void psneupd_c(MPI_Fint comm, bool rvec, char const* howmny, int const* select, float* dr, float* di, float* z, int ldz, float sigmar, float sigmai, char const* bmat, int n, char const* which, int nev, float tol, float* resid, int ncv, float* v, int ldv, int* iparam, int* ipntr, float* workd, float* workl, int lworkl, int* info);
-void pssaupd_c(MPI_Fint comm, int* ido, char const* bmat, int n, char const* which, int nev, float tol, float* resid, int ncv, float* v, int ldv, int* iparam, int* ipntr, float* workd, float* workl, int lworkl, int* info);
-void psseupd_c(MPI_Fint comm, bool rvec, char const* howmny, int const* select, float* d, float* z, int ldz, float sigma, char const* bmat, int n, char const* which, int nev, float tol, float* resid, int ncv, float* v, int ldv, int* iparam, int* ipntr, float* workd, float* workl, int lworkl, int* info);
-void pznaupd_c(MPI_Fint comm, int* ido, char const* bmat, int n, char const* which, int nev, double tol, double _Complex* resid, int ncv, double _Complex* v, int ldv, int* iparam, int* ipntr, double _Complex* workd, double _Complex* workl, int lworkl, double _Complex* rwork, int* info);
-void pzneupd_c(MPI_Fint comm, bool rvec, char const* howmny, int const* select, double _Complex* d, double _Complex* z, int ldz, double _Complex sigma, double _Complex* workev, char const* bmat, int n, char const* which, int nev, double tol, double _Complex* resid, int ncv, double _Complex* v, int ldv, int* iparam, int* ipntr, double _Complex* workd, double _Complex* workl, int lworkl, double _Complex* rwork, int* info);
-
-#ifdef  __cplusplus
-}
-#endif
-
-#endif

staini.f

@@ -0,0 +1,19 @@
+c     Initialisation of the stat common block.
+      block data staini
+      common /timing/
+     &     nopx, nbx, nrorth, nitref, nrstrt,
+     &     tsaupd, tsaup2, tsaitr, tseigt, tsgets, tsapps, tsconv,
+     &     tnaupd, tnaup2, tnaitr, tneigh, tngets, tnapps, tnconv,
+     &     tcaupd, tcaup2, tcaitr, tceigh, tcgets, tcapps, tcconv,
+     &     tmvopx, tmvbx, tgetv0, titref, trvec
+      data nopx, nbx, nrorth, nitref, nrstrt,
+     &     tsaupd, tsaup2, tsaitr, tseigt, tsgets, tsapps, tsconv,
+     &     tnaupd, tnaup2, tnaitr, tneigh, tngets, tnapps, tnconv,
+     &     tcaupd, tcaup2, tcaitr, tceigh, tcgets, tcapps, tcconv,
+     &     tmvopx, tmvbx, tgetv0, titref, trvec
+     &     / 0, 0, 0, 0, 0,
+     &     0., 0., 0., 0., 0., 0., 0.,
+     &     0., 0., 0., 0., 0., 0., 0.,
+     &     0., 0., 0., 0., 0., 0., 0.,
+     &     0., 0., 0., 0., 0. /
+      end