Update contribution

Without this line, pkg-config 0.29.2 refuses to parse the arpack.pc file:

# pkg-config --libs arpack
Variable 'includedir' not defined in '/usr/lib/pkgconfig/arpack.pc'

.clang-format

@@ -0,0 +1 @@
+BasedOnStyle:  Google

.gitignore

@@ -3,16 +3,11 @@
 # Generated by `autoreconf`
 Makefile.in
 aclocal.m4
+configure
 autom4te.cache/
+build-aux/
 m4/libtool.m4
 m4/lt*.m4
-configure
-config.guess
-config.sub
-depcomp
-install-sh
-ltmain.sh
-missing
 
 # Generated by `./configure`
 Makefile
@@ -24,15 +19,27 @@ arpack.pc
 
 # Generated by `make`
 *.o
+*.mod
 *.lo
 *.la
 .libs/
 
+# CMake build directories
+BUILD
+Build
+build
+
 # Generated by `make check`
+*.log
+*.trs
 TESTS/dnsimp
 TESTS/bug_1315_double
 TESTS/bug_1315_single
 TESTS/bug_1323
+TESTS/bug_58_double
+TESTS/bug_79_double_complex
+TESTS/icb_arpack_c
+TESTS/icb_arpack_cpp
 EXAMPLES/BAND/[sd][sn]bdr[123456]
 EXAMPLES/BAND/[cz]nbdr[1234]
 EXAMPLES/COMPLEX/[cz]ndrv[1234]
@@ -44,6 +51,9 @@ EXAMPLES/SYM/[sd]sdrv[123456]
 PARPACK/EXAMPLES/MPI/p[sd]ndrv[13]
 PARPACK/EXAMPLES/MPI/p[sd]sdrv1
 PARPACK/EXAMPLES/MPI/p[cz]ndrv1
+PARPACK/TESTS/MPI/icb_parpack_c
+PARPACK/TESTS/MPI/icb_parpack_cpp
+PARPACK/TESTS/MPI/issue46
 
 *.lib
 *.def

.travis.yml

@@ -9,14 +9,67 @@ addons:
   apt:
     packages:
       - gfortran
-      - libopenblas-dev
+      - gcc
+      - g++
+      - libblas-dev
       - liblapack-dev
       - libopenmpi-dev
       - openmpi-bin
 
+env:
+  - BUILD=cmake
+  - BUILD=configure
+
 script:
-  - ./bootstrap
-  - ./configure --enable-mpi
-  - make all
-  - make check
-  - make distcheck
+  # CMAKE
+  - if [ "$BUILD" == "cmake" ]; then mkdir build ; fi
+  - if [ "$BUILD" == "cmake" ]; then cd build ; fi
+  - if [ "$BUILD" == "cmake" ]; then cmake -D EXAMPLES=ON -D MPI=ON  -D ICB=OFF .. ; fi
+  - if [ "$BUILD" == "cmake" ]; then make all ; fi
+  - if [ "$BUILD" == "cmake" ]; then make test ; fi
+  - if [ "$BUILD" == "cmake" ]; then rm -fr * ; fi # Restart cmake from scratch.
+  - if [ "$BUILD" == "cmake" ]; then cmake -D EXAMPLES=ON -D MPI=ON  -D ICB=ON  .. ; fi
+  - if [ "$BUILD" == "cmake" ]; then make all ; fi
+  - if [ "$BUILD" == "cmake" ]; then make test ; fi
+  - if [ "$BUILD" == "cmake" ]; then rm -fr * ; fi # Restart cmake from scratch.
+  - if [ "$BUILD" == "cmake" ]; then cmake -D EXAMPLES=ON -D MPI=OFF -D ICB=OFF .. ; fi
+  - if [ "$BUILD" == "cmake" ]; then make all ; fi
+  - if [ "$BUILD" == "cmake" ]; then make test ; fi
+  - if [ "$BUILD" == "cmake" ]; then rm -fr * ; fi # Restart cmake from scratch.
+  - if [ "$BUILD" == "cmake" ]; then cmake -D EXAMPLES=ON -D MPI=OFF -D ICB=ON  .. ; fi
+  - if [ "$BUILD" == "cmake" ]; then make all ; fi
+  - if [ "$BUILD" == "cmake" ]; then make test ; fi
+  - if [ "$BUILD" == "cmake" ]; then rm -fr * ; fi # Restart cmake from scratch.
+  - if [ "$BUILD" == "cmake" ]; then cmake -D BUILD_SHARED_LIBS=ON -D MPI=ON  -D ICB=OFF .. ; fi
+  - if [ "$BUILD" == "cmake" ]; then make all ; fi
+  - if [ "$BUILD" == "cmake" ]; then make test ; fi
+  - if [ "$BUILD" == "cmake" ]; then rm -fr * ; fi # Restart cmake from scratch.
+  - if [ "$BUILD" == "cmake" ]; then cmake -D BUILD_SHARED_LIBS=ON -D MPI=ON  -D ICB=ON  .. ; fi
+  - if [ "$BUILD" == "cmake" ]; then make all ; fi
+  - if [ "$BUILD" == "cmake" ]; then make test ; fi
+  - if [ "$BUILD" == "cmake" ]; then rm -fr * ; fi # Restart cmake from scratch.
+  - if [ "$BUILD" == "cmake" ]; then cmake -D BUILD_SHARED_LIBS=ON -D MPI=OFF -D ICB=OFF .. ; fi
+  - if [ "$BUILD" == "cmake" ]; then make all ; fi
+  - if [ "$BUILD" == "cmake" ]; then make test ; fi
+  - if [ "$BUILD" == "cmake" ]; then rm -fr * ; fi # Restart cmake from scratch.
+  - if [ "$BUILD" == "cmake" ]; then cmake -D BUILD_SHARED_LIBS=ON -D MPI=OFF -D ICB=ON  .. ; fi
+  - if [ "$BUILD" == "cmake" ]; then make all ; fi
+  - if [ "$BUILD" == "cmake" ]; then make test ; fi
+  - if [ "$BUILD" == "cmake" ]; then rm -fr * ; fi # Restart cmake from scratch for code coverage
+  - if [ "$BUILD" == "cmake" ]; then cmake -D EXAMPLES=ON -DCOVERALLS=ON -D MPI=ON -D ICB=ON ..  ; fi
+  - if [ "$BUILD" == "cmake" ]; then make all test coveralls VERBOSE=1 ; fi
+  # CONFIGURE
+  - if [ "$BUILD" == "configure" ]; then ./bootstrap ; fi
+  - if [ "$BUILD" == "configure" ]; then ./configure --enable-mpi ; fi
+  - if [ "$BUILD" == "configure" ]; then make all ; fi
+  - if [ "$BUILD" == "configure" ]; then make check ; fi
+  - if [ "$BUILD" == "configure" ]; then make distcheck ; fi
+  - if [ "$BUILD" == "configure" ]; then ./configure --enable-mpi --enable-icb ; fi
+  - if [ "$BUILD" == "configure" ]; then make all ; fi
+  - if [ "$BUILD" == "configure" ]; then make check ; fi
+  - if [ "$BUILD" == "configure" ]; then ./configure --enable-icb ; fi
+  - if [ "$BUILD" == "configure" ]; then make all ; fi
+  - if [ "$BUILD" == "configure" ]; then make check ; fi
+  - if [ "$BUILD" == "configure" ]; then ./configure ; fi
+  - if [ "$BUILD" == "configure" ]; then make all ; fi
+  - if [ "$BUILD" == "configure" ]; then make check ; fi

CHANGES

@@ -1,4 +1,82 @@
+arpack-ng - 3.6.0
+
+ [ Franck Houssen ]
+ * Add support for ISO_C_BINDING (Fortran 2003) for ARPACK, PARPACK (Fortran <-> C/C++).
+   ARPACK:  example of   C/Fortran binding can be found in the         TESTS/icb_arpack_c.c          file.
+   ARPACK:  example of C++/Fortran binding can be found in the         TESTS/icb_arpack_cpp.cpp      file.
+   PARPACK: example of   C/Fortran binding can be found in the PARPACK/TESTS/MPI/icb_parpack_c.c     file.
+   PARPACK: example of C++/Fortran binding can be found in the PARPACK/TESTS/MPI/icb_parpack_cpp.cpp file.
+   DEBUG:   add support for debug.
+   STAT:    add support for statistics (timers, nb operations, ...).
+ * Provide tarball generation using cmake (cpack).
+ * Provide find_package for (cmake) users to find arpack-ng.
+
+ [ Denis Davydov ]
+ * Rename pslamch to pslamch10 to avoid symbol collision with Scalapack 2.0.2 in MPI context.
+
+ [ Kyle Guinn ]
+ * Autoconf cleanup; move generated files to the build-aux subdirectory.
+
+ [ Marco Caliari ]
+ * Force the initial residual to be in the range of the operator OP in the standard case, too (Closes: #79).
+
+ [ Sylvestre Ledru ]
+ * Add coverage information to improve testing: https://coveralls.io/github/opencollab/arpack-ng
+
+ [ Darcy Beurle]
+ * Add C++11 interface through arpack.hpp and parpack.hpp
+ * Rewrite C++ examples / tests demonstrating new C++11 interface
+ * Pre-C++11 interface available through arpack.h and parpack.h
+
+ -- Sylvestre Ledru <sylvestre@debian.org> Mon, 30 Oct 2017 14:21:48 +0200
+
+arpack-ng - 3.5.0
+
+ [ Julien Schueller ]
+ * Improve cmake build system: disable C++ detection, set default build type.
+
+ [ Marco Atzeri]
+ * Use AC_PROG_FC instead of AC_PROG_F77 for proper inizialization
+   for the usage of AC_FC_LINE_LENGTH. Noted on Cygwin.
+
+ [ Denis Davydov ]
+ * Improve cmake build system: add make install and fix shared libraries.
+
+ [ Zhang Z ]
+ * fix usages of DLACPY to not alias inputs
+   (patch from https://software.intel.com/en-us/articles/how-to-resolve-arpack-issues-with-intel-mkl-110-update-3)
+
+ [ Iskakov Sergei ]
+ *  Fix possible deadlock when PARPACK call uses communicator with a larger
+    number of CPUs than previous call
+
+ [ Kyle Guinn ]
+ * Portability improvements to the autotools build system.
+ * Let cmake guess the default installation directories.  Can be
+   overridden by changing CMAKE_INSTALL_LIBDIR and CMAKE_INSTALL_BINDIR.
+ * Shared libraries built by cmake now have their SONAME set identical to
+   those built by autotools.
+
+ [ Marco Caliari ]
+ * Avoid purification stage in [d,s]neupd.f if it requires division
+    by zero (Closes: #58)
+
+ -- Sylvestre Ledru <sylvestre@debian.org>  Mon, 15 May 2017 14:21:48 +0200
+
+arpack-ng - 3.4.0
+
+  [ Milan Bouchet-Valat ]
+  * Allow adding suffixes to symbols and library names to build ILP64 version
+    based on ILP64 BLAS/LAPACK with suffixes. This avoids conflicts when loading
+    libraries with different integer sizes in the same program.
+
+  [ Martin Reuter ]
+  * Add the support of cmake build system
+
+ -- Sylvestre Ledru <sylvestre@debian.org>  Sat, 02 Jul 2016 21:51:52 +0200
+
 arpack-ng - 3.3.0
+
  [ Denis Davydov ]
  * Rename pdlamch to pdlamch10 to avoid symbol collision with Scalapack 2.0.2 in MPI context.
 
@@ -6,7 +84,7 @@ arpack-ng - 3.3.0
  * General improvements on the build system
  * libparpack links against libarpack (instead of doing a static link)
 
- [Guillaume Horel]
+ [ Guillaume Horel ]
  * reverts using {d,s}lahqr from lapack 2
  * use dlahqr from lapack 3 instead of dlaqrb (credit to Marco Caliari)
 
@@ -85,7 +163,7 @@ arpack-ng - 3.1.2
 
    * Wrong call to pdlamch was causing segfaults
      Thanks to Kyrre Sjøbæk for finding the bug and the fix.
-   * Get rid of the mpif.h occurences in the source code (Closes: #782)
+   * Get rid of the mpif.h occurrences in the source code (Closes: #782)
    * Compile also PARPACK / MPI example (Closes: #783)
    * Configure detected built-in LAPACK and BLAS, but refused to use them
      (Closes: #784)
@@ -110,7 +188,7 @@ arpack-ng - 3.1.0
    * Change the bug report from arpack@caam.rice.edu to
      http://forge.scilab.org/index.php/p/arpack-ng/issues/
    * Provide a M4 macro (detect_arpack_bug.m4) to check if the underlying
-     arpack is buggy (ie not arpack-ng). This allows developper applications
+     arpack is buggy (ie not arpack-ng). This allows developer applications
      to perform the check in their autotools build system (configure).
    * Fixed a lack of appropriate bounds check in DNAUP2. Thanks to Pauli Virtanen
      for the patch (Closes: #632)

CMakeLists.txt

@@ -0,0 +1,573 @@
+cmake_minimum_required(VERSION 2.8.6)
+
+if (NOT DEFINED CMAKE_BUILD_TYPE)
+   set (CMAKE_BUILD_TYPE Release CACHE STRING "Build type")
+endif ()
+
+project(arpack C Fortran)
+
+set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${PROJECT_SOURCE_DIR}/cmake)
+
+# Adopted from https://github.com/feymark/arpack.git
+
+if (POLICY CMP0042)
+  # enable MACOSX_RPATH by default
+  cmake_policy (SET CMP0042 NEW)
+endif ()
+
+option(MPI "Enable parallel support" OFF)
+option(ICB "Enable support for *[ae]upd_c with ISO_C_BINDING" OFF)
+#option(SYSTEM_BLAS "Use system BLAS" ON)
+#option(SYSTEM_LAPACK "Use system LAPACK" ON)
+option(EXAMPLES "Compile ARPACK examples" OFF)
+set(LIBSUFFIX ""
+    CACHE STRING "suffix to add to ARPACK libraries names")
+set(SYMBOLSUFFIX ""
+    CACHE STRING "suffix to add to ARPACK, BLAS and LAPACK function names")
+option(INTERFACE64 "use the 64-bit integer interface (ILP64) for ARPACK, BLAS and LAPACK")
+
+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(EXTRA_LDLAGS "${EXTRA_LDLAGS}-lgcov")
+endif()
+
+function(prefixlist list_name prefix)
+    set(${list_name}_TMP)
+    foreach(l ${${list_name}})
+      list(APPEND ${list_name}_TMP ${prefix}${l} )
+    endforeach()
+    set(${list_name} "${${list_name}_TMP}" PARENT_SCOPE)
+endfunction(prefixlist)
+
+function(examples list_name)
+    foreach(l ${${list_name}})
+        get_filename_component(lwe ${l} NAME_WE)
+        add_executable(${lwe} ${arpackexample_DIR}/${l} ${examples_EXTRA_SRCS})
+        target_link_libraries(${lwe} arpack ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES} ${EXTRA_LDLAGS})
+        add_test(NAME "${lwe}_ex" COMMAND ${lwe} WORKING_DIRECTORY ${arpackexample_DIR})
+    endforeach()
+endfunction(examples)
+
+function(pexamples list_name)
+    foreach(l ${${list_name}})
+        get_filename_component(lwe ${l} NAME_WE)
+        add_executable(${lwe} ${parpackexample_DIR}/${l} )
+        target_link_libraries(${lwe} parpack arpack  ${MPI_Fortran_LIBRARIES})
+    endforeach()
+endfunction(pexamples)
+
+# 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++.
+
+    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=gnu++11")
+
+    file(WRITE ${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/CMakeTmp/PROG_ICB.f90
+         "
+          PROGRAM PROG_ICB
+          USE iso_c_binding
+          IMPLICIT NONE
+          INTEGER(C_INT) :: a
+          a = 1
+          END PROGRAM PROG_ICB
+         ")
+    try_compile(COMPILE_ICB ${CMAKE_BINARY_DIR} ${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/CMakeTmp/PROG_ICB.f90)
+    if(NOT ${COMPILE_ICB})
+      message(FATAL_ERROR "-- Fortran compiler does not support iso_c_binding.")
+    else()
+      message("-- Fortran compiler does support iso_c_binding.")
+    endif()
+endif()
+
+#if (SYSTEM_BLAS)
+    find_package(BLAS REQUIRED)
+#endif()
+#if (BLAS_LIBRARIES)
+#    set(SYSTEM_BLAS ON)
+#endif()
+
+if (MPI)
+    include(FindMPI)
+    find_package(MPI REQUIRED)
+    include_directories(${MPI_Fortran_INCLUDE_PATH})
+    set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} ${MPI_Fortran_COMPILE_FLAG}")
+
+    if(ICB)
+        include_directories(${MPI_C_INCLUDE_PATH})
+        check_symbol_exists(MPI_Comm_c2f "${MPI_C_INCLUDE_PATH}/mpi.h" MPI_Comm_c2f_FOUND)
+        if(NOT ${MPI_Comm_c2f_FOUND})
+          message(FATAL_ERROR "symbol MPI_Comm_c2f does not exist")
+        endif()
+    endif()
+endif()
+
+#if (SYSTEM_LAPACK)
+    find_package(LAPACK REQUIRED)
+#endif()
+#if (BLAS_LIBRARIES)
+#    set(SYSTEM_LAPACK ON)
+#endif()
+
+if (MPI)
+    set(parpackutil_DIR ${arpack_SOURCE_DIR}/PARPACK/UTIL/)
+    set(parpacksrc_DIR ${arpack_SOURCE_DIR}/PARPACK/SRC/)
+endif()
+
+#if (NOT SYSTEM_BLAS)
+#    file(GLOB arpackblas_STAT_SRCS blas/*.f)
+#endif()
+#if (NOT SYSTEM_LAPACK)
+#    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)
+
+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)
+endif()
+
+set(arpackutil_STAT_SRCS
+    ${arpack_SOURCE_DIR}/UTIL/icnteq.f
+    ${arpack_SOURCE_DIR}/UTIL/icopy.f
+    ${arpack_SOURCE_DIR}/UTIL/iset.f
+    ${arpack_SOURCE_DIR}/UTIL/iswap.f
+    ${arpack_SOURCE_DIR}/UTIL/ivout.f
+    ${arpack_SOURCE_DIR}/UTIL/second_NONE.f
+    ${arpack_SOURCE_DIR}/UTIL/svout.f
+    ${arpack_SOURCE_DIR}/UTIL/smout.f
+    ${arpack_SOURCE_DIR}/UTIL/dvout.f
+    ${arpack_SOURCE_DIR}/UTIL/dmout.f
+    ${arpack_SOURCE_DIR}/UTIL/cvout.f
+    ${arpack_SOURCE_DIR}/UTIL/cmout.f
+    ${arpack_SOURCE_DIR}/UTIL/zvout.f
+    ${arpack_SOURCE_DIR}/UTIL/zmout.f )
+
+
+if (MPI)
+    file(GLOB parpackutil_STAT_SRCS PARPACK/UTIL/MPI/*.f)
+    file(GLOB parpacksrc_STAT_SRCS 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)
+  install(FILES arpack.hpp DESTINATION include)
+  if (MPI)
+    install(FILES parpack.h DESTINATION include)
+    install(FILES parpack.hpp DESTINATION include)
+  endif()
+endif()
+
+install(FILES debug.h DESTINATION include)
+if(ICB)
+  install(FILES debug_c.h DESTINATION include)
+  install(FILES debug_c.hpp DESTINATION include)
+endif()
+
+install(FILES stat.h DESTINATION include)
+if(ICB)
+  install(FILES stat_c.h DESTINATION include)
+  install(FILES stat_c.hpp DESTINATION include)
+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 SOVERSION 2)
+
+if (MPI)
+    # use -DBUILD_SHARED_LIBS=ON|OFF to control static/shared
+    add_library(parpack ${parpacksrc_STAT_SRCS} ${parpackutil_STAT_SRCS} ${parpacksrc_ICB})
+
+    if (ICB)
+        target_include_directories(parpack PUBLIC ${MPI_C_INCLUDE_DIRS})
+        target_include_directories(parpack PUBLIC ${MPI_CXX_INCLUDE_DIRS})
+        target_link_libraries(parpack ${MPI_C_LIBRARIES})
+        target_link_libraries(parpack ${MPI_CXX_LIBRARIES})
+    endif()
+    target_include_directories(parpack PUBLIC ${MPI_Fortran_INCLUDE_DIRS})
+    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 SOVERSION 2)
+endif ()
+
+if (INTERFACE64)
+    if ("${CMAKE_Fortran_COMPILER_ID}" MATCHES "GNU")
+      set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -fdefault-integer-8")
+    elseif ("${CMAKE_Fortran_COMPILER_ID}" MATCHES "Intel")
+      set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -i8")
+    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")
+endif ()
+
+if (SYMBOLSUFFIX)
+    if ("${CMAKE_Fortran_COMPILER_ID}" MATCHES "GNU")
+      set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -cpp -ffixed-line-length-none")
+    elseif ("${CMAKE_Fortran_COMPILER_ID}" MATCHES "Intel")
+      set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -fpp -extend-source")
+    else ()
+      message(WARNING "build script does not know how to enable your Fortran compiler's preprocessor and support for lines longer than 72 characters: set them manually via FFLAGS.")
+    endif ()
+
+    set(SCALARFUNS naitr napps naup2 naupd nconv neigh neupd ngets statn saitr sapps saup2 saupd sconv seigt seupd sgets stats getv0 sortc sortr sesrt stqrb)
+    set(COMPLEXFUNS, naitr napps naup2 naupd neigh neupd ngets statn getv0 sortc)
+
+    set(BLASFUNS1 axpy copy gemv geqr2 lacpy lae2 lahqr lanhs larnv lartg lascl laset lasrt scal trevc trmm trsen gbmv gbtrf gbtrs gttrf gttrs pttrf pttrs)
+    set(BLASFUNS2 dot ger labad laev2 lamch lanst lanv2 lapy2 larf larfg lasr nrm2 orm2r rot steqr swap)
+    set(BLASFUNS3 dotc geru unm2r)
+    set(BLASFUNS4 COPY LABAD LAMCH LANHS LANV2 LARFG ROT GEMV)
+    set(BLASFUNS5 scnrm2 dznrm2 csscal zdscal)
+
+    foreach (f IN LISTS SCALARFUNS BLASFUNS1 BLASFUNS2)
+      set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -Ds${f}=s${f}${SYMBOLSUFFIX} -Dd${f}=d${f}${SYMBOLSUFFIX}")
+    endforeach ()
+
+    foreach (f IN LISTS COMPLEXFUNS BLASFUNS1 BLASFUNS3)
+      set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -Dc${f}=c${f}${SYMBOLSUFFIX} -Dz${f}=z${f}${SYMBOLSUFFIX}")
+    endforeach ()
+
+    foreach (f IN LISTS BLASFUNS4)
+      set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -DS${f}=S${f}${SYMBOLSUFFIX} -DD${f}=D${f}${SYMBOLSUFFIX}")
+    endforeach ()
+
+    foreach (f IN LISTS BLASFUNS5)
+      set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -D${f}=${f}${SYMBOLSUFFIX}")
+    endforeach ()
+endif ()
+
+set(CFUNS sgemm snaupd sneupd dnaupd dneupd cheev)
+foreach (f IN LISTS CFUNS)
+  set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -D${f}=${f}${SYMBOLSUFFIX}")
+  list(APPEND CFUNS_SUFFIXED ${f}${SYMBOLSUFFIX})
+endforeach ()
+
+include(FortranCInterface)
+FortranCInterface_HEADER(FCMangle.h SYMBOLS ${CFUNS_SUFFIXED})
+
+set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -DINCLUDE_FCMANGLE")
+set(CMAKE_INCLUDE_CURRENT_DIR ON)
+
+############################
+# EXAMPLES
+############################
+if (EXAMPLES)
+    ############################
+    # EXAMPLES/BAND
+    ############################
+    set(arpackexample_DIR ${arpack_SOURCE_DIR}/EXAMPLES/BAND/)
+    set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/Examples/band/)
+
+    set(examples_EXTRA_SRCS ${arpackexample_DIR}/cnband.f)
+    set(examples_STAT_SRCS
+        cnbdr1.f
+        cnbdr2.f
+        cnbdr3.f
+        cnbdr4.f)
+    examples(examples_STAT_SRCS)
+
+    set(examples_EXTRA_SRCS ${arpackexample_DIR}/dnband.f)
+    set(examples_STAT_SRCS
+        dnbdr1.f
+        dnbdr2.f
+        dnbdr3.f
+        dnbdr4.f
+        dnbdr5.f
+        dnbdr6.f)
+    examples(examples_STAT_SRCS)
+
+    set(examples_EXTRA_SRCS ${arpackexample_DIR}/dsband.f)
+    set(examples_STAT_SRCS
+        dsbdr1.f
+        dsbdr2.f
+        dsbdr3.f
+        dsbdr4.f
+        dsbdr5.f
+        dsbdr6.f)
+    examples(examples_STAT_SRCS)
+
+    set(examples_EXTRA_SRCS ${arpackexample_DIR}/snband.f)
+    set(examples_STAT_SRCS
+        snbdr1.f
+        snbdr2.f
+        snbdr3.f
+        snbdr4.f
+        snbdr5.f
+        snbdr6.f)
+    examples(examples_STAT_SRCS)
+
+    set(examples_EXTRA_SRCS ${arpackexample_DIR}/ssband.f)
+    set(examples_STAT_SRCS
+        ssbdr1.f
+        ssbdr2.f
+        ssbdr3.f
+        ssbdr4.f
+        ssbdr5.f
+        ssbdr6.f)
+    examples(examples_STAT_SRCS)
+
+    set(examples_EXTRA_SRCS ${arpackexample_DIR}/znband.f)
+    set(examples_STAT_SRCS
+        znbdr1.f
+        znbdr2.f
+        znbdr3.f
+        znbdr4.f)
+    examples(examples_STAT_SRCS)
+    set(examples_EXTRA_SRCS)
+
+    ############################
+    # EXAMPLES/COMPLEX
+    ############################
+    set(arpackexample_DIR ${arpack_SOURCE_DIR}/EXAMPLES/COMPLEX/)
+    set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/Examples/complex/)
+
+    set(examples_STAT_SRCS
+        cndrv1.f
+        cndrv2.f
+        cndrv3.f
+        cndrv4.f
+        zndrv1.f
+        zndrv2.f
+        zndrv3.f
+        zndrv4.f)
+    examples(examples_STAT_SRCS)
+
+    ############################
+    # examples/nonsym
+    ############################
+    set(arpackexample_DIR ${arpack_SOURCE_DIR}/EXAMPLES/NONSYM/)
+    set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/Examples/nonsym/)
+
+    set(examples_STAT_SRCS
+        dndrv1.f
+        dndrv2.f
+        dndrv3.f
+        dndrv4.f
+        dndrv5.f
+        dndrv6.f
+        sndrv1.f
+        sndrv2.f
+        sndrv3.f
+        sndrv4.f
+        sndrv5.f
+        sndrv6.f)
+    examples(examples_STAT_SRCS)
+
+    ############################
+    # examples/SIMPLE
+    ############################
+    set(arpackexample_DIR ${arpack_SOURCE_DIR}/EXAMPLES/SIMPLE/)
+    set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/Examples/simple/)
+
+    set(examples_STAT_SRCS
+        cnsimp.f
+        dnsimp.f
+        dssimp.f
+        snsimp.f
+        sssimp.f
+        znsimp.f)
+
+    examples(examples_STAT_SRCS)
+
+    ############################
+    # examples/svd
+    ############################
+    set(arpackexample_DIR ${arpack_SOURCE_DIR}/EXAMPLES/SVD/)
+    set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/Examples/svd/)
+
+    set(examples_STAT_SRCS
+        ssvd.f)
+
+    examples(examples_STAT_SRCS)
+
+    ############################
+    # examples/sym
+    ############################
+    set(arpackexample_DIR ${arpack_SOURCE_DIR}/EXAMPLES/SYM/)
+    set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/Examples/sym/)
+
+    set(examples_STAT_SRCS
+        dsdrv1.f
+        dsdrv2.f
+        dsdrv3.f
+        dsdrv4.f
+        dsdrv5.f
+        dsdrv6.f
+        ssdrv1.f
+        ssdrv2.f
+        ssdrv3.f
+        ssdrv4.f
+        ssdrv5.f
+        ssdrv6.f)
+
+    examples(examples_STAT_SRCS)
+
+    ############################
+    # parpack/examples
+    ############################
+    if (MPI)
+        set(parpackexample_DIR ${arpack_SOURCE_DIR}/PARPACK/EXAMPLES/MPI/)
+        set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/Examples/parpack)
+
+        set(pexamples_STAT_SRCS
+            pcndrv1.f
+            pdndrv1.f
+            pdndrv3.f
+            pdsdrv1.f
+            psndrv1.f
+            psndrv3.f
+            pssdrv1.f
+            pzndrv1.f)
+
+        pexamples(pexamples_STAT_SRCS)
+    endif()
+endif()
+
+
+############################
+# TEST
+############################
+
+enable_testing()
+
+set(CMAKE_CTEST_COMMAND ctest -V)
+
+set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/Tests)
+
+add_executable(dnsimp_test TESTS/dnsimp.f TESTS/mmio.f TESTS/debug.h)
+set_target_properties( dnsimp_test PROPERTIES OUTPUT_NAME  dnsimp )
+target_link_libraries(dnsimp_test arpack ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES} ${EXTRA_LDLAGS})
+add_custom_command(TARGET dnsimp_test POST_BUILD
+  COMMAND ${CMAKE_COMMAND} -E copy ${CMAKE_CURRENT_SOURCE_DIR}/TESTS/testA.mtx testA.mtx
+)
+add_test(dnsimp_tst Tests/dnsimp)
+
+add_executable(bug_1315_single TESTS/bug_1315_single.c)
+target_link_libraries(bug_1315_single arpack ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES} ${EXTRA_LDLAGS})
+add_test(bug_1315_single_tst Tests/bug_1315_single)
+
+add_executable(bug_1315_double TESTS/bug_1315_double.c)
+target_link_libraries(bug_1315_double arpack ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES} ${EXTRA_LDLAGS})
+add_test(bug_1315_double_tst Tests/bug_1315_double)
+
+add_executable(bug_1323 TESTS/bug_1323.f)
+target_link_libraries(bug_1323 arpack ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES} ${EXTRA_LDLAGS})
+add_test(bug_1323_tst Tests/bug_1323)
+
+add_executable(bug_58_double TESTS/bug_58_double.f)
+target_link_libraries(bug_58_double arpack ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES} ${EXTRA_LDLAGS})
+add_test(bug_58_double_tst Tests/bug_58_double)
+
+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)
+
+if(MPI)
+  add_executable(issue46 PARPACK/TESTS/MPI/issue46.f)
+  target_link_libraries(issue46 parpack arpack ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES} ${EXTRA_LDLAGS})
+  add_test(issue46_tst Tests/issue46)
+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_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_link_libraries(icb_arpack_cpp arpack ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES} ${EXTRA_LDLAGS})
+  add_test(icb_arpack_cpp_tst Tests/icb_arpack_cpp)
+
+  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_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_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()
+endif()
+
+############################
+# install
+############################
+# '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}
+    RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR})
+
+if (MPI)
+  install(TARGETS parpack
+      ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR}
+      LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}
+      RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR})
+endif ()
+
+# Provide find_package for arpack-ng to users.
+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_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)
+install(FILES "${PROJECT_BINARY_DIR}/arpack-ng-config-version.cmake" DESTINATION "${CMAKE_INSTALL_PREFIX}/lib/cmake")
+
+# Packaging: ease arpack-ng distribution (precompiled binaries and sources tarballs).
+set(CPACK_VERSION_MAJOR "${arpack_ng_MAJOR_VERSION}")
+set(CPACK_VERSION_MINOR "${arpack_ng_MINOR_VERSION}")
+set(CPACK_VERSION_PATCH "${arpack_ng_PATCH_VERSION}")
+set(CPACK_SOURCE_IGNORE_FILES "/BUILD/" "/Build/" "/build/" "/local/") # Do not ship (eventual) build or install directories in tarballs.
+set(CPACK_SOURCE_PACKAGE_FILE_NAME "arpack-ng-${CPACK_VERSION_MAJOR}.${CPACK_VERSION_MINOR}.${CPACK_VERSION_PATCH}")
+include(CPack)
+
+
+# We don't want this to run on every build.
+option(COVERALLS "Generate coveralls data" OFF)
+
+
+if (COVERALLS)
+    set(COVERAGE_SRCS awesome.c code.c files.c)
+
+    # Create the coveralls target.
+    coveralls_setup(
+        "${arpackutil_STAT_SRCS} ${arpacksrc_STAT_SRCS} ${arpacksrc_ICB}" # The source files.
+        ON                 # If we should upload.
+        "${PROJECT_SOURCE_DIR}/cmake/") # (Optional) Alternate project cmake module path.
+endif()
+

DOCUMENTS/debug.doc

@@ -219,7 +219,7 @@ c     %---------------------------------%
               of the Hessenberg matrix and their the corresponding 
               error bounds, KEV, NP.
 
- mnapps > 0: *Print information about where deflation occured.
+ mnapps > 0: *Print information about where deflation occurred.
  mnapps > 1: *Print sigmak, betak, order of the final Hessenberg matrix,
               and the final compressed upper Hessenberg matrix.
  mnapps > 2: *Print implicit application of shift number, real and imaginary 
@@ -303,7 +303,7 @@ c     %---------------------------------%
               of the Hessenberg matrix and their the corresponding 
               error bounds, KEV, NP.
 
- mcapps > 0: *Print information about where deflation occured.
+ mcapps > 0: *Print information about where deflation occurred.
  mcapps > 1: *Print sigmak, betak, order of the final Hessenberg matrix,
               and the final compressed upper Hessenberg matrix.
  mcapps > 2: *Print implicit application of shift number, the shift.

EXAMPLES/BAND/Makefile.am

@@ -1,4 +1,4 @@
-LDADD = $(top_builddir)/libarpack.la $(BLAS_LIBS) $(LAPACK_LIBS)
+LDADD = $(top_builddir)/SRC/libarpack$(LIBSUFFIX).la $(LAPACK_LIBS) $(BLAS_LIBS)
 
 SSBDR = ssbdr1 ssbdr2 ssbdr3 ssbdr4 ssbdr5 ssbdr6
 DSBDR = dsbdr1 dsbdr2 dsbdr3 dsbdr4 dsbdr5 dsbdr6
@@ -13,7 +13,7 @@ TESTS = $(check_PROGRAMS)
 
 EXTRA_DIST = README
 
-# Simple symetric problem using BAND solver (single precision)
+# Simple symmetric problem using BAND solver (single precision)
 ssbdr1_SOURCES = ssbdr1.f ssband.f
 ssbdr2_SOURCES = ssbdr2.f ssband.f
 ssbdr3_SOURCES = ssbdr3.f ssband.f

EXAMPLES/BAND/cnband.f

@@ -158,7 +158,7 @@ c          Stopping criteria: the relative accuracy of the Ritz value
 c          is considered acceptable if BOUNDS(I) .LE. TOL*ABS(RITZ(I))
 c          where ABS(RITZ(I)) is the magnitude when RITZ(I) is complex.
 c          DEFAULT = slamch('EPS')  (machine precision as computed
-c                    by the LAPACK auxilliary subroutine slamch).
+c                    by the LAPACK auxiliary subroutine slamch).
 c
 c  RESID   Complex  array of length N.  (INPUT/OUTPUT)
 c          On INPUT:

EXAMPLES/BAND/znband.f

@@ -158,7 +158,7 @@ c          Stopping criteria: the relative accuracy of the Ritz value
 c          is considered acceptable if BOUNDS(I) .LE. TOL*ABS(RITZ(I))
 c          where ABS(RITZ(I)) is the magnitude when RITZ(I) is complex.
 c          DEFAULT = dlamch ('EPS')  (machine precision as computed
-c                    by the LAPACK auxilliary subroutine dlamch ).
+c                    by the LAPACK auxiliary subroutine dlamch ).
 c
 c  RESID   Complex*16  array of length N.  (INPUT/OUTPUT)
 c          On INPUT:

EXAMPLES/COMPLEX/Makefile.am

@@ -1,4 +1,4 @@
-LDADD = $(top_builddir)/libarpack.la $(BLAS_LIBS) $(LAPACK_LIBS)
+LDADD = $(top_builddir)/SRC/libarpack$(LIBSUFFIX).la $(LAPACK_LIBS) $(BLAS_LIBS)
 
 CNDRV = cndrv1 cndrv2 cndrv3 cndrv4
 ZNDRV = zndrv1 zndrv2 zndrv3 zndrv4

EXAMPLES/NONSYM/Makefile.am

@@ -1,4 +1,4 @@
-LDADD = $(top_builddir)/libarpack.la $(BLAS_LIBS) $(LAPACK_LIBS)
+LDADD = $(top_builddir)/SRC/libarpack$(LIBSUFFIX).la $(LAPACK_LIBS) $(BLAS_LIBS)
 
 SNDRV = sndrv1 sndrv2 sndrv3 sndrv4 sndrv5 sndrv6
 DNDRV = dndrv1 dndrv2 dndrv3 dndrv4 dndrv5 dndrv6

EXAMPLES/README.CALLING-ARPACK-FROM-C-OR-CPP

@@ -0,0 +1,2 @@
+In ../TESTS, the file icb_arpack_c.c     is an example of how to call arpack from C.
+In ../TESTS, the file icb_arpack_cpp.cpp is an example of how to call arpack from C++.

EXAMPLES/SIMPLE/Makefile.am

@@ -1,4 +1,4 @@
-LDADD = $(top_builddir)/libarpack.la $(BLAS_LIBS) $(LAPACK_LIBS)
+LDADD = $(top_builddir)/SRC/libarpack$(LIBSUFFIX).la $(LAPACK_LIBS) $(BLAS_LIBS)
 
 SIMPLE = sssimp dssimp snsimp dnsimp cnsimp znsimp
 

EXAMPLES/SVD/Makefile.am

@@ -1,4 +1,4 @@
-LDADD = $(top_builddir)/libarpack.la $(BLAS_LIBS) $(LAPACK_LIBS)
+LDADD = $(top_builddir)/SRC/libarpack$(LIBSUFFIX).la $(LAPACK_LIBS) $(BLAS_LIBS)
 
 SVD = ssvd dsvd
 

EXAMPLES/SYM/Makefile.am

@@ -1,4 +1,4 @@
-LDADD = $(top_builddir)/libarpack.la $(BLAS_LIBS) $(LAPACK_LIBS)
+LDADD = $(top_builddir)/SRC/libarpack$(LIBSUFFIX).la $(LAPACK_LIBS) $(BLAS_LIBS)
 
 SSDRV = ssdrv1 ssdrv2 ssdrv3 ssdrv4 ssdrv5 ssdrv6
 DSDRV = dsdrv1 dsdrv2 dsdrv3 dsdrv4 dsdrv5 dsdrv6

Makefile.am

@@ -1,32 +1,34 @@
-SUBDIRS = UTIL SRC . TESTS EXAMPLES
-
-if MPI
-SUBDIRS += PARPACK
-AM_DISTCHECK_CONFIGURE_FLAGS = --enable-mpi
-endif
-
-lib_LTLIBRARIES = libarpack.la
-
 ACLOCAL_AMFLAGS = -I m4
 
-libarpack_la_SOURCES =
+AM_DISTCHECK_CONFIGURE_FLAGS =
+if MPI
+AM_DISTCHECK_CONFIGURE_FLAGS += --enable-mpi
+endif
+include_HEADERS = debug.h stat.h
+if ICB
+AM_DISTCHECK_CONFIGURE_FLAGS += --enable-icb
+include_HEADERS += debug_c.h debug_c.hpp
+include_HEADERS += stat_c.h stat_c.hpp
+include_HEADERS += arpack.h arpack.hpp
+if MPI
+include_HEADERS += parpack.h parpack.hpp
+endif
+endif
 
-# Force libarpack to be linked with the Fortran compiler.
-# The file dummy.f does not need to exist in the source tree.
-nodist_EXTRA_libarpack_la_SOURCES = dummy.f
+SUBDIRS = UTIL SRC . TESTS EXAMPLES PARPACK
 
-libarpack_la_LDFLAGS = -no-undefined -version-info 2:0
+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
 
-libarpack_la_LIBADD = \
-	$(top_builddir)/SRC/libarpacksrc.la \
-	$(top_builddir)/UTIL/libarpackutil.la \
-	$(LAPACK_LIBS) $(BLAS_LIBS)
+pkgconfig_DATA = arpack@LIBSUFFIX@.pc
 
-EXTRA_DIST = README PARPACK_CHANGES CHANGES DOCUMENTS VISUAL_STUDIO \
-detect_arpack_bug.m4
-
-# Pkgconfig directory
-pkgconfigdir = $(libdir)/pkgconfig
-
-# Files to install in Pkgconfig directory
-pkgconfig_DATA = arpack.pc
+DISTCLEANFILES = $(pkgconfig_DATA)

PARPACK/EXAMPLES/BLACS/Makefile.am

@@ -1,43 +1,38 @@
-#LDADD = $(top_builddir)/PARPACK/SRC/BLACS/libparpacksrc.la $(BLAS_LIBS) $(LAPACK_LIBS)
+LDADD = $(top_builddir)/PARPACK/SRC/BLACS/libparpack$(LIBSUFFIX).la $(LAPACK_LIBS) $(BLAS_LIBS)
 
-#SNDRV = psndrv1 psndrv3
-#DNDRV = pdndrv1 pdndrv3
-#SSDRV = pssdrv1
-#DSDRV = pdsdrv1
-#CNDRV = pcndrv1
-#ZNDRV = pzndrv1
-#NTEST = psntest1 pdntest1
+SNDRV = psndrv1 psndrv3
+DNDRV = pdndrv1 pdndrv3
+SSDRV = pssdrv1
+DSDRV = pdsdrv1
+CNDRV = pcndrv1
+ZNDRV = pzndrv1
+NTEST = psntest1 pdntest1
 
-#check_PROGRAMS = $(SNDRV) $(DNDRV) $(SSDRV) $(DSDRV) $(CNDRV) $(ZNDRV) $(NTEST)
-#TESTS = $(check_PROGRAMS)
+check_PROGRAMS = $(SNDRV) $(DNDRV) $(SSDRV) $(DSDRV) $(CNDRV) $(ZNDRV) $(NTEST)
+TESTS = $(check_PROGRAMS)
 
-EXTRA_DIST = debug.h stat.h \
-             psndrv1.f psndrv3.f \
-             pdndrv1.f pdndrv3.f \
-             pssdrv1.f pdsdrv1.f \
-             pcndrv1.f pzndrv1.f \
-             psntest1.f pdntest1.f
+EXTRA_DIST = debug.h stat.h
 
 # Simple nonsymmetric problem using single precision
-#psndrv1_SOURCES = psndrv1.f
-#psndrv3_SOURCES = psndrv3.f
+psndrv1_SOURCES = psndrv1.f
+psndrv3_SOURCES = psndrv3.f
 
 # Simple nonsymmetric problem using double precision
-#pdndrv1_SOURCES = pdndrv1.f
-#pdndrv3_SOURCES = pdndrv3.f
+pdndrv1_SOURCES = pdndrv1.f
+pdndrv3_SOURCES = pdndrv3.f
 
 # Simple symmetric problem using single precision
-#pssdrv1_SOURCES = pssdrv1.f
+pssdrv1_SOURCES = pssdrv1.f
 
 # Simple symmetric problem using double precision
-#pdsdrv1_SOURCES = pdsdrv1.f
+pdsdrv1_SOURCES = pdsdrv1.f
 
 # Complex problem using single complex
-#pcndrv1_SOURCES = pcndrv1.f
+pcndrv1_SOURCES = pcndrv1.f
 
 # Complex problem using double complex
-#pzndrv1_SOURCES = pzndrv1.f
+pzndrv1_SOURCES = pzndrv1.f
 
 # Test routines for timing
-#psntest1_SOURCES = psntest1.f
-#pdntest1_SOURCES = pdntest1.f
+psntest1_SOURCES = psntest1.f
+pdntest1_SOURCES = pdntest1.f

PARPACK/EXAMPLES/MPI/Makefile.am

@@ -1,5 +1,5 @@
 F77 = $(MPIF77)
-LDADD = $(top_builddir)/PARPACK/libparpack.la $(BLAS_LIBS) $(LAPACK_LIBS)
+LDADD = $(top_builddir)/PARPACK/SRC/MPI/libparpack$(LIBSUFFIX).la $(LAPACK_LIBS) $(BLAS_LIBS)
 
 SNDRV = psndrv1 psndrv3
 DNDRV = pdndrv1 pdndrv3

PARPACK/EXAMPLES/MPI/README.CALLING-PARPACK-FROM-C-OR-CPP

@@ -0,0 +1,2 @@
+In PARPACK/TESTS/MPI, the file icb_parpack_c.c     is an example of how to call parpack from C.
+In PARPACK/TESTS/MPI, the file icb_parpack_cpp.cpp is an example of how to call parpack from C++.

PARPACK/Makefile.am

@@ -1,18 +1,7 @@
-SUBDIRS = UTIL SRC . EXAMPLES/MPI EXAMPLES/BLACS
+if MPI
+SUBDIRS = UTIL/MPI SRC/MPI . TESTS/MPI EXAMPLES/MPI
+endif
 
-lib_LTLIBRARIES = libparpack.la
-
-libparpack_la_SOURCES =
-
-# Force libparpack to be linked with the MPI Fortran compiler.
-# The file dummy.f does not need to exist in the source tree.
-F77 = $(MPIF77)
-nodist_EXTRA_libparpack_la_SOURCES = dummy.f
-
-libparpack_la_LDFLAGS = -no-undefined -version-info 2:0
-
-libparpack_la_LIBADD = \
-	$(top_builddir)/PARPACK/SRC/MPI/libparpacksrcmpi.la \
-	$(top_builddir)/PARPACK/UTIL/MPI/libparpackutilmpi.la \
-	$(top_builddir)/libarpack.la \
-	$(LAPACK_LIBS) $(BLAS_LIBS) $(MPILIBS)
+if BLACS
+SUBDIRS = UTIL/BLACS SRC/BLACS . EXAMPLES/BLACS
+endif

PARPACK/SRC/BLACS/Makefile.am

@@ -14,5 +14,10 @@ ZSRC = pznaitr.f pznapps.f pznaup2.f pznaupd.f pzneigh.f pzneupd.f pzngets.f \
 
 EXTRA_DIST = debug.h stat.h
 
-noinst_LTLIBRARIES = libparpacksrc.la
-libparpacksrc_la_SOURCES = $(SSRC) $(DSRC) $(CSRC) $(ZSRC)
+lib_LTLIBRARIES = libparpack@LIBSUFFIX@.la
+libparpack@LIBSUFFIX@_la_SOURCES = $(SSRC) $(DSRC) $(CSRC) $(ZSRC)
+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

PARPACK/SRC/BLACS/pcnapps.f

@@ -225,7 +225,7 @@ c
       cabs1( cdum ) = abs( real( cdum ) ) + abs( aimag( cdum ) )
 c
 c     %----------------%
-c     | Data statments |
+c     | Data statements |
 c     %----------------%
 c
       data       first / .true. /

PARPACK/SRC/BLACS/pcnaup2.f

@@ -28,7 +28,7 @@ c          products (involving the operator OP) per Arnoldi iteration.
 c          The logic for adjusting is contained within the current
 c          subroutine.
 c          If ISHIFT=0, NP is the number of shifts the user needs
-c          to provide via reverse comunication. 0 < NP < NCV-NEV.
+c          to provide via reverse communication. 0 < NP < NCV-NEV.
 c          NP may be less than NCV-NEV since a leading block of the current
 c          upper Hessenberg matrix has split off and contains "unwanted"
 c          Ritz values.

PARPACK/SRC/BLACS/pcnaupd.f

@@ -258,7 +258,7 @@ c          = -7: Length of private work array is not sufficient.
 c          = -8: Error return from LAPACK eigenvalue calculation;
 c          = -9: Starting vector is zero.
 c          = -10: IPARAM(7) must be 1,2,3.
-c          = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatable.
+c          = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatible.
 c          = -12: IPARAM(1) must be equal to 0 or 1.
 c          = -9999: Could not build an Arnoldi factorization.
 c                   User input error highly likely.  Please

PARPACK/SRC/BLACS/pdlamch.f

@@ -1,6 +1,6 @@
       DOUBLE PRECISION   FUNCTION PDLAMCH( ICTXT, CMACH )
 *
-*  -- ScaLAPACK auxilliary routine (version 1.0) --
+*  -- ScaLAPACK auxiliary routine (version 1.0) --
 *     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
 *     and University of California, Berkeley.
 *     February 28, 1995

PARPACK/SRC/BLACS/pdnaitr.f

@@ -398,7 +398,7 @@ c
 c 
 c        %---------------------------------------------------%
 c        | STEP 1: Check if the B norm of j-th residual      |
-c        | vector is zero. Equivalent to determing whether   |
+c        | vector is zero. Equivalent to determining whether |
 c        | an exact j-step Arnoldi factorization is present. |
 c        %---------------------------------------------------%
 c

PARPACK/SRC/BLACS/pdnapps.f

@@ -34,7 +34,7 @@ c
 c  KEV     Integer.  (INPUT/OUTPUT)
 c          KEV+NP is the size of the input matrix H.
 c          KEV is the size of the updated matrix HNEW.  KEV is only 
-c          updated on ouput when fewer than NP shifts are applied in
+c          updated on output when fewer than NP shifts are applied in
 c          order to keep the conjugate pair together.
 c
 c  NP      Integer.  (INPUT)
@@ -215,7 +215,7 @@ c
       intrinsic  abs, max, min
 c
 c     %----------------%
-c     | Data statments |
+c     | Data statements |
 c     %----------------%
 c
       data       first / .true. /
@@ -395,7 +395,7 @@ c
             do 80 i = istart, iend-1
 c
 c              %-----------------------------------------------------%
-c              | Contruct the plane rotation G to zero out the bulge |
+c              | Construct the plane rotation G to zero out the bulge |
 c              %-----------------------------------------------------%
 c
                call dlartg (f, g, c, s, r)

PARPACK/SRC/BLACS/pdnaup2.f

@@ -28,7 +28,7 @@ c          products (involving the operator OP) per Arnoldi iteration.
 c          The logic for adjusting is contained within the current
 c          subroutine.
 c          If ISHIFT=0, NP is the number of shifts the user needs
-c          to provide via reverse comunication. 0 < NP < NCV-NEV.
+c          to provide via reverse communication. 0 < NP < NCV-NEV.
 c          NP may be less than NCV-NEV for two reasons. The first, is
 c          to keep complex conjugate pairs of "wanted" Ritz values 
 c          together. The second, is that a leading block of the current
@@ -715,7 +715,7 @@ c
          if (ishift .eq. 0) then
 c
 c           %-------------------------------------------------------%
-c           | User specified shifts: reverse comminucation to       |
+c           | User specified shifts: reverse communication to       |
 c           | compute the shifts. They are returned in the first    |
 c           | 2*NP locations of WORKL.                              |
 c           %-------------------------------------------------------%

PARPACK/SRC/BLACS/pdnaupd.f

@@ -276,7 +276,7 @@ c          = -7: Length of private work array is not sufficient.
 c          = -8: Error return from LAPACK eigenvalue calculation;
 c          = -9: Starting vector is zero.
 c          = -10: IPARAM(7) must be 1,2,3,4.
-c          = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatable.
+c          = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatible.
 c          = -12: IPARAM(1) must be equal to 0 or 1.
 c          = -9999: Could not build an Arnoldi factorization.
 c                   IPARAM(5) returns the size of the current Arnoldi

PARPACK/SRC/BLACS/pdsaitr.f

@@ -396,8 +396,8 @@ c
          end if
 c
 c        %---------------------------------------------------------%
-c        | Check for exact zero. Equivalent to determing whether a |
-c        | j-step Arnoldi factorization is present.                |
+c        | Check for exact zero. Equivalent to determining whether |
+c        | a j-step Arnoldi factorization is present.              |
 c        %---------------------------------------------------------%
 c
          if (rnorm .gt. zero) go to 40

PARPACK/SRC/BLACS/pdsapps.f

@@ -204,7 +204,7 @@ c
       intrinsic  abs
 c
 c     %----------------%
-c     | Data statments |
+c     | Data statements |
 c     %----------------%
 c
       data       first / .true. /
@@ -275,7 +275,7 @@ c
                   call pivout (comm, logfil, 1, i, ndigit, 
      &                 '_sapps: deflation at row/column no.')
                   call pivout (comm, logfil, 1, jj, ndigit, 
-     &                 '_sapps: occured before shift number.')
+     &                 '_sapps: occurred before shift number.')
                   call pdvout (comm, logfil, 1, h(i+1,1), ndigit, 
      &                 '_sapps: the corresponding off diagonal element')
                end if

PARPACK/SRC/BLACS/pdsaup2.f

@@ -29,7 +29,7 @@ c          products (involving the operator OP) per Arnoldi iteration.
 c          The logic for adjusting is contained within the current
 c          subroutine.
 c          If ISHIFT=0, NP is the number of shifts the user needs
-c          to provide via reverse comunication. 0 < NP < NCV-NEV.
+c          to provide via reverse communication. 0 < NP < NCV-NEV.
 c          NP may be less than NCV-NEV since a leading block of the current
 c          upper Tridiagonal matrix has split off and contains "unwanted"
 c          Ritz values.

PARPACK/SRC/BLACS/pdsaupd.f

@@ -271,9 +271,9 @@ c          = -8: Error return from trid. eigenvalue calculation;
 c                Informatinal error from LAPACK routine dsteqr .
 c          = -9: Starting vector is zero.
 c          = -10: IPARAM(7) must be 1,2,3,4,5.
-c          = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatable.
+c          = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatible.
 c          = -12: IPARAM(1) must be equal to 0 or 1.
-c          = -13: NEV and WHICH = 'BE' are incompatable.
+c          = -13: NEV and WHICH = 'BE' are incompatible.
 c          = -9999: Could not build an Arnoldi factorization.
 c                   IPARAM(5) returns the size of the current Arnoldi
 c                   factorization. The user is advised to check that
@@ -310,7 +310,7 @@ c     also increases the work and storage required to maintain the orthogonal
 c     basis vectors.   The optimal "cross-over" with respect to CPU time
 c     is problem dependent and must be determined empirically.
 c
-c  5. If IPARAM(7) = 2 then in the Reverse commuication interface the user
+c  5. If IPARAM(7) = 2 then in the Reverse communication interface the user
 c     must do the following. When IDO = 1, Y = OP * X is to be computed.
 c     When IPARAM(7) = 2 OP = inv(B)*A. After computing A*X the user
 c     must overwrite X with A*X. Y is then the solution to the linear set

PARPACK/SRC/BLACS/pslamch.f

@@ -1,6 +1,6 @@
       REAL               FUNCTION PSLAMCH( ICTXT, CMACH )
 *
-*  -- ScaLAPACK auxilliary routine (version 1.0) --
+*  -- ScaLAPACK auxiliary routine (version 1.0) --
 *     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
 *     and University of California, Berkeley.
 *     February 28, 1995

PARPACK/SRC/BLACS/psnaitr.f

@@ -398,7 +398,7 @@ c
 c 
 c        %---------------------------------------------------%
 c        | STEP 1: Check if the B norm of j-th residual      |
-c        | vector is zero. Equivalent to determing whether   |
+c        | vector is zero. Equivalent to determining whether |
 c        | an exact j-step Arnoldi factorization is present. |
 c        %---------------------------------------------------%
 c

PARPACK/SRC/BLACS/psnapps.f

@@ -34,7 +34,7 @@ c
 c  KEV     Integer.  (INPUT/OUTPUT)
 c          KEV+NP is the size of the input matrix H.
 c          KEV is the size of the updated matrix HNEW.  KEV is only 
-c          updated on ouput when fewer than NP shifts are applied in
+c          updated on output when fewer than NP shifts are applied in
 c          order to keep the conjugate pair together.
 c
 c  NP      Integer.  (INPUT)
@@ -215,7 +215,7 @@ c
       intrinsic  abs, max, min
 c
 c     %----------------%
-c     | Data statments |
+c     | Data statements |
 c     %----------------%
 c
       data       first / .true. /
@@ -395,7 +395,7 @@ c
             do 80 i = istart, iend-1
 c
 c              %-----------------------------------------------------%
-c              | Contruct the plane rotation G to zero out the bulge |
+c              | Construct the plane rotation G to zero out the bulge |
 c              %-----------------------------------------------------%
 c
                call slartg (f, g, c, s, r)

PARPACK/SRC/BLACS/psnaup2.f

@@ -28,7 +28,7 @@ c          products (involving the operator OP) per Arnoldi iteration.
 c          The logic for adjusting is contained within the current
 c          subroutine.
 c          If ISHIFT=0, NP is the number of shifts the user needs
-c          to provide via reverse comunication. 0 < NP < NCV-NEV.
+c          to provide via reverse communication. 0 < NP < NCV-NEV.
 c          NP may be less than NCV-NEV for two reasons. The first, is
 c          to keep complex conjugate pairs of "wanted" Ritz values 
 c          together. The second, is that a leading block of the current
@@ -715,7 +715,7 @@ c
          if (ishift .eq. 0) then
 c
 c           %-------------------------------------------------------%
-c           | User specified shifts: reverse comminucation to       |
+c           | User specified shifts: reverse communication to       |
 c           | compute the shifts. They are returned in the first    |
 c           | 2*NP locations of WORKL.                              |
 c           %-------------------------------------------------------%

PARPACK/SRC/BLACS/psnaupd.f

@@ -276,7 +276,7 @@ c          = -7: Length of private work array is not sufficient.
 c          = -8: Error return from LAPACK eigenvalue calculation;
 c          = -9: Starting vector is zero.
 c          = -10: IPARAM(7) must be 1,2,3,4.
-c          = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatable.
+c          = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatible.
 c          = -12: IPARAM(1) must be equal to 0 or 1.
 c          = -9999: Could not build an Arnoldi factorization.
 c                   IPARAM(5) returns the size of the current Arnoldi

PARPACK/SRC/BLACS/pssaitr.f

@@ -396,8 +396,8 @@ c
          end if
 c
 c        %---------------------------------------------------------%
-c        | Check for exact zero. Equivalent to determing whether a |
-c        | j-step Arnoldi factorization is present.                |
+c        | Check for exact zero. Equivalent to determining whether |
+c        | a j-step Arnoldi factorization is present.              |
 c        %---------------------------------------------------------%
 c
          if (rnorm .gt. zero) go to 40

PARPACK/SRC/BLACS/pssapps.f

@@ -204,7 +204,7 @@ c
       intrinsic  abs
 c
 c     %----------------%
-c     | Data statments |
+c     | Data statements |
 c     %----------------%
 c
       data       first / .true. /
@@ -275,7 +275,7 @@ c
                   call pivout (comm, logfil, 1, i, ndigit, 
      &                 '_sapps: deflation at row/column no.')
                   call pivout (comm, logfil, 1, jj, ndigit, 
-     &                 '_sapps: occured before shift number.')
+     &                 '_sapps: occurred before shift number.')
                   call psvout (comm, logfil, 1, h(i+1,1), ndigit, 
      &                 '_sapps: the corresponding off diagonal element')
                end if

PARPACK/SRC/BLACS/pssaup2.f

@@ -29,7 +29,7 @@ c          products (involving the operator OP) per Arnoldi iteration.
 c          The logic for adjusting is contained within the current
 c          subroutine.
 c          If ISHIFT=0, NP is the number of shifts the user needs
-c          to provide via reverse comunication. 0 < NP < NCV-NEV.
+c          to provide via reverse communication. 0 < NP < NCV-NEV.
 c          NP may be less than NCV-NEV since a leading block of the current
 c          upper Tridiagonal matrix has split off and contains "unwanted"
 c          Ritz values.

PARPACK/SRC/BLACS/pssaupd.f

@@ -271,9 +271,9 @@ c          = -8: Error return from trid. eigenvalue calculation;
 c                Informatinal error from LAPACK routine ssteqr.
 c          = -9: Starting vector is zero.
 c          = -10: IPARAM(7) must be 1,2,3,4,5.
-c          = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatable.
+c          = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatible.
 c          = -12: IPARAM(1) must be equal to 0 or 1.
-c          = -13: NEV and WHICH = 'BE' are incompatable.
+c          = -13: NEV and WHICH = 'BE' are incompatible.
 c          = -9999: Could not build an Arnoldi factorization.
 c                   IPARAM(5) returns the size of the current Arnoldi
 c                   factorization. The user is advised to check that
@@ -310,7 +310,7 @@ c     also increases the work and storage required to maintain the orthogonal
 c     basis vectors.   The optimal "cross-over" with respect to CPU time
 c     is problem dependent and must be determined empirically.
 c
-c  5. If IPARAM(7) = 2 then in the Reverse commuication interface the user
+c  5. If IPARAM(7) = 2 then in the Reverse communication interface the user
 c     must do the following. When IDO = 1, Y = OP * X is to be computed.
 c     When IPARAM(7) = 2 OP = inv(B)*A. After computing A*X the user
 c     must overwrite X with A*X. Y is then the solution to the linear set

PARPACK/SRC/BLACS/pznapps.f

@@ -225,7 +225,7 @@ c
       cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( cdum ) )
 c
 c     %----------------%
-c     | Data statments |
+c     | Data statements |
 c     %----------------%
 c
       data       first / .true. /

PARPACK/SRC/BLACS/pznaup2.f

@@ -28,7 +28,7 @@ c          products (involving the operator OP) per Arnoldi iteration.
 c          The logic for adjusting is contained within the current
 c          subroutine.
 c          If ISHIFT=0, NP is the number of shifts the user needs
-c          to provide via reverse comunication. 0 < NP < NCV-NEV.
+c          to provide via reverse communication. 0 < NP < NCV-NEV.
 c          NP may be less than NCV-NEV since a leading block of the current
 c          upper Hessenberg matrix has split off and contains "unwanted"
 c          Ritz values.

PARPACK/SRC/BLACS/pznaupd.f

@@ -258,7 +258,7 @@ c          = -7: Length of private work array is not sufficient.
 c          = -8: Error return from LAPACK eigenvalue calculation;
 c          = -9: Starting vector is zero.
 c          = -10: IPARAM(7) must be 1,2,3.
-c          = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatable.
+c          = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatible.
 c          = -12: IPARAM(1) must be equal to 0 or 1.
 c          = -9999: Could not build an Arnoldi factorization.
 c                   User input error highly likely.  Please

PARPACK/SRC/MPI/Makefile.am

@@ -1,7 +1,10 @@
+AUTOMAKE_OPTIONS = subdir-objects # Needed as debug_init/icb.f90 are not in current directory.
+
 F77 = $(MPIF77)
 FFLAGS_SAV = @FFLAGS@
 FFLAGS =
 
+PSRC = pcontext.f
 SSRC = psnaitr.f psnapps.f psnaup2.f psnaupd.f psneigh.f psneupd.f psngets.f \
        pssaitr.f pssapps.f pssaup2.f pssaupd.f psseigt.f psseupd.f pssgets.f \
        psgetv0.f pslarnv.f psnorm2.f
@@ -16,13 +19,29 @@ CSRC = pcnaitr.f pcnapps.f pcnaup2.f pcnaupd.f pcneigh.f pcneupd.f pcngets.f \
 ZSRC = pznaitr.f pznapps.f pznaup2.f pznaupd.f pzneigh.f pzneupd.f pzngets.f \
        pzgetv0.f pzlarnv.f pdznorm2.f
 
-EXTRA_DIST = debug.h stat.h
+if ICB
+SSRC += icbpss.f90 icbpsn.f90
+DSRC += icbpds.f90 icbpdn.f90
+CSRC += icbpcn.f90
+ZSRC += icbpzn.f90
+endif
 
-noinst_LTLIBRARIES = libparpack_noopt.la libparpacksrcmpi.la
+EXTRA_DIST = debug.h stat.h pcontext.h
 
-libparpack_noopt_la_SOURCES = pslamch.f pdlamch10.f
-libparpack_noopt_la_FFLAGS = -O0
+noinst_LTLIBRARIES = libparpack@LIBSUFFIX@_noopt.la
+libparpack@LIBSUFFIX@_noopt_la_SOURCES = pslamch10.f pdlamch10.f
+libparpack@LIBSUFFIX@_noopt_la_FFLAGS = -O0
 
-libparpacksrcmpi_la_SOURCES = $(SSRC) $(DSRC) $(CSRC) $(ZSRC)
-libparpacksrcmpi_la_FFLAGS = $(FFLAGS_SAV)
-libparpacksrcmpi_la_LIBADD = libparpack_noopt.la
+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_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

PARPACK/SRC/MPI/icbpcn.f90

@@ -0,0 +1,64 @@
+! icbp : iso_c_binding for parpack
+
+subroutine pcnaupd_c(comm, ido, bmat, n, which, nev, tol, resid, ncv, v, ldv,&
+                     iparam, ipntr, workd, workl, lworkl, rwork, info)       &
+                     bind(c, name="pcnaupd_c")
+  use :: iso_c_binding
+  implicit none
+  integer(kind=c_int),         value,               intent(in)    :: comm
+  integer(kind=c_int),                              intent(inout) :: ido
+  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
+  integer(kind=c_int),         value,               intent(in)    :: nev
+  real(kind=c_float_complex),  value,               intent(in)    :: tol
+  real(kind=c_float_complex),  dimension(n),        intent(inout) :: resid
+  integer(kind=c_int),         value,               intent(in)    :: ncv
+  real(kind=c_float_complex),  dimension(ldv, ncv), intent(out)   :: v
+  integer(kind=c_int),         value,               intent(in)    :: ldv
+  integer(kind=c_int),         dimension(11),       intent(inout) :: iparam
+  integer(kind=c_int),         dimension(11),       intent(out)   :: ipntr
+  real(kind=c_float_complex),  dimension(3*n),      intent(out)   :: workd
+  real(kind=c_float_complex),  dimension(lworkl),   intent(out)   :: workl
+  integer(kind=c_int),         value,               intent(in)    :: lworkl
+  real(kind=c_float_complex),  dimension(ncv),      intent(out)   :: rwork
+  integer(kind=c_int),                              intent(inout) :: info
+  call pcnaupd(comm, ido, bmat, n, which, nev, tol, resid, ncv, v, ldv,&
+               iparam, ipntr, workd, workl, lworkl, rwork, info)
+end subroutine pcnaupd_c
+
+subroutine pcneupd_c(comm, rvec, howmny, select, d, z, ldz, sigma, workev,&
+                     bmat, n, which, nev, tol, resid, ncv, v, ldv,        &
+                     iparam, ipntr, workd, workl, lworkl, rwork, info)    &
+                     bind(c, name="pcneupd_c")
+  use :: iso_c_binding
+  implicit none
+  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
+  logical(kind=c_bool),        dimension(ncv),      intent(in)    :: select
+  real(kind=c_float_complex),  dimension(nev),      intent(out)   :: d
+  real(kind=c_float_complex),  dimension(n, nev),   intent(out)   :: z
+  integer(kind=c_int),         value,               intent(in)    :: ldz
+  real(kind=c_float_complex),  value,               intent(in)    :: sigma
+  real(kind=c_float_complex),  dimension(2*ncv),    intent(in)    :: 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
+  integer(kind=c_int),         value,               intent(in)    :: nev
+  real(kind=c_float_complex),  value,               intent(in)    :: tol
+  real(kind=c_float_complex),  dimension(n),        intent(inout) :: resid
+  integer(kind=c_int),         value,               intent(in)    :: ncv
+  real(kind=c_float_complex),  dimension(ldv, ncv), intent(out)   :: v
+  integer(kind=c_int),         value,               intent(in)    :: ldv
+  integer(kind=c_int),         dimension(11),       intent(inout) :: iparam
+  integer(kind=c_int),         dimension(11),       intent(out)   :: ipntr
+  real(kind=c_float_complex),  dimension(3*n),      intent(out)   :: workd
+  real(kind=c_float_complex),  dimension(lworkl),   intent(out)   :: workl
+  integer(kind=c_int),         value,               intent(in)    :: lworkl
+  real(kind=c_float_complex),  dimension(ncv),      intent(out)   :: rwork
+  integer(kind=c_int),                              intent(inout) :: info
+  call pcneupd(comm, rvec, howmny, select, d, z, ldz, sigma, workev,&
+               bmat, n, which, nev, tol, resid, ncv, v, ldv,        &
+               iparam, ipntr, workd, workl, lworkl, rwork, info)
+end subroutine pcneupd_c

PARPACK/SRC/MPI/icbpdn.f90

@@ -0,0 +1,65 @@
+! icbp : iso_c_binding for parpack
+
+subroutine pdnaupd_c(comm, ido, bmat, n, which, nev, tol, resid, ncv, v, ldv,&
+                     iparam, ipntr, workd, workl, lworkl, info)              &
+                     bind(c, name="pdnaupd_c")
+  use :: iso_c_binding
+  implicit none
+  integer(kind=c_int),    value,               intent(in)    :: comm
+  integer(kind=c_int),                         intent(inout) :: ido
+  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
+  integer(kind=c_int),    value,               intent(in)    :: nev
+  real(kind=c_double),    value,               intent(in)    :: tol
+  real(kind=c_double),    dimension(n),        intent(inout) :: resid
+  integer(kind=c_int),    value,               intent(in)    :: ncv
+  real(kind=c_double),    dimension(ldv, ncv), intent(out)   :: v
+  integer(kind=c_int),    value,               intent(in)    :: ldv
+  integer(kind=c_int),    dimension(11),       intent(inout) :: iparam
+  integer(kind=c_int),    dimension(11),       intent(out)   :: ipntr
+  real(kind=c_double),    dimension(3*n),      intent(out)   :: workd
+  real(kind=c_double),    dimension(lworkl),   intent(out)   :: workl
+  integer(kind=c_int),    value,               intent(in)    :: lworkl
+  integer(kind=c_int),                         intent(inout) :: info
+  call pdnaupd(comm, ido, bmat, n, which, nev, tol, resid, ncv, v, ldv,&
+               iparam, ipntr, workd, workl, lworkl, info)
+end subroutine pdnaupd_c
+
+subroutine pdneupd_c(comm, rvec, howmny, select,                  &
+                     dr, di, z, ldz, sigmar, sigmai,              &
+                     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
+  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
+  logical(kind=c_bool),   dimension(ncv),      intent(in)    :: select
+  real(kind=c_double),    dimension(nev+1),    intent(out)   :: dr
+  real(kind=c_double),    dimension(nev+1),    intent(out)   :: di
+  real(kind=c_double),    dimension(n, nev+1), intent(out)   :: z
+  integer(kind=c_int),    value,               intent(in)    :: ldz
+  real(kind=c_double),    value,               intent(in)    :: sigmar
+  real(kind=c_double),    value,               intent(in)    :: sigmai
+  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
+  integer(kind=c_int),    value,               intent(in)    :: nev
+  real(kind=c_double),    value,               intent(in)    :: tol
+  real(kind=c_double),    dimension(n),        intent(inout) :: resid
+  integer(kind=c_int),    value,               intent(in)    :: ncv
+  real(kind=c_double),    dimension(ldv, ncv), intent(out)   :: v
+  integer(kind=c_int),    value,               intent(in)    :: ldv
+  integer(kind=c_int),    dimension(11),       intent(inout) :: iparam
+  integer(kind=c_int),    dimension(11),       intent(out)   :: ipntr
+  real(kind=c_double),    dimension(3*n),      intent(out)   :: workd
+  real(kind=c_double),    dimension(lworkl),   intent(out)   :: workl
+  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,              &
+               bmat, n, which, nev, tol, resid, ncv, v, ldv,&
+               iparam, ipntr, workd, workl, lworkl, info)
+end subroutine pdneupd_c

PARPACK/SRC/MPI/icbpds.f90

@@ -0,0 +1,61 @@
+! icbp : iso_c_binding for parpack
+
+subroutine pdsaupd_c(comm, ido, bmat, n, which, nev, tol, resid, ncv, v, ldv,&
+                     iparam, ipntr, workd, workl, lworkl, info)              &
+                     bind(c, name="pdsaupd_c")
+  use :: iso_c_binding
+  implicit none
+  integer(kind=c_int),    value,               intent(in)    :: comm
+  integer(kind=c_int),                         intent(inout) :: ido
+  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
+  integer(kind=c_int),    value,               intent(in)    :: nev
+  real(kind=c_double),    value,               intent(in)    :: tol
+  real(kind=c_double),    dimension(n),        intent(inout) :: resid
+  integer(kind=c_int),    value,               intent(in)    :: ncv
+  real(kind=c_double),    dimension(ldv, ncv), intent(out)   :: v
+  integer(kind=c_int),    value,               intent(in)    :: ldv
+  integer(kind=c_int),    dimension(11),       intent(inout) :: iparam
+  integer(kind=c_int),    dimension(11),       intent(out)   :: ipntr
+  real(kind=c_double),    dimension(3*n),      intent(out)   :: workd
+  real(kind=c_double),    dimension(lworkl),   intent(out)   :: workl
+  integer(kind=c_int),    value,               intent(in)    :: lworkl
+  integer(kind=c_int),                         intent(inout) :: info
+  call pdsaupd(comm, ido, bmat, n, which, nev, tol, resid, ncv, v, ldv,&
+               iparam, ipntr, workd, workl, lworkl, info)
+end subroutine pdsaupd_c
+
+subroutine pdseupd_c(comm, rvec, howmny, select, d, z, ldz, sigma,&
+                     bmat, n, which, nev, tol, resid, ncv, v, ldv,&
+                     iparam, ipntr, workd, workl, lworkl, info)   &
+                     bind(c, name="pdseupd_c")
+  use :: iso_c_binding
+  implicit none
+  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
+  logical(kind=c_bool),   dimension(ncv),      intent(in)    :: select
+  real(kind=c_double),    dimension(nev),      intent(out)   :: d
+  real(kind=c_double),    dimension(n, nev),   intent(out)   :: z
+  integer(kind=c_int),    value,               intent(in)    :: ldz
+  real(kind=c_double),    value,               intent(in)    :: sigma
+  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
+  integer(kind=c_int),    value,               intent(in)    :: nev
+  real(kind=c_double),    value,               intent(in)    :: tol
+  real(kind=c_double),    dimension(n),        intent(inout) :: resid
+  integer(kind=c_int),    value,               intent(in)    :: ncv
+  real(kind=c_double),    dimension(ldv, ncv), intent(out)   :: v
+  integer(kind=c_int),    value,               intent(in)    :: ldv
+  integer(kind=c_int),    dimension(11),       intent(inout) :: iparam
+  integer(kind=c_int),    dimension(11),       intent(out)   :: ipntr
+  real(kind=c_double),    dimension(3*n),      intent(out)   :: workd
+  real(kind=c_double),    dimension(lworkl),   intent(out)   :: workl
+  integer(kind=c_int),    value,               intent(in)    :: lworkl
+  integer(kind=c_int),                         intent(inout) :: info
+  call pdseupd(comm, rvec, howmny, select, d, z, ldz, sigma,&
+               bmat, n, which, nev, tol, resid, ncv, v, ldv,&
+               iparam, ipntr, workd, workl, lworkl, info)
+end subroutine pdseupd_c

PARPACK/SRC/MPI/icbpsn.f90

@@ -0,0 +1,65 @@
+! icbp : iso_c_binding for parpack
+
+subroutine psnaupd_c(comm, ido, bmat, n, which, nev, tol, resid, ncv, v, ldv,&
+                     iparam, ipntr, workd, workl, lworkl, info)              &
+                     bind(c, name="psnaupd_c")
+  use :: iso_c_binding
+  implicit none
+  integer(kind=c_int),    value,               intent(in)    :: comm
+  integer(kind=c_int),                         intent(inout) :: ido
+  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
+  integer(kind=c_int),    value,               intent(in)    :: nev
+  real(kind=c_float),     value,               intent(in)    :: tol
+  real(kind=c_float),     dimension(n),        intent(inout) :: resid
+  integer(kind=c_int),    value,               intent(in)    :: ncv
+  real(kind=c_float),     dimension(ldv, ncv), intent(out)   :: v
+  integer(kind=c_int),    value,               intent(in)    :: ldv
+  integer(kind=c_int),    dimension(11),       intent(inout) :: iparam
+  integer(kind=c_int),    dimension(11),       intent(out)   :: ipntr
+  real(kind=c_float),     dimension(3*n),      intent(out)   :: workd
+  real(kind=c_float),     dimension(lworkl),   intent(out)   :: workl
+  integer(kind=c_int),    value,               intent(in)    :: lworkl
+  integer(kind=c_int),                         intent(inout) :: info
+  call psnaupd(comm, ido, bmat, n, which, nev, tol, resid, ncv, v, ldv,&
+               iparam, ipntr, workd, workl, lworkl, info)
+end subroutine psnaupd_c
+
+subroutine psneupd_c(comm, rvec, howmny, select,                  &
+                     dr, di, z, ldz, sigmar, sigmai,              &
+                     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
+  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
+  logical(kind=c_bool),   dimension(ncv),      intent(in)    :: select
+  real(kind=c_float),     dimension(nev+1),    intent(out)   :: dr
+  real(kind=c_float),     dimension(nev+1),    intent(out)   :: di
+  real(kind=c_float),     dimension(n, nev+1), intent(out)   :: z
+  integer(kind=c_int),    value,               intent(in)    :: ldz
+  real(kind=c_float),     value,               intent(in)    :: sigmar
+  real(kind=c_float),     value,               intent(in)    :: sigmai
+  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
+  integer(kind=c_int),    value,               intent(in)    :: nev
+  real(kind=c_float),     value,               intent(in)    :: tol
+  real(kind=c_float),     dimension(n),        intent(inout) :: resid
+  integer(kind=c_int),    value,               intent(in)    :: ncv
+  real(kind=c_float),     dimension(ldv, ncv), intent(out)   :: v
+  integer(kind=c_int),    value,               intent(in)    :: ldv
+  integer(kind=c_int),    dimension(11),       intent(inout) :: iparam
+  integer(kind=c_int),    dimension(11),       intent(out)   :: ipntr
+  real(kind=c_float),     dimension(3*n),      intent(out)   :: workd
+  real(kind=c_float),     dimension(lworkl),   intent(out)   :: workl
+  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,              &
+               bmat, n, which, nev, tol, resid, ncv, v, ldv,&
+               iparam, ipntr, workd, workl, lworkl, info)
+end subroutine psneupd_c

PARPACK/SRC/MPI/icbpss.f90

@@ -0,0 +1,61 @@
+! icbp : iso_c_binding for parpack
+
+subroutine pssaupd_c(comm, ido, bmat, n, which, nev, tol, resid, ncv, v, ldv,&
+                     iparam, ipntr, workd, workl, lworkl, info)              &
+                     bind(c, name="pssaupd_c")
+  use :: iso_c_binding
+  implicit none
+  integer(kind=c_int),    value,               intent(in)    :: comm
+  integer(kind=c_int),                         intent(inout) :: ido
+  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
+  integer(kind=c_int),    value,               intent(in)    :: nev
+  real(kind=c_float),     value,               intent(in)    :: tol
+  real(kind=c_float),     dimension(n),        intent(inout) :: resid
+  integer(kind=c_int),    value,               intent(in)    :: ncv
+  real(kind=c_float),     dimension(ldv, ncv), intent(out)   :: v
+  integer(kind=c_int),    value,               intent(in)    :: ldv
+  integer(kind=c_int),    dimension(11),       intent(inout) :: iparam
+  integer(kind=c_int),    dimension(11),       intent(out)   :: ipntr
+  real(kind=c_float),     dimension(3*n),      intent(out)   :: workd
+  real(kind=c_float),     dimension(lworkl),   intent(out)   :: workl
+  integer(kind=c_int),    value,               intent(in)    :: lworkl
+  integer(kind=c_int),                         intent(inout) :: info
+  call pssaupd(comm, ido, bmat, n, which, nev, tol, resid, ncv, v, ldv,&
+               iparam, ipntr, workd, workl, lworkl, info)
+end subroutine pssaupd_c
+
+subroutine psseupd_c(comm, rvec, howmny, select, d, z, ldz, sigma,&
+                     bmat, n, which, nev, tol, resid, ncv, v, ldv,&
+                     iparam, ipntr, workd, workl, lworkl, info)   &
+                     bind(c, name="psseupd_c")
+  use :: iso_c_binding
+  implicit none
+  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
+  logical(kind=c_bool),   dimension(ncv),      intent(in)    :: select
+  real(kind=c_float),     dimension(nev),      intent(out)   :: d
+  real(kind=c_float),     dimension(n, nev),   intent(out)   :: z
+  integer(kind=c_int),    value,               intent(in)    :: ldz
+  real(kind=c_float),     value,               intent(in)    :: sigma
+  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
+  integer(kind=c_int),    value,               intent(in)    :: nev
+  real(kind=c_float),     value,               intent(in)    :: tol
+  real(kind=c_float),     dimension(n),        intent(inout) :: resid
+  integer(kind=c_int),    value,               intent(in)    :: ncv
+  real(kind=c_float),     dimension(ldv, ncv), intent(out)   :: v
+  integer(kind=c_int),    value,               intent(in)    :: ldv
+  integer(kind=c_int),    dimension(11),       intent(inout) :: iparam
+  integer(kind=c_int),    dimension(11),       intent(out)   :: ipntr
+  real(kind=c_float),     dimension(3*n),      intent(out)   :: workd
+  real(kind=c_float),     dimension(lworkl),   intent(out)   :: workl
+  integer(kind=c_int),    value,               intent(in)    :: lworkl
+  integer(kind=c_int),                         intent(inout) :: info
+  call psseupd(comm, rvec, howmny, select, d, z, ldz, sigma,&
+               bmat, n, which, nev, tol, resid, ncv, v, ldv,&
+               iparam, ipntr, workd, workl, lworkl, info)
+end subroutine psseupd_c

PARPACK/SRC/MPI/icbpzn.f90

@@ -0,0 +1,64 @@
+! icbp : iso_c_binding for parpack
+
+subroutine pznaupd_c(comm, ido, bmat, n, which, nev, tol, resid, ncv, v, ldv,&
+                     iparam, ipntr, workd, workl, lworkl, rwork, info)       &
+                     bind(c, name="pznaupd_c")
+  use :: iso_c_binding
+  implicit none
+  integer(kind=c_int),         value,               intent(in)    :: comm
+  integer(kind=c_int),                              intent(inout) :: ido
+  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
+  integer(kind=c_int),         value,               intent(in)    :: nev
+  real(kind=c_double_complex), value,               intent(in)    :: tol
+  real(kind=c_double_complex), dimension(n),        intent(inout) :: resid
+  integer(kind=c_int),         value,               intent(in)    :: ncv
+  real(kind=c_double_complex), dimension(ldv, ncv), intent(out)   :: v
+  integer(kind=c_int),         value,               intent(in)    :: ldv
+  integer(kind=c_int),         dimension(11),       intent(inout) :: iparam
+  integer(kind=c_int),         dimension(11),       intent(out)   :: ipntr
+  real(kind=c_double_complex), dimension(3*n),      intent(out)   :: workd
+  real(kind=c_double_complex), dimension(lworkl),   intent(out)   :: workl
+  integer(kind=c_int),         value,               intent(in)    :: lworkl
+  real(kind=c_double_complex), dimension(ncv),      intent(out)   :: rwork
+  integer(kind=c_int),                              intent(inout) :: info
+  call pznaupd(comm, ido, bmat, n, which, nev, tol, resid, ncv, v, ldv,&
+               iparam, ipntr, workd, workl, lworkl, rwork, info)
+end subroutine pznaupd_c
+
+subroutine pzneupd_c(comm, rvec, howmny, select, d, z, ldz, sigma, workev,&
+                     bmat, n, which, nev, tol, resid, ncv, v, ldv,        &
+                     iparam, ipntr, workd, workl, lworkl, rwork, info)    &
+                     bind(c, name="pzneupd_c")
+  use :: iso_c_binding
+  implicit none
+  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
+  logical(kind=c_bool),        dimension(ncv),      intent(in)    :: select
+  real(kind=c_double_complex), dimension(nev),      intent(out)   :: d
+  real(kind=c_double_complex), dimension(n, nev),   intent(out)   :: z
+  integer(kind=c_int),         value,               intent(in)    :: ldz
+  real(kind=c_double_complex), value,               intent(in)    :: sigma
+  real(kind=c_double_complex), dimension(2*ncv),    intent(in)    :: 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
+  integer(kind=c_int),         value,               intent(in)    :: nev
+  real(kind=c_double_complex), value,               intent(in)    :: tol
+  real(kind=c_double_complex), dimension(n),        intent(inout) :: resid
+  integer(kind=c_int),         value,               intent(in)    :: ncv
+  real(kind=c_double_complex), dimension(ldv, ncv), intent(out)   :: v
+  integer(kind=c_int),         value,               intent(in)    :: ldv
+  integer(kind=c_int),         dimension(11),       intent(inout) :: iparam
+  integer(kind=c_int),         dimension(11),       intent(out)   :: ipntr
+  real(kind=c_double_complex), dimension(3*n),      intent(out)   :: workd
+  real(kind=c_double_complex), dimension(lworkl),   intent(out)   :: workl
+  integer(kind=c_int),         value,               intent(in)    :: lworkl
+  real(kind=c_double_complex), dimension(ncv),      intent(out)   :: rwork
+  integer(kind=c_int),                              intent(inout) :: info
+  call pzneupd(comm, rvec, howmny, select, d, z, ldz, sigma, workev,&
+               bmat, n, which, nev, tol, resid, ncv, v, ldv,        &
+               iparam, ipntr, workd, workl, lworkl, rwork, info)
+end subroutine pzneupd_c

PARPACK/SRC/MPI/pcnaitr.f

@@ -144,7 +144,7 @@ c     clanhs    LAPACK routine that computes various norms of a matrix.
 c     clascl    LAPACK routine for careful scaling of a matrix.
 c     slabad    LAPACK routine for defining the underflow and overflow
 c               limits
-c     pslamch   ScaLAPACK routine that determines machine constants.
+c     pslamch10   ScaLAPACK routine that determines machine constants.
 c     slapy2    LAPACK routine to compute sqrt(x**2+y**2) carefully.
 c     cgemv     Level 2 BLAS routine for matrix vector multiplication.
 c     caxpy     Level 1 BLAS that computes a vector triad.
@@ -223,6 +223,7 @@ c
      &   (comm, ido, bmat, n, k, np, nb, resid, rnorm, v, ldv, h, ldh, 
      &    ipntr, workd, workl, info)
 c
+      include   'pcontext.h'
       include   'mpif.h'
 c
 c     %---------------%
@@ -278,7 +279,7 @@ c     %---------------%
 c     | Local Scalars |
 c     %---------------%
 c
-      logical    first, orth1, orth2, rstart, step3, step4
+      logical    orth1, orth2, rstart, step3, step4
       integer    ierr, i, infol, ipj, irj, ivj, iter, itry, j, msglvl,
      &           jj
       Real
@@ -287,7 +288,7 @@ c
       Complex
      &           cnorm
 c
-      save       first, orth1, orth2, rstart, step3, step4,
+      save       orth1, orth2, rstart, step3, step4,
      &           ierr, ipj, irj, ivj, iter, itry, j, msglvl, ovfl,
      &           betaj, rnorm1, smlnum, ulp, unfl, wnorm
 c
@@ -308,8 +309,8 @@ c
       Complex
      &           cdotc 
       Real            
-     &           pslamch, pscnorm2, clanhs, slapy2
-      external   cdotc, pscnorm2, clanhs, pslamch, slapy2
+     &           pslamch10, pscnorm2, clanhs, slapy2
+      external   cdotc, pscnorm2, clanhs, pslamch10, slapy2
 c
 c     %---------------------%
 c     | Intrinsic Functions |
@@ -321,13 +322,12 @@ c     %-----------------%
 c     | Data statements |
 c     %-----------------%
 c
-      data       first / .true. /
 c
 c     %-----------------------%
 c     | Executable Statements |
 c     %-----------------------%
 c
-      if (first) then
+      if (aitr_first) then
 c
 c        %-----------------------------------------%
 c        | Set machine-dependent constants for the |
@@ -337,12 +337,12 @@ c        | overflow should not occur.              |
 c        | REFERENCE: LAPACK subroutine clahqr     |
 c        %-----------------------------------------%
 c
-         unfl = pslamch(comm,  'safe minimum' )
+         unfl = pslamch10(comm,  'safe minimum' )
          ovfl = real(one / unfl)
          call slabad( unfl, ovfl )
-         ulp = pslamch( comm, 'precision' )
+         ulp = pslamch10( comm, 'precision' )
          smlnum = unfl*( n / ulp )
-         first = .false.
+         aitr_first = .false.
       end if
 c
       if (ido .eq. 0) then

PARPACK/SRC/MPI/pcnapps.f

@@ -105,7 +105,7 @@ c     clartg  LAPACK Givens rotation construction routine.
 c     claset  LAPACK matrix initialization routine.
 c     slabad  LAPACK routine for defining the underflow and overflow
 c             limits.
-c     pslamch ScaLAPACK routine that determines machine constants.
+c     pslamch10 ScaLAPACK routine that determines machine constants.
 c     slapy2  LAPACK routine to compute sqrt(x**2+y**2) carefully.
 c     cgemv   Level 2 BLAS routine for matrix vector multiplication.
 c     caxpy   Level 1 BLAS that computes a vector triad.
@@ -149,6 +149,7 @@ c     %--------------------%
 c     | MPI Communicator |
 c     %--------------------%
 c
+      include   'pcontext.h'
       integer   comm
 c
 c     %----------------------------------------------------%
@@ -188,12 +189,11 @@ c     | Local Scalars & Arrays |
 c     %------------------------%
 c
       integer    i, iend, istart, j, jj, kplusp, msglvl
-      logical    first
       Complex
      &           cdum, f, g, h11, h21, r, s, sigma, t
       Real             
      &           c,  ovfl, smlnum, ulp, unfl, tst1
-      save       first, ovfl, smlnum, ulp, unfl 
+      save       ovfl, smlnum, ulp, unfl 
 c
 c     %----------------------%
 c     | External Subroutines |
@@ -207,8 +207,8 @@ c     | External Functions |
 c     %--------------------%
 c
       Real                 
-     &           clanhs, pslamch, slapy2
-      external   clanhs, pslamch, slapy2
+     &           clanhs, pslamch10, slapy2
+      external   clanhs, pslamch10, slapy2
 c
 c     %----------------------%
 c     | Intrinsics Functions |
@@ -225,16 +225,15 @@ c
       cabs1( cdum ) = abs( real( cdum ) ) + abs( aimag( cdum ) )
 c
 c     %----------------%
-c     | Data statments |
+c     | Data statements |
 c     %----------------%
 c
-      data       first / .true. /
 c
 c     %-----------------------%
 c     | Executable Statements |
 c     %-----------------------%
 c
-      if (first) then
+      if (apps_first) then
 c
 c        %-----------------------------------------------%
 c        | Set machine-dependent constants for the       |
@@ -243,12 +242,12 @@ c        | overflow should not occur.                    |
 c        | REFERENCE: LAPACK subroutine clahqr           |
 c        %-----------------------------------------------%
 c
-         unfl = pslamch( comm, 'safe minimum' )
+         unfl = pslamch10( comm, 'safe minimum' )
          ovfl = real(one / unfl)
          call slabad( unfl, ovfl )
-         ulp = pslamch( comm,  'precision' )
+         ulp = pslamch10( comm,  'precision' )
          smlnum = unfl*( n / ulp )
-         first = .false.
+         apps_first = .false.
       end if
 c
 c     %-------------------------------%

PARPACK/SRC/MPI/pcnaup2.f

@@ -28,7 +28,7 @@ c          products (involving the operator OP) per Arnoldi iteration.
 c          The logic for adjusting is contained within the current
 c          subroutine.
 c          If ISHIFT=0, NP is the number of shifts the user needs
-c          to provide via reverse comunication. 0 < NP < NCV-NEV.
+c          to provide via reverse communication. 0 < NP < NCV-NEV.
 c          NP may be less than NCV-NEV since a leading block of the current
 c          upper Hessenberg matrix has split off and contains "unwanted"
 c          Ritz values.
@@ -141,7 +141,7 @@ c     arscnd   ARPACK utility routine for timing.
 c     pcmout   Parallel ARPACK utility routine that prints matrices
 c     pcvout   Parallel ARPACK utility routine that prints vectors.
 c     psvout   ARPACK utility routine that prints vectors.
-c     pslamch  ScaLAPACK routine that determines machine constants.
+c     pslamch10  ScaLAPACK routine that determines machine constants.
 c     slapy2   LAPACK routine to compute sqrt(x**2+y**2) carefully.
 c     ccopy    Level 1 BLAS that copies one vector to another .
 c     cdotc    Level 1 BLAS that computes the scalar product of two vectors. 
@@ -259,8 +259,8 @@ c
       Complex
      &           cdotc
       Real  
-     &           pscnorm2, pslamch, slapy2
-      external   cdotc, pscnorm2, pslamch, slapy2
+     &           pscnorm2, pslamch10, slapy2
+      external   cdotc, pscnorm2, pslamch10, slapy2
 c
 c     %---------------------%
 c     | Intrinsic Functions |
@@ -298,7 +298,7 @@ c        %---------------------------------%
 c        | Get machine dependent constant. |
 c        %---------------------------------%
 c
-         eps23 = pslamch(comm, 'Epsilon-Machine')
+         eps23 = pslamch10(comm, 'Epsilon-Machine')
          eps23 = eps23**(2.0 / 3.0)
 c
 c        %---------------------------------------%

PARPACK/SRC/MPI/pcnaupd.f

@@ -108,8 +108,8 @@ c  TOL     Real   scalar.  (INPUT)
 c          Stopping criteria: the relative accuracy of the Ritz value 
 c          is considered acceptable if BOUNDS(I) .LE. TOL*ABS(RITZ(I))
 c          where ABS(RITZ(I)) is the magnitude when RITZ(I) is complex.
-c          DEFAULT = pslamch(comm, 'EPS')  (machine precision as computed
-c                    by the ScaLAPACK auxiliary subroutine pslamch).
+c          DEFAULT = pslamch10(comm, 'EPS')  (machine precision as computed
+c                    by the ScaLAPACK auxiliary subroutine pslamch10).
 c
 c  RESID   Complex  array of length N.  (INPUT/OUTPUT)
 c          On INPUT: 
@@ -258,7 +258,7 @@ c          = -7: Length of private work array is not sufficient.
 c          = -8: Error return from LAPACK eigenvalue calculation;
 c          = -9: Starting vector is zero.
 c          = -10: IPARAM(7) must be 1,2,3.
-c          = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatable.
+c          = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatible.
 c          = -12: IPARAM(1) must be equal to 0 or 1.
 c          = -9999: Could not build an Arnoldi factorization.
 c                   User input error highly likely.  Please
@@ -360,7 +360,7 @@ c     cstatn   ARPACK routine that initializes the timing variables.
 c     pivout   Parallel ARPACK utility routine that prints integers.
 c     pcvout   Parallel ARPACK utility routine that prints vectors.
 c     arscnd   ARPACK utility routine for timing.
-c     pslamch  ScaLAPACK routine that determines machine constants.
+c     pslamch10  ScaLAPACK routine that determines machine constants.
 c
 c\Author
 c     Danny Sorensen               Phuong Vu
@@ -389,6 +389,7 @@ c
      &   ( comm, ido, bmat, n, which, nev, tol, resid, ncv, v, ldv, 
      &     iparam, ipntr, workd, workl, lworkl, rwork, info )
 c
+      include   'pcontext.h'
       include  'mpif.h'
 c
 c     %------------------%
@@ -446,21 +447,28 @@ c     %----------------------%
 c     | External Subroutines |
 c     %----------------------%
 c
-      external   pcnaup2, pcvout, pivout, arscnd, cstatn
+      external   pcnaup2, pcvout, pivout, arscnd, cstatn, pcontext
 c
 c     %--------------------%
 c     | External Functions |
 c     %--------------------%
 c
       Real  
-     &           pslamch
-      external   pslamch
+     &           pslamch10
+      external   pslamch10
 c
 c     %-----------------------%
 c     | Executable Statements |
 c     %-----------------------%
 c 
       if (ido .eq. 0) then
+c
+c        %-------------------------------%
+c        | Initialize parallel execution |
+c        | context                       |
+c        %-------------------------------%
+c
+      call pcontext
 c 
 c        %-------------------------------%
 c        | Initialize timing statistics  |
@@ -529,7 +537,7 @@ c        | Set default parameters |
 c        %------------------------%
 c
          if (nb .le. 0)	nb = 1
-         if (tol .le. 0.0  ) tol = pslamch(comm, 'EpsMach')
+         if (tol .le. 0.0  ) tol = pslamch10(comm, 'EpsMach')
          if (ishift .ne. 0  .and.  
      &       ishift .ne. 1  .and.
      &       ishift .ne. 2)	ishift = 1

PARPACK/SRC/MPI/pcneupd.f

@@ -213,7 +213,7 @@ c             in upper triangular form.
 c     ctrsen  LAPACK routine that re-orders the Schur form.
 c     cunm2r  LAPACK routine that applies an orthogonal matrix in 
 c             factored form.
-c     pslamch ScaLAPACK routine that determines machine constants.
+c     pslamch10 ScaLAPACK routine that determines machine constants.
 c     ctrmm   Level 3 BLAS matrix times an upper triangular matrix.
 c     cgeru   Level 2 BLAS rank one update to a matrix.
 c     ccopy   Level 1 BLAS that copies one vector to another .
@@ -336,8 +336,8 @@ c     | External Functions |
 c     %--------------------%
 c
       Real
-     &           scnrm2,pslamch,slapy2
-      external   scnrm2,pslamch,slapy2
+     &           scnrm2,pslamch10,slapy2
+      external   scnrm2,pslamch10,slapy2
 c
       Complex
      &           cdotc
@@ -367,7 +367,7 @@ c     %---------------------------------%
 c     | Get machine dependent constant. |
 c     %---------------------------------%
 c
-      eps23 = pslamch(comm, 'Epsilon-Machine')
+      eps23 = pslamch10(comm, 'Epsilon-Machine')
       eps23 = eps23**(2.0 / 3.0)
 c
 c     %-------------------------------%

PARPACK/SRC/MPI/pcontext.f

@@ -0,0 +1,10 @@
+c
+c   Flags for parallel execution
+c   to be cleared on each new execution of parpack
+c   
+c
+      subroutine pcontext
+      include 'pcontext.h'
+      apps_first = .true.
+      aitr_first = .true.
+      end subroutine

PARPACK/SRC/MPI/pcontext.h

@@ -0,0 +1,8 @@
+c
+c   Flags for parallel execution
+c   to be cleared on each new execution of parpack
+c   
+c
+
+      logical    apps_first, aitr_first
+      common/pcntxt/  apps_first, aitr_first

PARPACK/SRC/MPI/pdlamch10.f

@@ -1,12 +1,12 @@
       DOUBLE PRECISION   FUNCTION PDLAMCH10( ICTXT, CMACH )
       include "mpif.h"
 *
-*  -- ScaLAPACK auxilliary routine (version 1.0) --
+*  -- ScaLAPACK auxiliary routine (version 1.0) --
 *     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
 *     and University of California, Berkeley.
 *     February 28, 1995
 *
-*     The name has been changed in order to avoid symbol collision with 
+*     The name has been changed in order to avoid symbol collision with
 *     the newer version of PDLAMCH in Scalapack which can not be used
 *     with MPI context.
 *
@@ -86,8 +86,8 @@
           TEMP = TEMP1
       END IF
 *
-      PDLAMCH = TEMP
+      PDLAMCH10 = TEMP
 *
-*     End of PDLAMCH
+*     End of PDLAMCH10
 *
       END

PARPACK/SRC/MPI/pdnaitr.f

@@ -220,6 +220,7 @@ c
      &    ipntr, workd, workl, info)
 c
       include   'mpif.h'
+      include   'pcontext.h'
 c
 c     %---------------%
 c     | MPI Variables |
@@ -310,13 +311,12 @@ c     %-----------------%
 c     | Data statements |
 c     %-----------------%
 c
-      data      first / .true. /
 c
 c     %-----------------------%
 c     | Executable Statements |
 c     %-----------------------%
 c
-      if (first) then
+      if (aitr_first) then
 c
 c        %-----------------------------------------%
 c        | Set machine-dependent constants for the |
@@ -401,7 +401,7 @@ c
 c 
 c        %---------------------------------------------------%
 c        | STEP 1: Check if the B norm of j-th residual      |
-c        | vector is zero. Equivalent to determing whether   |
+c        | vector is zero. Equivalent to determining whether |
 c        | an exact j-step Arnoldi factorization is present. |
 c        %---------------------------------------------------%
 c

PARPACK/SRC/MPI/pdnapps.f

@@ -34,7 +34,7 @@ c
 c  KEV     Integer.  (INPUT/OUTPUT)
 c          KEV+NP is the size of the input matrix H.
 c          KEV is the size of the updated matrix HNEW.  KEV is only 
-c          updated on ouput when fewer than NP shifts are applied in
+c          updated on output when fewer than NP shifts are applied in
 c          order to keep the conjugate pair together.
 c
 c  NP      Integer.  (INPUT)
@@ -151,6 +151,7 @@ c     %--------------------%
 c     | MPI Communicator |
 c     %--------------------%
 c
+      include   'pcontext.h'
       integer   comm
 c
 c     %----------------------------------------------------%
@@ -187,11 +188,11 @@ c     | Local Scalars & Arrays |
 c     %------------------------%
 c
       integer    i, iend, ir, istart, j, jj, kplusp, msglvl, nr
-      logical    cconj, first
+      logical    cconj
       Double precision
      &           c, f, g, h11, h12, h21, h22, h32, ovfl, r, s, sigmai, 
      &           sigmar, smlnum, ulp, unfl, u(3), t, tau, tst1
-      save       first, ovfl, smlnum, ulp, unfl 
+      save       ovfl, smlnum, ulp, unfl 
 c
 c     %----------------------%
 c     | External Subroutines |
@@ -215,16 +216,15 @@ c
       intrinsic  abs, max, min
 c
 c     %----------------%
-c     | Data statments |
+c     | Data statements |
 c     %----------------%
 c
-      data       first / .true. /
 c
 c     %-----------------------%
 c     | Executable Statements |
 c     %-----------------------%
 c
-      if (first) then
+      if (apps_first) then
 c
 c        %-----------------------------------------------%
 c        | Set machine-dependent constants for the       |
@@ -238,7 +238,7 @@ c
          call dlabad( unfl, ovfl )
          ulp = pdlamch10( comm, 'precision' )
          smlnum = unfl*( n / ulp )
-         first = .false.
+         apps_first = .false.
       end if
 c
 c     %-------------------------------%
@@ -395,7 +395,7 @@ c
             do 80 i = istart, iend-1
 c
 c              %-----------------------------------------------------%
-c              | Contruct the plane rotation G to zero out the bulge |
+c              | Construct the plane rotation G to zero out the bulge |
 c              %-----------------------------------------------------%
 c
                call dlartg (f, g, c, s, r)

PARPACK/SRC/MPI/pdnaup2.f

@@ -28,7 +28,7 @@ c          products (involving the operator OP) per Arnoldi iteration.
 c          The logic for adjusting is contained within the current
 c          subroutine.
 c          If ISHIFT=0, NP is the number of shifts the user needs
-c          to provide via reverse comunication. 0 < NP < NCV-NEV.
+c          to provide via reverse communication. 0 < NP < NCV-NEV.
 c          NP may be less than NCV-NEV for two reasons. The first, is
 c          to keep complex conjugate pairs of "wanted" Ritz values 
 c          together. The second, is that a leading block of the current
@@ -718,7 +718,7 @@ c
          if (ishift .eq. 0) then
 c
 c           %-------------------------------------------------------%
-c           | User specified shifts: reverse comminucation to       |
+c           | User specified shifts: reverse communication to       |
 c           | compute the shifts. They are returned in the first    |
 c           | 2*NP locations of WORKL.                              |
 c           %-------------------------------------------------------%

PARPACK/SRC/MPI/pdnaupd.f

@@ -276,7 +276,7 @@ c          = -7: Length of private work array is not sufficient.
 c          = -8: Error return from LAPACK eigenvalue calculation;
 c          = -9: Starting vector is zero.
 c          = -10: IPARAM(7) must be 1,2,3,4.
-c          = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatable.
+c          = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatible.
 c          = -12: IPARAM(1) must be equal to 0 or 1.
 c          = -9999: Could not build an Arnoldi factorization.
 c                   IPARAM(5) returns the size of the current Arnoldi
@@ -414,6 +414,7 @@ c
      &     iparam, ipntr, workd, workl, lworkl, info )
 c
       include  'mpif.h'
+      include   'pcontext.h'
 c
 c     %------------------%
 c     | MPI Variables    |
@@ -468,7 +469,7 @@ c     %----------------------%
 c     | External Subroutines |
 c     %----------------------%
 c
-      external   pdnaup2 , pdvout , pivout, arscnd, dstatn 
+      external   pdnaup2 , pdvout , pivout, arscnd, dstatn, pcontext
 c
 c     %--------------------%
 c     | External Functions |
@@ -483,6 +484,13 @@ c     | Executable Statements |
 c     %-----------------------%
 c 
       if (ido .eq. 0) then
+c
+c        %-------------------------------%
+c        | Initialize parallel execution |
+c        | context                       |
+c        %-------------------------------%
+c
+      call pcontext
 c 
 c        %-------------------------------%
 c        | Initialize timing statistics  |

PARPACK/SRC/MPI/pdsaitr.f

@@ -215,6 +215,7 @@ c
      &    ipntr, workd, workl, info)
 c
       include   'mpif.h'
+      include   'pcontext.h'
 c
 c     %---------------%
 c     | MPI Variables |
@@ -259,7 +260,7 @@ c     %---------------%
 c     | Local Scalars |
 c     %---------------%
 c
-      logical    first, orth1, orth2, rstart, step3, step4
+      logical    orth1, orth2, rstart, step3, step4
       integer    i, ierr, ipj, irj, ivj, iter, itry, j, msglvl, infol,
      &           jj
       Double precision
@@ -297,7 +298,6 @@ c     %-----------------%
 c     | Data statements |
 c     %-----------------%
 c
-      data      first / .true. /
 c
 c
 c     %---------------------%
@@ -310,8 +310,8 @@ c     %-----------------------%
 c     | Executable Statements |
 c     %-----------------------%
 c
-      if (first) then
-         first = .false.
+      if (aitr_first) then
+         aitr_first = .false.
 c
 c        %--------------------------------%
 c        | safmin = safe minimum is such  |
@@ -399,8 +399,8 @@ c
          end if
 c
 c        %---------------------------------------------------------%
-c        | Check for exact zero. Equivalent to determing whether a |
-c        | j-step Arnoldi factorization is present.                |
+c        | Check for exact zero. Equivalent to determining whether |
+c        | a j-step Arnoldi factorization is present.              |
 c        %---------------------------------------------------------%
 c
          if (rnorm .gt. zero) go to 40

PARPACK/SRC/MPI/pdsapps.f

@@ -140,6 +140,8 @@ c     %--------------------%
 c     | MPI Communicator |
 c     %--------------------%
 c
+      include  'mpif.h'
+      include  'pcontext.h'
       integer   comm
 c
 c     %----------------------------------------------------%
@@ -176,7 +178,6 @@ c     | Local Scalars |
 c     %---------------%
 c
       integer    i, iend, istart, itop, j, jj, kplusp, msglvl
-      logical    first
       Double precision
      &           a1, a2, a3, a4, big, c, epsmch, f, g, r, s
       save       epsmch, first
@@ -204,18 +205,17 @@ c
       intrinsic  abs
 c
 c     %----------------%
-c     | Data statments |
+c     | Data statements |
 c     %----------------%
 c
-      data       first / .true. /
 c
 c     %-----------------------%
 c     | Executable Statements |
 c     %-----------------------%
 c
-      if (first) then
+      if (apps_first) then
          epsmch = pdlamch10(comm, 'Epsilon-Machine')
-         first = .false.
+         apps_first = .false.
       end if
       itop = 1
 c
@@ -275,7 +275,7 @@ c
                   call pivout (comm, logfil, 1, i, ndigit, 
      &                 '_sapps: deflation at row/column no.')
                   call pivout (comm, logfil, 1, jj, ndigit, 
-     &                 '_sapps: occured before shift number.')
+     &                 '_sapps: occurred before shift number.')
                   call pdvout (comm, logfil, 1, h(i+1,1), ndigit, 
      &                 '_sapps: the corresponding off diagonal element')
                end if

PARPACK/SRC/MPI/pdsaup2.f

@@ -29,7 +29,7 @@ c          products (involving the operator OP) per Arnoldi iteration.
 c          The logic for adjusting is contained within the current
 c          subroutine.
 c          If ISHIFT=0, NP is the number of shifts the user needs
-c          to provide via reverse comunication. 0 < NP < NCV-NEV.
+c          to provide via reverse communication. 0 < NP < NCV-NEV.
 c          NP may be less than NCV-NEV since a leading block of the current
 c          upper Tridiagonal matrix has split off and contains "unwanted"
 c          Ritz values.

PARPACK/SRC/MPI/pdsaupd.f

@@ -271,9 +271,9 @@ c          = -8: Error return from trid. eigenvalue calculation;
 c                Informatinal error from LAPACK routine dsteqr .
 c          = -9: Starting vector is zero.
 c          = -10: IPARAM(7) must be 1,2,3,4,5.
-c          = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatable.
+c          = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatible.
 c          = -12: IPARAM(1) must be equal to 0 or 1.
-c          = -13: NEV and WHICH = 'BE' are incompatable.
+c          = -13: NEV and WHICH = 'BE' are incompatible.
 c          = -9999: Could not build an Arnoldi factorization.
 c                   IPARAM(5) returns the size of the current Arnoldi
 c                   factorization. The user is advised to check that
@@ -310,7 +310,7 @@ c     also increases the work and storage required to maintain the orthogonal
 c     basis vectors.   The optimal "cross-over" with respect to CPU time
 c     is problem dependent and must be determined empirically.
 c
-c  5. If IPARAM(7) = 2 then in the Reverse commuication interface the user
+c  5. If IPARAM(7) = 2 then in the Reverse communication interface the user
 c     must do the following. When IDO = 1, Y = OP * X is to be computed.
 c     When IPARAM(7) = 2 OP = inv(B)*A. After computing A*X the user
 c     must overwrite X with A*X. Y is then the solution to the linear set
@@ -418,6 +418,7 @@ c
      &     iparam, ipntr, workd, workl, lworkl, info )
 c
       include  'mpif.h'
+      include  'pcontext.h'
 c
 c     %------------------%
 c     | MPI Variables    |
@@ -472,7 +473,7 @@ c     %----------------------%
 c     | External Subroutines |
 c     %----------------------%
 c
-      external   pdsaup2 , pdvout , pivout, arscnd, dstats 
+      external   pdsaup2 , pdvout , pivout, arscnd, dstats, pcontext
 c
 c     %--------------------%
 c     | External Functions |
@@ -489,6 +490,13 @@ c
       if (ido .eq. 0) then
 c
 c        %-------------------------------%
+c        | Initialize parallel execution |
+c        | context                       |
+c        %-------------------------------%
+c
+      call pcontext
+c
+c        %-------------------------------%
 c        | Initialize timing statistics  |
 c        | & message level for debugging |
 c        %-------------------------------%

PARPACK/SRC/MPI/pslamch10.f

@@ -1,7 +1,7 @@
-      REAL               FUNCTION PSLAMCH( ICTXT, CMACH )
+      REAL               FUNCTION PSLAMCH10( ICTXT, CMACH )
 *
       include "mpif.h"
-*  -- ScaLAPACK auxilliary routine (version 1.0) --
+*  -- ScaLAPACK auxiliary routine (version 1.0) --
 *     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
 *     and University of California, Berkeley.
 *     February 28, 1995
@@ -82,8 +82,8 @@
           TEMP = TEMP1
       END IF
 *
-      PSLAMCH = TEMP
+      PSLAMCH10 = TEMP
 *
-*     End of PSLAMCH
+*     End of PSLAMCH10
 *
       END

PARPACK/SRC/MPI/psnaitr.f

@@ -142,7 +142,7 @@ c     arscnd   ARPACK utility routine for timing.
 c     psmout   Parallel ARPACK utility routine that prints matrices
 c     psvout   Parallel ARPACK utility routine that prints vectors.
 c     slabad   LAPACK routine that computes machine constants.
-c     pslamch  ScaLAPACK routine that determines machine constants.
+c     pslamch10  ScaLAPACK routine that determines machine constants.
 c     slascl   LAPACK routine for careful scaling of a matrix.
 c     slanhs   LAPACK routine that computes various norms of a matrix.
 c     sgemv    Level 2 BLAS routine for matrix vector multiplication.
@@ -220,6 +220,7 @@ c
      &    ipntr, workd, workl, info)
 c
       include   'mpif.h'
+      include   'pcontext.h'
 c
 c     %---------------%
 c     | MPI Variables |
@@ -264,13 +265,13 @@ c     %---------------%
 c     | Local Scalars |
 c     %---------------%
 c
-      logical    first, orth1, orth2, rstart, step3, step4
+      logical    orth1, orth2, rstart, step3, step4
       integer    ierr, i, infol, ipj, irj, ivj, iter, itry, j, msglvl,
      &           jj
       Real
      &           betaj, ovfl, temp1, rnorm1, smlnum, tst1, ulp, unfl, 
      &           wnorm
-      save       first, orth1, orth2, rstart, step3, step4,
+      save       orth1, orth2, rstart, step3, step4,
      &           ierr, ipj, irj, ivj, iter, itry, j, msglvl, ovfl,
      &           betaj, rnorm1, smlnum, ulp, unfl, wnorm
 c
@@ -297,8 +298,8 @@ c     | External Functions |
 c     %--------------------%
 c
       Real
-     &           sdot, psnorm2, slanhs, pslamch
-      external   sdot, psnorm2, slanhs, pslamch
+     &           sdot, psnorm2, slanhs, pslamch10
+      external   sdot, psnorm2, slanhs, pslamch10
 c
 c     %---------------------%
 c     | Intrinsic Functions |
@@ -310,13 +311,12 @@ c     %-----------------%
 c     | Data statements |
 c     %-----------------%
 c
-      data      first / .true. /
 c
 c     %-----------------------%
 c     | Executable Statements |
 c     %-----------------------%
 c
-      if (first) then
+      if (aitr_first) then
 c
 c        %-----------------------------------------%
 c        | Set machine-dependent constants for the |
@@ -326,12 +326,12 @@ c        | overflow should not occur.              |
 c        | REFERENCE: LAPACK subroutine slahqr     |
 c        %-----------------------------------------%
 c
-         unfl = pslamch(comm, 'safe minimum' )
+         unfl = pslamch10(comm, 'safe minimum' )
          ovfl = one / unfl
          call slabad( unfl, ovfl )
-         ulp = pslamch( comm, 'precision' )
+         ulp = pslamch10( comm, 'precision' )
          smlnum = unfl*( n / ulp )
-         first = .false.
+         aitr_first = .false.
       end if
 c
       if (ido .eq. 0) then
@@ -401,7 +401,7 @@ c
 c 
 c        %---------------------------------------------------%
 c        | STEP 1: Check if the B norm of j-th residual      |
-c        | vector is zero. Equivalent to determing whether   |
+c        | vector is zero. Equivalent to determining whether   |
 c        | an exact j-step Arnoldi factorization is present. |
 c        %---------------------------------------------------%
 c

PARPACK/SRC/MPI/psnapps.f

@@ -34,7 +34,7 @@ c
 c  KEV     Integer.  (INPUT/OUTPUT)
 c          KEV+NP is the size of the input matrix H.
 c          KEV is the size of the updated matrix HNEW.  KEV is only 
-c          updated on ouput when fewer than NP shifts are applied in
+c          updated on output when fewer than NP shifts are applied in
 c          order to keep the conjugate pair together.
 c
 c  NP      Integer.  (INPUT)
@@ -102,7 +102,7 @@ c
 c\Routines called:
 c     slabad  LAPACK routine that computes machine constants.
 c     slacpy  LAPACK matrix copy routine.
-c     pslamch ScaLAPACK routine that determines machine constants.
+c     pslamch10 ScaLAPACK routine that determines machine constants.
 c     slanhs  LAPACK routine that computes various norms of a matrix.
 c     slapy2  LAPACK routine to compute sqrt(x**2+y**2) carefully.
 c     slarf   LAPACK routine that applies Householder reflection to
@@ -151,6 +151,7 @@ c     %--------------------%
 c     | MPI Communicator |
 c     %--------------------%
 c
+      include   'pcontext.h'
       integer   comm
 c
 c     %----------------------------------------------------%
@@ -187,11 +188,11 @@ c     | Local Scalars & Arrays |
 c     %------------------------%
 c
       integer    i, iend, ir, istart, j, jj, kplusp, msglvl, nr
-      logical    cconj, first
+      logical    cconj
       Real
      &           c, f, g, h11, h12, h21, h22, h32, ovfl, r, s, sigmai, 
      &           sigmar, smlnum, ulp, unfl, u(3), t, tau, tst1
-      save       first, ovfl, smlnum, ulp, unfl 
+      save       ovfl, smlnum, ulp, unfl 
 c
 c     %----------------------%
 c     | External Subroutines |
@@ -205,8 +206,8 @@ c     | External Functions |
 c     %--------------------%
 c
       Real
-     &           pslamch, slanhs, slapy2
-      external   pslamch, slanhs, slapy2
+     &           pslamch10, slanhs, slapy2
+      external   pslamch10, slanhs, slapy2
 c
 c     %----------------------%
 c     | Intrinsics Functions |
@@ -215,16 +216,15 @@ c
       intrinsic  abs, max, min
 c
 c     %----------------%
-c     | Data statments |
+c     | Data statements |
 c     %----------------%
 c
-      data       first / .true. /
 c
 c     %-----------------------%
 c     | Executable Statements |
 c     %-----------------------%
 c
-      if (first) then
+      if (apps_first) then
 c
 c        %-----------------------------------------------%
 c        | Set machine-dependent constants for the       |
@@ -233,12 +233,12 @@ c        | overflow should not occur.                    |
 c        | REFERENCE: LAPACK subroutine slahqr           |
 c        %-----------------------------------------------%
 c
-         unfl = pslamch( comm, 'safe minimum' )
+         unfl = pslamch10( comm, 'safe minimum' )
          ovfl = one / unfl
          call slabad( unfl, ovfl )
-         ulp = pslamch( comm, 'precision' )
+         ulp = pslamch10( comm, 'precision' )
          smlnum = unfl*( n / ulp )
-         first = .false.
+         apps_first = .false.
       end if
 c
 c     %-------------------------------%
@@ -395,7 +395,7 @@ c
             do 80 i = istart, iend-1
 c
 c              %-----------------------------------------------------%
-c              | Contruct the plane rotation G to zero out the bulge |
+c              | Construct the plane rotation G to zero out the bulge |
 c              %-----------------------------------------------------%
 c
                call slartg (f, g, c, s, r)

PARPACK/SRC/MPI/psnaup2.f

@@ -28,7 +28,7 @@ c          products (involving the operator OP) per Arnoldi iteration.
 c          The logic for adjusting is contained within the current
 c          subroutine.
 c          If ISHIFT=0, NP is the number of shifts the user needs
-c          to provide via reverse comunication. 0 < NP < NCV-NEV.
+c          to provide via reverse communication. 0 < NP < NCV-NEV.
 c          NP may be less than NCV-NEV for two reasons. The first, is
 c          to keep complex conjugate pairs of "wanted" Ritz values 
 c          together. The second, is that a leading block of the current
@@ -148,7 +148,7 @@ c     pivout   Parallel ARPACK utility routine that prints integers.
 c     arscnd   ARPACK utility routine for timing.
 c     psmout   Parallel ARPACK utility routine that prints matrices
 c     psvout   ARPACK utility routine that prints vectors.
-c     pslamch  ScaLAPACK routine that determines machine constants.
+c     pslamch10  ScaLAPACK routine that determines machine constants.
 c     slapy2   LAPACK routine to compute sqrt(x**2+y**2) carefully.
 c     scopy    Level 1 BLAS that copies one vector to another .
 c     sdot     Level 1 BLAS that computes the scalar product of two vectors. 
@@ -262,8 +262,8 @@ c     | External Functions |
 c     %--------------------%
 c
       Real 
-     &           sdot, psnorm2, slapy2, pslamch
-      external   sdot, psnorm2, slapy2, pslamch
+     &           sdot, psnorm2, slapy2, pslamch10
+      external   sdot, psnorm2, slapy2, pslamch10
 c
 c     %---------------------%
 c     | Intrinsic Functions |
@@ -285,7 +285,7 @@ c        %-------------------------------------%
 c        | Get the machine dependent constant. |
 c        %-------------------------------------%
 c
-         eps23 = pslamch(comm, 'Epsilon-Machine')
+         eps23 = pslamch10(comm, 'Epsilon-Machine')
          eps23 = eps23**(2.0  / 3.0 )
 c
          nev0   = nev
@@ -718,7 +718,7 @@ c
          if (ishift .eq. 0) then
 c
 c           %-------------------------------------------------------%
-c           | User specified shifts: reverse comminucation to       |
+c           | User specified shifts: reverse communication to       |
 c           | compute the shifts. They are returned in the first    |
 c           | 2*NP locations of WORKL.                              |
 c           %-------------------------------------------------------%

PARPACK/SRC/MPI/psnaupd.f

@@ -276,7 +276,7 @@ c          = -7: Length of private work array is not sufficient.
 c          = -8: Error return from LAPACK eigenvalue calculation;
 c          = -9: Starting vector is zero.
 c          = -10: IPARAM(7) must be 1,2,3,4.
-c          = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatable.
+c          = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatible.
 c          = -12: IPARAM(1) must be equal to 0 or 1.
 c          = -9999: Could not build an Arnoldi factorization.
 c                   IPARAM(5) returns the size of the current Arnoldi
@@ -384,7 +384,7 @@ c              Arnoldi Iteration.
 c     pivout   Parallel ARPACK utility routine that prints integers.
 c     arscnd   ARPACK utility routine for timing.
 c     psvout   Parallel ARPACK utility routine that prints vectors.
-c     pslamch  ScaLAPACK routine that determines machine constants.
+c     pslamch10  ScaLAPACK routine that determines machine constants.
 c
 c\Author
 c     Danny Sorensen               Phuong Vu
@@ -414,6 +414,7 @@ c
      &     iparam, ipntr, workd, workl, lworkl, info )
 c
       include  'mpif.h'
+      include   'pcontext.h'
 c
 c     %------------------%
 c     | MPI Variables    |
@@ -468,21 +469,28 @@ c     %----------------------%
 c     | External Subroutines |
 c     %----------------------%
 c
-      external   psnaup2, psvout, pivout, arscnd, sstatn
+      external   psnaup2, psvout, pivout, arscnd, sstatn, pcontext
 c
 c     %--------------------%
 c     | External Functions |
 c     %--------------------%
 c
       Real 
-     &           pslamch
-      external   pslamch
+     &           pslamch10
+      external   pslamch10
 c
 c     %-----------------------%
 c     | Executable Statements |
 c     %-----------------------%
 c 
       if (ido .eq. 0) then
+c
+c        %-------------------------------%
+c        | Initialize parallel execution |
+c        | context                       |
+c        %-------------------------------%
+c
+      call pcontext
 c 
 c        %-------------------------------%
 c        | Initialize timing statistics  |
@@ -553,7 +561,7 @@ c        | Set default parameters |
 c        %------------------------%
 c
          if (nb .le. 0)	nb = 1
-         if (tol .le. zero)	tol = pslamch(comm, 'EpsMach')
+         if (tol .le. zero)	tol = pslamch10(comm, 'EpsMach')
 c
 c        %----------------------------------------------%
 c        | NP is the number of additional steps to      |

PARPACK/SRC/MPI/psneupd.f

@@ -238,7 +238,7 @@ c             a matrix.
 c     slacpy  LAPACK matrix copy routine.
 c     slahqr  LAPACK routine to compute the real Schur form of an
 c             upper Hessenberg matrix.
-c     pslamch ScaLAPACK routine that determines machine constants.
+c     pslamch10 ScaLAPACK routine that determines machine constants.
 c     slapy2  LAPACK routine to compute sqrt(x**2+y**2) carefully.
 c     slaset  LAPACK matrix initialization routine.
 c     sorm2r  LAPACK routine that applies an orthogonal matrix in 
@@ -388,8 +388,8 @@ c     | External Functions |
 c     %--------------------%
 c
       Real 
-     &           slapy2, snrm2, pslamch
-      external   slapy2, snrm2, pslamch
+     &           slapy2, snrm2, pslamch10
+      external   slapy2, snrm2, pslamch10
 c
 c     %---------------------%
 c     | Intrinsic Functions |
@@ -414,7 +414,7 @@ c     %---------------------------------%
 c     | Get machine dependent constant. |
 c     %---------------------------------%
 c
-      eps23 = pslamch(comm, 'Epsilon-Machine')
+      eps23 = pslamch10(comm, 'Epsilon-Machine')
       eps23 = eps23**(2.0  / 3.0 )
 c
 c     %--------------%

PARPACK/SRC/MPI/pssaitr.f

@@ -139,7 +139,7 @@ c     psgetv0  Parallel ARPACK routine to generate the initial vector.
 c     pivout   Parallel ARPACK utility routine that prints integers.
 c     psmout   Parallel ARPACK utility routine that prints matrices.
 c     psvout   Parallel ARPACK utility routine that prints vectors.
-c     pslamch  ScaLAPACK routine that determines machine constants.
+c     pslamch10  ScaLAPACK routine that determines machine constants.
 c     slascl   LAPACK routine for careful scaling of a matrix.
 c     sgemv    Level 2 BLAS routine for matrix vector multiplication.
 c     saxpy    Level 1 BLAS that computes a vector triad.
@@ -214,6 +214,7 @@ c
      &   (comm, ido, bmat, n, k, np, mode, resid, rnorm, v, ldv, h, ldh, 
      &    ipntr, workd, workl, info)
 c
+      include   'pcontext.h'
       include   'mpif.h'
 c
 c     %---------------%
@@ -259,7 +260,7 @@ c     %---------------%
 c     | Local Scalars |
 c     %---------------%
 c
-      logical    first, orth1, orth2, rstart, step3, step4
+      logical    orth1, orth2, rstart, step3, step4
       integer    i, ierr, ipj, irj, ivj, iter, itry, j, msglvl, infol,
      &           jj
       Real
@@ -290,14 +291,13 @@ c     | External Functions |
 c     %--------------------%
 c
       Real
-     &           sdot, psnorm2, pslamch
-      external   sdot, psnorm2, pslamch
+     &           sdot, psnorm2, pslamch10
+      external   sdot, psnorm2, pslamch10
 c
 c     %-----------------%
 c     | Data statements |
 c     %-----------------%
 c
-      data      first / .true. /
 c
 c
 c     %---------------------%
@@ -310,15 +310,15 @@ c     %-----------------------%
 c     | Executable Statements |
 c     %-----------------------%
 c
-      if (first) then
-         first = .false.
+      if (aitr_first) then
+         aitr_first = .false.
 c
 c        %--------------------------------%
 c        | safmin = safe minimum is such  |
 c        | that 1/sfmin does not overflow |
 c        %--------------------------------%
 c
-         safmin = pslamch(comm,'safmin')
+         safmin = pslamch10(comm,'safmin')
       end if
 c
       if (ido .eq. 0) then
@@ -399,7 +399,7 @@ c
          end if
 c
 c        %---------------------------------------------------------%
-c        | Check for exact zero. Equivalent to determing whether a |
+c        | Check for exact zero. Equivalent to determining whether a |
 c        | j-step Arnoldi factorization is present.                |
 c        %---------------------------------------------------------%
 c

PARPACK/SRC/MPI/pssapps.f

@@ -95,7 +95,7 @@ c\Routines called:
 c     pivout  Parallel ARPACK utility routine that prints integers. 
 c     arscnd  ARPACK utility routine for timing.
 c     psvout  Parallel ARPACK utility routine that prints vectors.
-c     pslamch ScaLAPACK routine that determines machine constants.
+c     pslamch10 ScaLAPACK routine that determines machine constants.
 c     slartg  LAPACK Givens rotation construction routine.
 c     slacpy  LAPACK matrix copy routine.
 c     slaset  LAPACK matrix initialization routine.
@@ -140,6 +140,7 @@ c     %--------------------%
 c     | MPI Communicator |
 c     %--------------------%
 c
+      include   'pcontext.h'
       integer   comm
 c
 c     %----------------------------------------------------%
@@ -176,10 +177,9 @@ c     | Local Scalars |
 c     %---------------%
 c
       integer    i, iend, istart, itop, j, jj, kplusp, msglvl
-      logical    first
       Real
      &           a1, a2, a3, a4, big, c, epsmch, f, g, r, s
-      save       epsmch, first
+      save       epsmch
 c
 c
 c     %----------------------%
@@ -194,8 +194,8 @@ c     | External Functions |
 c     %--------------------%
 c
       Real
-     &           pslamch
-      external   pslamch
+     &           pslamch10
+      external   pslamch10
 c
 c     %----------------------%
 c     | Intrinsics Functions |
@@ -204,18 +204,17 @@ c
       intrinsic  abs
 c
 c     %----------------%
-c     | Data statments |
+c     | Data statements |
 c     %----------------%
 c
-      data       first / .true. /
 c
 c     %-----------------------%
 c     | Executable Statements |
 c     %-----------------------%
 c
-      if (first) then
-         epsmch = pslamch(comm, 'Epsilon-Machine')
-         first = .false.
+      if (apps_first) then
+         epsmch = pslamch10(comm, 'Epsilon-Machine')
+         apps_first = .false.
       end if
       itop = 1
 c
@@ -275,7 +274,7 @@ c
                   call pivout (comm, logfil, 1, i, ndigit, 
      &                 '_sapps: deflation at row/column no.')
                   call pivout (comm, logfil, 1, jj, ndigit, 
-     &                 '_sapps: occured before shift number.')
+     &                 '_sapps: occurred before shift number.')
                   call psvout (comm, logfil, 1, h(i+1,1), ndigit, 
      &                 '_sapps: the corresponding off diagonal element')
                end if

PARPACK/SRC/MPI/pssaup2.f

@@ -29,7 +29,7 @@ c          products (involving the operator OP) per Arnoldi iteration.
 c          The logic for adjusting is contained within the current
 c          subroutine.
 c          If ISHIFT=0, NP is the number of shifts the user needs
-c          to provide via reverse comunication. 0 < NP < NCV-NEV.
+c          to provide via reverse communication. 0 < NP < NCV-NEV.
 c          NP may be less than NCV-NEV since a leading block of the current
 c          upper Tridiagonal matrix has split off and contains "unwanted"
 c          Ritz values.
@@ -155,7 +155,7 @@ c              tridiagonal matrix using the implicit QL or QR method.
 c     pivout   Parallel ARPACK utility routine that prints integers.
 c     arscnd   ARPACK utility routine for timing.
 c     psvout   Parallel ARPACK utility routine that prints vectors.
-c     pslamch  ScaLAPACK routine that determines machine constants.
+c     pslamch10  ScaLAPACK routine that determines machine constants.
 c     scopy    Level 1 BLAS that copies one vector to another.
 c     sdot     Level 1 BLAS that computes the scalar product of two vectors.
 c     psnorm2  Parallel version of Level 1 BLAS that computes the norm of a vector.
@@ -263,8 +263,8 @@ c     | External Functions |
 c     %--------------------%
 c
       Real
-     &           sdot, psnorm2, pslamch
-      external   sdot, psnorm2, pslamch
+     &           sdot, psnorm2, pslamch10
+      external   sdot, psnorm2, pslamch10
 c
 c     %---------------------%
 c     | Intrinsic Functions |
@@ -290,7 +290,7 @@ c        %---------------------------------%
 c        | Set machine dependent constant. |
 c        %---------------------------------%
 c
-         eps23 = pslamch(comm, 'Epsilon-Machine')
+         eps23 = pslamch10(comm, 'Epsilon-Machine')
          eps23 = eps23**(2.0/3.0)
 c
 c        %-------------------------------------%

PARPACK/SRC/MPI/pssaupd.f

@@ -271,9 +271,9 @@ c          = -8: Error return from trid. eigenvalue calculation;
 c                Informatinal error from LAPACK routine ssteqr.
 c          = -9: Starting vector is zero.
 c          = -10: IPARAM(7) must be 1,2,3,4,5.
-c          = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatable.
+c          = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatible.
 c          = -12: IPARAM(1) must be equal to 0 or 1.
-c          = -13: NEV and WHICH = 'BE' are incompatable.
+c          = -13: NEV and WHICH = 'BE' are incompatible.
 c          = -9999: Could not build an Arnoldi factorization.
 c                   IPARAM(5) returns the size of the current Arnoldi
 c                   factorization. The user is advised to check that
@@ -310,7 +310,7 @@ c     also increases the work and storage required to maintain the orthogonal
 c     basis vectors.   The optimal "cross-over" with respect to CPU time
 c     is problem dependent and must be determined empirically.
 c
-c  5. If IPARAM(7) = 2 then in the Reverse commuication interface the user
+c  5. If IPARAM(7) = 2 then in the Reverse communication interface the user
 c     must do the following. When IDO = 1, Y = OP * X is to be computed.
 c     When IPARAM(7) = 2 OP = inv(B)*A. After computing A*X the user
 c     must overwrite X with A*X. Y is then the solution to the linear set
@@ -386,7 +386,7 @@ c              variables.
 c     pivout   Parallel ARPACK utility routine that prints integers.
 c     arscnd   ARPACK utility routine for timing.
 c     psvout   Parallel ARPACK utility routine that prints vectors.
-c     pslamch  ScaLAPACK routine that determines machine constants.
+c     pslamch10  ScaLAPACK routine that determines machine constants.
 c
 c\Authors
 c     Kristi Maschhoff ( Parallel Code )
@@ -418,6 +418,7 @@ c
      &     iparam, ipntr, workd, workl, lworkl, info )
 c
       include  'mpif.h'
+      include   'pcontext.h'
 c
 c     %------------------%
 c     | MPI Variables    |
@@ -472,15 +473,15 @@ c     %----------------------%
 c     | External Subroutines |
 c     %----------------------%
 c
-      external   pssaup2, psvout, pivout, arscnd, sstats
+      external   pssaup2, psvout, pivout, arscnd, sstats, pcontext
 c
 c     %--------------------%
 c     | External Functions |
 c     %--------------------%
 c
       Real 
-     &           pslamch
-      external   pslamch
+     &           pslamch10
+      external   pslamch10
 c
 c     %-----------------------%
 c     | Executable Statements |
@@ -489,6 +490,13 @@ c
       if (ido .eq. 0) then
 c
 c        %-------------------------------%
+c        | Initialize parallel execution |
+c        | context                       |
+c        %-------------------------------%
+c
+      call pcontext
+c
+c        %-------------------------------%
 c        | Initialize timing statistics  |
 c        | & message level for debugging |
 c        %-------------------------------%
@@ -563,7 +571,7 @@ c        | Set default parameters |
 c        %------------------------%
 c
          if (nb .le. 0) nb = 1
-         if (tol .le. zero) tol = pslamch(comm, 'EpsMach')
+         if (tol .le. zero) tol = pslamch10(comm, 'EpsMach')
 c
 c        %----------------------------------------------%
 c        | NP is the number of additional steps to      |

PARPACK/SRC/MPI/psseupd.f

@@ -191,7 +191,7 @@ c     psvout  Parallel ARPACK utility routine that prints vectors.
 c     sgeqr2  LAPACK routine that computes the QR factorization of
 c             a matrix.
 c     slacpy  LAPACK matrix copy routine.
-c     pslamch ScaLAPACK routine that determines machine constants.
+c     pslamch10 ScaLAPACK routine that determines machine constants.
 c     sorm2r  LAPACK routine that applies an orthogonal matrix in
 c             factored form.
 c     ssteqr  LAPACK routine that computes eigenvalues and eigenvectors
@@ -295,8 +295,8 @@ c     | External Functions |
 c     %--------------------%
 c
       Real
-     &           psnorm2, pslamch
-      external   psnorm2, pslamch
+     &           psnorm2, pslamch10
+      external   psnorm2, pslamch10
 c
 c     %---------------------%
 c     | Intrinsic Functions |
@@ -439,7 +439,7 @@ c     %---------------------------------%
 c     | Set machine dependent constant. |
 c     %---------------------------------%
 c
-      eps23 = pslamch(comm, 'Epsilon-Machine') 
+      eps23 = pslamch10(comm, 'Epsilon-Machine') 
       eps23 = eps23**(2.0 / 3.0)
 c
 c     %---------------------------------------%

PARPACK/SRC/MPI/pznaitr.f

@@ -223,6 +223,7 @@ c
      &   (comm, ido, bmat, n, k, np, nb, resid, rnorm, v, ldv, h, ldh, 
      &    ipntr, workd, workl, info)
 c
+      include   'pcontext.h'
       include   'mpif.h'
 c
 c     %---------------%
@@ -278,7 +279,7 @@ c     %---------------%
 c     | Local Scalars |
 c     %---------------%
 c
-      logical    first, orth1, orth2, rstart, step3, step4
+      logical    orth1, orth2, rstart, step3, step4
       integer    ierr, i, infol, ipj, irj, ivj, iter, itry, j, msglvl,
      &           jj
       Double precision
@@ -287,7 +288,7 @@ c
       Complex*16
      &           cnorm
 c
-      save       first, orth1, orth2, rstart, step3, step4,
+      save       orth1, orth2, rstart, step3, step4,
      &           ierr, ipj, irj, ivj, iter, itry, j, msglvl, ovfl,
      &           betaj, rnorm1, smlnum, ulp, unfl, wnorm
 c
@@ -321,13 +322,12 @@ c     %-----------------%
 c     | Data statements |
 c     %-----------------%
 c
-      data       first / .true. /
 c
 c     %-----------------------%
 c     | Executable Statements |
 c     %-----------------------%
 c
-      if (first) then
+      if (aitr_first) then
 c
 c        %-----------------------------------------%
 c        | Set machine-dependent constants for the |
@@ -342,7 +342,7 @@ c
          call dlabad( unfl, ovfl )
          ulp = pdlamch10( comm, 'precision' )
          smlnum = unfl*( n / ulp )
-         first = .false.
+         aitr_first = .false.
       end if
 c
       if (ido .eq. 0) then

PARPACK/SRC/MPI/pznapps.f

@@ -149,6 +149,7 @@ c     %--------------------%
 c     | MPI Communicator |
 c     %--------------------%
 c
+      include   'pcontext.h'
       integer   comm
 c
 c     %----------------------------------------------------%
@@ -188,12 +189,11 @@ c     | Local Scalars & Arrays |
 c     %------------------------%
 c
       integer    i, iend, istart, j, jj, kplusp, msglvl
-      logical    first
       Complex*16
      &           cdum, f, g, h11, h21, r, s, sigma, t
       Double precision             
      &           c,  ovfl, smlnum, ulp, unfl, tst1
-      save       first, ovfl, smlnum, ulp, unfl 
+      save       ovfl, smlnum, ulp, unfl 
 c
 c     %----------------------%
 c     | External Subroutines |
@@ -225,16 +225,15 @@ c
       cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( cdum ) )
 c
 c     %----------------%
-c     | Data statments |
+c     | Data statements |
 c     %----------------%
 c
-      data       first / .true. /
 c
 c     %-----------------------%
 c     | Executable Statements |
 c     %-----------------------%
 c
-      if (first) then
+      if (apps_first) then
 c
 c        %-----------------------------------------------%
 c        | Set machine-dependent constants for the       |
@@ -248,7 +247,7 @@ c
          call dlabad( unfl, ovfl )
          ulp = pdlamch10( comm, 'precision' )
          smlnum = unfl*( n / ulp )
-         first = .false.
+         apps_first = .false.
       end if
 c
 c     %-------------------------------%

PARPACK/SRC/MPI/pznaup2.f

@@ -28,7 +28,7 @@ c          products (involving the operator OP) per Arnoldi iteration.
 c          The logic for adjusting is contained within the current
 c          subroutine.
 c          If ISHIFT=0, NP is the number of shifts the user needs
-c          to provide via reverse comunication. 0 < NP < NCV-NEV.
+c          to provide via reverse communication. 0 < NP < NCV-NEV.
 c          NP may be less than NCV-NEV since a leading block of the current
 c          upper Hessenberg matrix has split off and contains "unwanted"
 c          Ritz values.

PARPACK/SRC/MPI/pznaupd.f

@@ -258,7 +258,7 @@ c          = -7: Length of private work array is not sufficient.
 c          = -8: Error return from LAPACK eigenvalue calculation;
 c          = -9: Starting vector is zero.
 c          = -10: IPARAM(7) must be 1,2,3.
-c          = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatable.
+c          = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatible.
 c          = -12: IPARAM(1) must be equal to 0 or 1.
 c          = -9999: Could not build an Arnoldi factorization.
 c                   User input error highly likely.  Please
@@ -389,6 +389,7 @@ c
      &   ( comm, ido, bmat, n, which, nev, tol, resid, ncv, v, ldv, 
      &     iparam, ipntr, workd, workl, lworkl, rwork, info )
 c
+      include  'pcontext.h'
       include  'mpif.h'
 c
 c     %------------------%
@@ -446,7 +447,7 @@ c     %----------------------%
 c     | External Subroutines |
 c     %----------------------%
 c
-      external   pznaup2 , pzvout , pivout, arscnd, zstatn 
+      external   pznaup2 , pzvout , pivout, arscnd, zstatn, pcontext
 c
 c     %--------------------%
 c     | External Functions |
@@ -461,6 +462,13 @@ c     | Executable Statements |
 c     %-----------------------%
 c 
       if (ido .eq. 0) then
+c
+c        %-------------------------------%
+c        | Initialize parallel execution |
+c        | context                       |
+c        %-------------------------------%
+c
+      call pcontext
 c 
 c        %-------------------------------%
 c        | Initialize timing statistics  |

PARPACK/SRC/Makefile.am

@@ -1 +0,0 @@
-SUBDIRS = MPI BLACS

PARPACK/TESTS/MPI/Makefile.am

@@ -0,0 +1,30 @@
+F77 = $(MPIF77)
+LDADD = $(top_builddir)/PARPACK/SRC/MPI/libparpack$(LIBSUFFIX).la $(LAPACK_LIBS) $(BLAS_LIBS)
+
+# Run MPI tests with "mpirun -n 2"
+LOG_COMPILER = mpirun
+LOG_FLAGS = -n 2
+
+SISS = issue46
+
+check_PROGRAMS = $(SISS)
+if ICB
+check_PROGRAMS += icb_parpack_c
+check_PROGRAMS += icb_parpack_cpp
+endif
+
+EXTRA_DIST = debug.h stat.h
+
+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_cpp_SOURCES = icb_parpack_cpp.cpp
+icb_parpack_cpp_LDADD = $(LDADD) $(MPI_CXX_LIBS)
+icb_parpack_cpp_CPPFLAGS = -I$(top_builddir)
+endif
+
+TESTS = $(check_PROGRAMS)

PARPACK/TESTS/MPI/debug.h

@@ -0,0 +1,16 @@
+c
+c\SCCS Information: @(#) 
+c FILE: debug.h   SID: 2.3   DATE OF SID: 11/16/95   RELEASE: 2 
+c
+c     %---------------------------------%
+c     | See debug.doc for documentation |
+c     %---------------------------------%
+      integer  logfil, ndigit, mgetv0,
+     &         msaupd, msaup2, msaitr, mseigt, msapps, msgets, mseupd,
+     &         mnaupd, mnaup2, mnaitr, mneigh, mnapps, mngets, mneupd,
+     &         mcaupd, mcaup2, mcaitr, mceigh, mcapps, mcgets, mceupd
+      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

PARPACK/TESTS/MPI/icb_parpack_c.c

@@ -0,0 +1,202 @@
+/*
+ * This example demonstrates the use of ISO_C_BINDING to call arpack (portability).
+ * IMPORTANT: MPI communicators MUST be passed from C to Fortran using MPI_Comm_c2f.
+ *
+ * Just use arpack as you would have normally done, but, use *[ae]upd_c instead of *[ae]upd_.
+ * The main advantage is that compiler checks (arguments) are performed at build time.
+ * Note: to debug parpack, call debug_c.
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <stdbool.h> // bool.
+#include "mpi.h"
+#include "parpack.h"
+#include <complex.h> // creal, cimag.
+#include "debug_c.h" // debug parpack.
+#include "stat_c.h" // arpack statistics.
+
+/* test program to solve for the 9 largest eigenvalues of
+ * 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;
+  for ( i = 0; i < 1000; ++i)
+    y[i] = ((double) (i+1))*x[i];
+};
+
+int ds() {
+  BLASINT ido = 0;
+  char bmat[] = "I";
+  BLASINT N = 1000;
+  char which[] = "LM";
+  BLASINT nev = 3;
+  double tol = 0;
+  double resid[N];
+  BLASINT ncv = 2*nev+1;
+  double V[ncv*N];
+  BLASINT ldv = N;
+  BLASINT iparam[11];
+  BLASINT ipntr[14];
+  double workd[3*N];
+  bool rvec = true;
+  char howmny[] = "A";
+  double* d = (double*) malloc((nev+1)*sizeof(double));
+  int select[ncv];
+  double z[(N+1)*(nev+1)];
+  BLASINT ldz = N+1;
+  double sigma=0;
+  int k;
+  for (k=0; k < 3*N; ++k )
+    workd[k] = 0;
+  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;
+  int rank; MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+
+  iparam[0] = 1;
+  iparam[2] = 10*N;
+  iparam[3] = 1;
+  iparam[4] = 0; // number of ev found by arpack.
+  iparam[6] = 1;
+
+  MPI_Fint MCW = MPI_Comm_c2f(MPI_COMM_WORLD);
+  while(ido != 99) {
+    /* call arpack like you would have, but, use dsaupd_c instead of dsaupd_ */
+    pdsaupd_c(MCW, &ido, bmat, N, which, nev, tol, resid, ncv, V, ldv, iparam, ipntr,
+              workd, workl, lworkl, &info);
+
+    dMatVec(&(workd[ipntr[0]-1]), &(workd[ipntr[1]-1]));
+  }
+  if (iparam[4] != nev) return 1; // check number of ev found by arpack.
+
+  /* call arpack like you would have, but, use dseupd_c instead of dseupd_ */
+  pdseupd_c(MCW, rvec, howmny, select, d, z, ldz, sigma,
+            bmat, N, which, nev, tol, resid, ncv, V, ldv, iparam, ipntr,
+            workd, workl, lworkl, &info);
+  int i;
+  for (i = 0; i < nev; ++i) {
+    printf("rank %d - %f\n", rank, d[i]);
+    if(fabs(d[i] - (double)(1000-(nev-1)+i))>1e-6){
+      free(d);
+      return 1;
+    }
+  }
+  free(d);
+  return 0;
+}
+
+void zMatVec(double _Complex * x, double _Complex * y) {
+  int i;
+  for (i = 0; i < 1000; ++i)
+    y[i] = x[i] * (i+1.0 + _Complex_I * (i+1.0));
+};
+
+int zn() {
+  BLASINT ido = 0;
+  char bmat[] = "I";
+  BLASINT N = 1000;
+  char which[] = "LM";
+  BLASINT nev = 1;
+  double tol = 0;
+  double _Complex resid[N];
+  BLASINT ncv = 2*nev+1;
+  double _Complex V[ncv*N];
+  BLASINT ldv = N;
+  BLASINT iparam[11];
+  BLASINT 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];
+  double _Complex z[(N+1)*(nev+1)];
+  BLASINT ldz = N+1;
+  double sigma=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;
+  double _Complex rwork[ncv];
+  double _Complex workev[2*ncv];
+  BLASINT info = 0;
+  int rank; MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+
+  iparam[0] = 1;
+  iparam[2] = 10*N;
+  iparam[3] = 1;
+  iparam[4] = 0; // number of ev found by arpack.
+  iparam[6] = 1;
+
+  MPI_Fint MCW = MPI_Comm_c2f(MPI_COMM_WORLD);
+  while(ido != 99) {
+    /* call arpack like you would have, but, use znaupd_c instead of znaupd_ */
+    pznaupd_c(MCW, &ido, bmat, N, which, nev, tol, resid, ncv, V, ldv, iparam, ipntr,
+              workd, workl, lworkl, rwork, &info);
+
+    zMatVec(&(workd[ipntr[0]-1]), &(workd[ipntr[1]-1]));
+  }
+  if (iparam[4] != nev) return 1; // check number of ev found by arpack.
+
+  /* call arpack like you would have, but, use zneupd_c instead of zneupd_ */
+  pzneupd_c(MCW, rvec, howmny, select, d, z, ldz, sigma, workev,
+            bmat, N, which, nev, tol, resid, ncv, V, ldv, iparam, ipntr,
+            workd, workl, lworkl, rwork, &info);
+  int i;
+  for (i = 0; i < nev; ++i) {
+    printf("rank %d - %f %f\n", rank, creal(d[i]), cimag(d[i]));
+    if(fabs(creal(d[i]) - (double)(1000-i))>1e-6 || fabs(cimag(d[i]) - (double)(1000-i))>1e-6){
+      free(d);
+      return 1;
+    }
+  }
+  free(d);
+  return 0;
+}
+
+int main() {
+  MPI_Init(NULL, NULL);
+
+  sstats_c();
+  int rc = ds(); // parpack without debug.
+  fflush(stdout);
+  MPI_Barrier(MPI_COMM_WORLD);
+  if (rc != 0) return rc;
+  int nopx_c, nbx_c, nrorth_c, nitref_c, nrstrt_c;
+  float tsaupd_c, tsaup2_c, tsaitr_c, tseigt_c, tsgets_c, tsapps_c, tsconv_c;
+  float tnaupd_c, tnaup2_c, tnaitr_c, tneigt_c, tngets_c, tnapps_c, tnconv_c;
+  float tcaupd_c, tcaup2_c, tcaitr_c, tceigt_c, tcgets_c, tcapps_c, tcconv_c;
+  float tmvopx_c, tmvbx_c, tgetv0_c, titref_c, trvec_c;
+  stat_c(  &nopx_c,    &nbx_c, &nrorth_c, &nitref_c, &nrstrt_c,
+         &tsaupd_c, &tsaup2_c, &tsaitr_c, &tseigt_c, &tsgets_c, &tsapps_c, &tsconv_c,
+         &tnaupd_c, &tnaup2_c, &tnaitr_c, &tneigt_c, &tngets_c, &tnapps_c, &tnconv_c,
+         &tcaupd_c, &tcaup2_c, &tcaitr_c, &tceigt_c, &tcgets_c, &tcapps_c, &tcconv_c,
+         &tmvopx_c,  &tmvbx_c, &tgetv0_c, &titref_c,  &trvec_c);
+  printf("Timers : nopx %d, tmvopx %f - nbx %d, tmvbx %f\n", nopx_c, tmvopx_c, nbx_c, tmvbx_c);
+
+  int rank = 0; MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+  if (rank == 0) printf("------\n");
+
+  debug_c(6, -6, 1,
+          1, 1, 1, 1, 1, 1, 1,
+          1, 1, 1, 1, 1, 1, 1,
+          1, 1, 1, 1, 1, 1, 1); // set debug flags.
+  rc = zn(); // parpack with debug.
+  fflush(stdout);
+  MPI_Barrier(MPI_COMM_WORLD);
+  if (rc != 0) return rc;
+
+  MPI_Finalize();
+  return 0;
+}

PARPACK/TESTS/MPI/icb_parpack_cpp.cpp

@@ -0,0 +1,231 @@
+/*
+ * This example demonstrates the use of ISO_C_BINDING to call parpack
+ * (portability). IMPORTANT: MPI communicators MUST be passed from C to Fortran
+ * using MPI_Comm_c2f.
+ *
+ * Just use arpack as you would have normally done but use [ae]upd instead
+ * of *[ae]upd_. The main advantage is that checks of the arguments are
+ * performed at compile time. Note: to debug parpack, call debug_c.
+ * This test program solves for the 9 largest eigenvalues of
+ * A*x = lambda*x where A is the diagonal matrix
+ * with entries 1000, 999, ... , 2, 1 on the diagonal.
+ */
+
+#include "parpack.hpp"
+
+#include <array>
+#include <cmath>
+#include <iostream>
+#include <vector>
+
+#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];
+  }
+}
+
+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;
+
+  bool rvec = true;
+  float tol = 0.0f;
+  float sigma = 0.0f;
+
+  std::array<BLASINT, 14> ipntr;
+
+  std::vector<float> workd(3 * N, 0.0f);
+  std::vector<float> workl(3 * (ncv * ncv) + 6 * ncv, 0.0f);
+  std::vector<float> V(ncv * N);
+  std::vector<float> d(nev + 1);
+  std::vector<float> z((N + 1) * (nev + 1));
+  std::vector<float> resid(N);
+  std::vector<int> select(ncv);
+
+  BLASINT info = 0;
+
+  int rank;
+  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+
+  std::array<BLASINT, 11> iparam;
+  iparam[0] = 1;
+  iparam[2] = 10 * N;
+  iparam[3] = 1;
+  iparam[4] = 0;  // number of ev found by arpack.
+  iparam[6] = 1;
+
+  MPI_Fint MCW = MPI_Comm_c2f(MPI_COMM_WORLD);
+
+  BLASINT ido = 0;
+
+  while (ido != 99) {
+    arpack::saupd(MCW, ido, arpack::bmat::identity, N,
+                  arpack::which::largest_magnitude, nev, tol, resid.data(), ncv,
+                  V.data(), ldv, iparam.data(), ipntr.data(), workd.data(),
+                  workl.data(), lworkl, info);
+
+    diagonal_matrix_vector_product(&(workd[ipntr[0] - 1]),
+                                   &(workd[ipntr[1] - 1]));
+  }
+  // check number of ev found by arpack.
+  if (iparam[4] != nev || info != 0) {
+    throw std::domain_error("Error inside ARPACK routines");
+  }
+
+  arpack::seupd(MCW, rvec, arpack::howmny::ritz_vectors, select.data(),
+                d.data(), z.data(), ldz, sigma, arpack::bmat::identity, N,
+                arpack::which::largest_magnitude, nev, tol, resid.data(), ncv,
+                V.data(), ldv, iparam.data(), ipntr.data(), workd.data(),
+                workl.data(), lworkl, info);
+
+  for (int i = 0; i < nev; ++i) {
+    std::cout << "rank " << rank << " - " << d[i] << std::endl;
+    if (std::abs(d[i] - static_cast<float>(1000 - (nev - 1) + i)) > 1e-1) {
+      throw std::domain_error("Correct eigenvalues not computed");
+    }
+  }
+}
+
+void diagonal_matrix_vector_product(std::complex<float> const* const x,
+                                    std::complex<float>* const y) {
+  for (int i = 0; i < 1000; ++i) {
+    y[i] = x[i] * std::complex<float>{i + 1.0f, i + 1.0f};
+  }
+}
+
+void complex_symmetric_runner() {
+  BLASINT N = 1000;
+  BLASINT nev = 1;
+  BLASINT ncv = 2 * nev + 1;
+  BLASINT ldv = N;
+  BLASINT ldz = N + 1;
+
+  float tol = 0.0f;
+  bool rvec = true;
+  float sigma = 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);
+
+  BLASINT 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;
+
+  int rank;
+  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+
+  std::array<BLASINT, 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;
+
+  MPI_Fint MCW = MPI_Comm_c2f(MPI_COMM_WORLD);
+
+  BLASINT ido = 0;
+
+  while (ido != 99) {
+    arpack::naupd(MCW, ido, arpack::bmat::identity, N,
+                  arpack::which::largest_magnitude, nev, tol, resid.data(), ncv,
+                  V.data(), ldv, iparam.data(), ipntr.data(), workd.data(),
+                  workl.data(), lworkl, rwork.data(), info);
+
+    diagonal_matrix_vector_product(&(workd[ipntr[0] - 1]),
+                                   &(workd[ipntr[1] - 1]));
+  }
+
+  // check number of ev found by arpack
+  if (iparam[4] != nev || info != 0) {
+    throw std::domain_error("Error inside ARPACK routines");
+  }
+
+  arpack::neupd(MCW, rvec, arpack::howmny::ritz_vectors, select.data(),
+                d.data(), z.data(), ldz, sigma, workev.data(),
+                arpack::bmat::identity, N, arpack::which::largest_magnitude,
+                nev, tol, resid.data(), ncv, V.data(), ldv, iparam.data(),
+                ipntr.data(), workd.data(), workl.data(), lworkl, rwork.data(),
+                info);
+
+  for (int i = 0; i < nev; ++i) {
+    std::cout << "rank " << rank << " - " << std::real(d[i]) << " "
+              << std::imag(d[i]) << '\n';
+    if (std::abs(std::real(d[i]) - static_cast<float>(1000 - i)) > 1e-1 ||
+        std::abs(std::imag(d[i]) - static_cast<float>(1000 - i)) > 1e-1) {
+      throw std::domain_error("Correct eigenvalues not computed");
+    }
+  }
+}
+
+int main() {
+  MPI_Init(NULL, NULL);
+
+  sstats_c();
+
+  try {
+    // parpack without debug
+    real_symmetric_runner();
+  } catch (std::domain_error& e) {
+    std::cout << e.what() << '\n';
+    MPI_Abort(MPI_COMM_WORLD, 1);
+  }
+
+  MPI_Barrier(MPI_COMM_WORLD);
+
+  int nopx_c, nbx_c, nrorth_c, nitref_c, nrstrt_c;
+  float tsaupd_c, tsaup2_c, tsaitr_c, tseigt_c, tsgets_c, tsapps_c, tsconv_c;
+  float tnaupd_c, tnaup2_c, tnaitr_c, tneigt_c, tngets_c, tnapps_c, tnconv_c;
+  float tcaupd_c, tcaup2_c, tcaitr_c, tceigt_c, tcgets_c, tcapps_c, tcconv_c;
+  float tmvopx_c, tmvbx_c, tgetv0_c, titref_c, trvec_c;
+
+  stat_c(nopx_c, nbx_c, nrorth_c, nitref_c, nrstrt_c, tsaupd_c, tsaup2_c,
+         tsaitr_c, tseigt_c, tsgets_c, tsapps_c, tsconv_c, tnaupd_c, tnaup2_c,
+         tnaitr_c, tneigt_c, tngets_c, tnapps_c, tnconv_c, tcaupd_c, tcaup2_c,
+         tcaitr_c, tceigt_c, tcgets_c, tcapps_c, tcconv_c, tmvopx_c, tmvbx_c,
+         tgetv0_c, titref_c, trvec_c);
+
+  std::cout << "Timers : nopx " << nopx_c << ", tmvopx " << tmvopx_c;
+  std::cout << " - nbx " << nbx_c << ", tmvbx " << tmvbx_c << std::endl;
+
+  int rank = 0;
+  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+  if (rank == 0) std::cout << "------" << std::endl;
+
+  // set debug flags.
+  debug_c(6, -6, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
+          1);
+
+  try {
+    complex_symmetric_runner();
+  } catch (std::domain_error& e) {
+    std::cout << e.what() << '\n';
+    MPI_Abort(MPI_COMM_WORLD, 1);
+  }
+
+  MPI_Barrier(MPI_COMM_WORLD);
+
+  MPI_Finalize();
+
+  return 0;
+}

PARPACK/TESTS/MPI/issue46.f

@@ -0,0 +1,468 @@
+c
+c----------- Example to reproduce issue#46 ----------------------------------------------
+c
+      program issue46
+      include 'mpif.h'
+      integer ierr, comm, color, key, myid, nprocs, cnprocs
+      call MPI_INIT( ierr )
+      call MPI_COMM_RANK( MPI_COMM_WORLD, myid, ierr )
+      call MPI_COMM_SIZE( MPI_COMM_WORLD, nprocs, ierr )
+c     Create commuticator to run arnoldi only on the first CPU
+      if(myid .lt. nprocs/2) then
+      color = 1
+      else
+      color = 0
+      endif
+      write(*,*) "Hello", color
+      call MPI_Comm_split(MPI_COMM_WORLD, color, myid, comm, ierr)
+      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)
+      call MPI_FINALIZE(ierr)
+      end
+
+
+      subroutine parnoldi(comm)
+      include 'mpif.h'
+      include 'debug.h'
+      include 'stat.h'
+
+      integer           comm, myid, nprocs, rc, nloc
+c
+c     %-----------------------------%
+c     | Define leading 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          maxnloc, maxnev, maxncv, ldv
+      parameter       (maxnloc=256, maxnev=10, maxncv=25, 
+     &                 ldv=maxnloc )
+c
+c     %--------------%
+c     | Local Arrays |
+c     %--------------%
+c
+      Double precision
+     &                 v(ldv,maxncv), workl(maxncv*(maxncv+8)),
+     &                 workd(3*maxnloc), d(maxncv,2), resid(maxnloc),
+     &                 ax(maxnloc)
+      logical          select(maxncv)
+      integer          iparam(11), ipntr(11)
+c
+c     %---------------%
+c     | Local Scalars |
+c     %---------------%
+c
+      character        bmat*1, which*2
+      integer          ido, n, nev, ncv, lworkl, info, ierr, j, 
+     &                 nx, nconv, maxitr, mode, ishfts
+      logical          rvec
+      Double precision      
+     &                 tol, sigma
+c
+c     %----------------------------------------------%
+c     | Local Buffers needed for MPI communication |
+c     %----------------------------------------------%
+c
+      Double precision
+     &                  mv_buf(maxnloc)
+c
+c     %------------%
+c     | Parameters |
+c     %------------%
+c
+      Double precision
+     &                 zero
+      parameter        (zero = 0.0)
+c  
+c     %-----------------------------%
+c     | BLAS & LAPACK routines used |
+c     %-----------------------------%
+c
+      Double precision           
+     &                 pdnorm2 
+      external         pdnorm2, daxpy
+c
+c     %---------------------%
+c     | Intrinsic Functions |
+c     %---------------------%
+c
+      intrinsic         abs
+c
+c     %-----------------------%
+c     | Executable Statements |
+c     %-----------------------%
+c
+      call MPI_COMM_RANK( comm, myid, ierr )
+      call MPI_COMM_SIZE( comm, nprocs, ierr )
+c
+      ndigit = -3
+      logfil = 6
+      msaupd = 1
+c
+c     %----------------------------------------------------%
+c     | The number NX is the number of interior points     |
+c     | in the discretization of the 2-dimensional         |
+c     | Laplacian on the unit square with zero Dirichlet   |
+c     | boundary condition.  The number N(=NX*NX) is the   |
+c     | dimension of the matrix.  A standard eigenvalue    |
+c     | problem is solved (BMAT = 'I'). NEV is the number  |
+c     | of eigenvalues to be approximated.  The user can   |
+c     | modify NEV, NCV, WHICH to solve problems of        |
+c     | different sizes, and to get different parts of the |
+c     | spectrum.  However, The following conditions must  |
+c     | be satisfied:                                      |
+c     |                   N <= MAXN,                       | 
+c     |                 NEV <= MAXNEV,                     |
+c     |             NEV + 2 <= NCV <= MAXNCV               | 
+c     %----------------------------------------------------% 
+c
+      nx = 10
+      n = nx*nx
+      nev =  4 
+      ncv =  20 
+c
+c     %--------------------------------------%
+c     | Set up distribution of data to nodes |
+c     %--------------------------------------%
+c
+      nloc = (nx / nprocs)*nx
+      if ( mod(nx, nprocs) .gt. myid ) nloc = nloc + nx
+c
+      if ( nloc .gt. maxnloc ) then
+         print *, ' ERROR with _SDRV1: NLOC is greater than MAXNLOC '
+         go to 9000
+      else if ( nev .gt. maxnev ) then
+         print *, ' ERROR with _SDRV1: NEV is greater than MAXNEV '
+         go to 9000
+      else if ( ncv .gt. maxncv ) then
+         print *, ' ERROR with _SDRV1: NCV is greater than MAXNCV '
+         go to 9000
+      end if
+      bmat = 'I'
+      which = 'SM'
+c
+c     %--------------------------------------------------%
+c     | The work array WORKL is used in PSSAUPD 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 PSSAUPD to start the Arnoldi        |
+c     | iteration.                                       |
+c     %--------------------------------------------------%
+c
+      lworkl = ncv*(ncv+8)
+      tol = zero 
+      info = 0
+      ido = 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 PSSAUPD is used    |
+c     | (IPARAM(7) = 1).  All these options may be        |
+c     | changed by the user. For details, see the         |
+c     | documentation in PSSAUPD.                         |
+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 PSSAUPD and take| 
+c        | actions indicated by parameter IDO until    |
+c        | either convergence is indicated or maxitr   |
+c        | has been exceeded.                          |
+c        %---------------------------------------------%
+c
+         call pdsaupd ( comm, ido, bmat, nloc, 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 |
+c           | here that takes workd(ipntr(1)) as   |
+c           | the input, and return the result to  |
+c           | workd(ipntr(2)).                     |
+c           %--------------------------------------%
+c
+            call av ( comm, nloc, nx, mv_buf, 
+     &               workd(ipntr(1)), workd(ipntr(2)))
+c
+c           %-----------------------------------------%
+c           | L O O P   B A C K to call PSSAUPD again.|
+c           %-----------------------------------------%
+c
+            go to 10
+c
+         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 PSSAUPD.|
+c        %--------------------------%
+c
+         if ( myid .eq. 0 ) then
+            print *, ' '
+            print *, ' Error with _saupd, info = ', info
+            print *, ' Check documentation in _saupd '
+            print *, ' '
+         endif
+c
+      else 
+c
+c        %-------------------------------------------%
+c        | No fatal errors occurred.                 |
+c        | Post-Process using PSSEUPD.               |
+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 pdseupd ( comm, rvec, 'All', select, 
+     &        d, v, ldv, sigma, 
+     &        bmat, nloc, which, nev, tol, resid, ncv, v, ldv, 
+     &        iparam, ipntr, workd, workl, lworkl, ierr )
+c        %----------------------------------------------%
+c        | Eigenvalues are returned in the first column |
+c        | of the two dimensional array D and the       |
+c        | corresponding eigenvectors are returned in   |
+c        | the first NEV columns of the two dimensional |
+c        | array V if requested.  Otherwise, an         |
+c        | orthogonal basis for the invariant subspace  |
+c        | corresponding to the eigenvalues in D is     |
+c        | returned in V.                               |
+c        %----------------------------------------------%
+c
+         if ( ierr .ne. 0) then
+c
+c            %------------------------------------%
+c            | Error condition:                   |
+c            | Check the documentation of PSSEUPD.|
+c            %------------------------------------%
+c
+c
+            if ( myid .eq. 0 ) then
+             	print *, ' '
+             	print *, ' Error with _seupd, info = ', ierr
+             	print *, ' Check the documentation of _seupd. '
+             	print *, ' '
+            endif
+c
+         else
+c
+             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
+                call av(comm, nloc, nx, mv_buf, v(1,j), ax)
+                call daxpy(nloc, -d(j,1), v(1,j), 1, ax, 1)
+                d(j,2) = pdnorm2( comm, nloc, ax, 1 )
+c
+ 20          continue
+c
+c            %-------------------------------%
+c            | Display computed residuals    |
+c            %-------------------------------%
+c
+             call pdmout(comm, 6, nconv, 2, d, maxncv, -6,
+     &            'Ritz values and direct residuals')
+         end if
+c
+c        %------------------------------------------%
+c        | Print additional convergence information |
+c        %------------------------------------------%
+c
+         if (myid .eq. 0)then
+         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 *, '_SDRV1 '
+         print *, '====== '
+         print *, ' '
+         print *, ' Size of the matrix is ', n
+         print *, ' The number of processors is ', nprocs
+         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 *, ' '
+         endif
+c
+      end if
+c
+c     %---------------------------%
+c     | Done with program pdsdrv1.|
+c     %---------------------------%
+c
+ 9000 continue
+c
+c
+      end
+c 
+c ------------------------------------------------------------------
+c     parallel matrix vector subroutine
+c
+c     The matrix used is the 2 dimensional discrete Laplacian on unit
+c     square with zero Dirichlet boundary condition.
+c
+c     Computes w <--- OP*v, where OP is the nx*nx by nx*nx block
+c     tridiagonal matrix
+c
+c                  | T -I          |
+c                  |-I  T -I       |
+c             OP = |   -I  T       |
+c                  |        ...  -I|
+c                  |           -I T|
+c
+c     The subroutine TV is called to computed y<---T*x.
+c-------------------------------------------------------------------
+c
+      subroutine av (comm, nloc, nx, mv_buf, v, w)
+c
+c     .. MPI Declarations ...
+      include           'mpif.h'
+      integer           comm, nprocs, myid, ierr,
+     &                  status(MPI_STATUS_SIZE)
+      integer           nloc, nx, np, j, lo, next, prev
+      Double precision
+     &                  v(nloc), w(nloc), mv_buf(nx), one
+      parameter         (one = 1.0 )
+      external          daxpy
+ 
+      call MPI_COMM_RANK( comm, myid, ierr )
+      call MPI_COMM_SIZE( comm, nprocs, ierr )
+c
+      np = nloc/nx
+      call tv(nx,v(1),w(1))
+      call daxpy(nx, -one, v(nx+1), 1, w(1), 1)
+c
+      if ( np .gt. 2) then
+         do 10 j = 2, np-1
+            lo = (j-1)*nx
+            call tv(nx, v(lo+1), w(lo+1))
+            call daxpy(nx, -one, v(lo-nx+1), 1, w(lo+1), 1)
+            call daxpy(nx, -one, v(lo+nx+1), 1, w(lo+1), 1)
+  10     continue
+      end if
+c
+      if ( np .gt. 1) then
+         lo = (np-1)*nx
+         call tv(nx, v(lo+1), w(lo+1))
+         call daxpy(nx, -one, v(lo-nx+1), 1, w(lo+1), 1)
+      end if
+c
+      next = myid + 1
+      prev = myid - 1
+      if ( myid .lt. nprocs-1 ) then
+         call mpi_send( v((np-1)*nx+1), nx, MPI_DOUBLE_PRECISION,
+     &                  next, myid+1, comm, ierr )
+      endif
+      if ( myid .gt. 0 ) then
+         call mpi_recv( mv_buf, nx, MPI_DOUBLE_PRECISION, prev, myid,
+     &                  comm, status, ierr )
+         call daxpy( nx, -one, mv_buf, 1, w(1), 1 )
+      endif
+c
+      if ( myid .gt. 0 ) then
+         call mpi_send( v(1), nx, MPI_DOUBLE_PRECISION,
+     &                  prev, myid-1, comm, ierr )
+      endif
+      if ( myid .lt. nprocs-1 ) then
+         call mpi_recv( mv_buf, nx, MPI_DOUBLE_PRECISION, next, myid,
+     &                  comm, status, ierr )
+         call daxpy( nx, -one, mv_buf, 1, w(lo+1), 1 )
+      endif
+c
+      return
+      end
+c=========================================================================
+      subroutine tv (nx, x, y)
+c
+      integer           nx, j 
+      Double precision
+     &                  x(nx), y(nx), dd, dl, du
+c
+      Double precision
+     &                 one
+      parameter        (one = 1.0 )
+c
+c     Compute the matrix vector multiplication y<---T*x
+c     where T is a nx by nx tridiagonal matrix with DD on the 
+c     diagonal, DL on the subdiagonal, and DU on the superdiagonal.
+c     
+c
+      dd  = 4.0
+      dl  = -one 
+      du  = -one
+c 
+      y(1) =  dd*x(1) + du*x(2)
+      do 10 j = 2,nx-1
+         y(j) = dl*x(j-1) + dd*x(j) + du*x(j+1) 
+ 10   continue 
+      y(nx) =  dl*x(nx-1) + dd*x(nx) 
+      return
+      end

PARPACK/TESTS/MPI/stat.h

@@ -0,0 +1,21 @@
+c     %--------------------------------%
+c     | See stat.doc for documentation |
+c     %--------------------------------%
+c
+c\SCCS Information: @(#) 
+c FILE: stat.h   SID: 2.2   DATE OF SID: 11/16/95   RELEASE: 2 
+c
+      real       t0, t1, t2, t3, t4, t5
+      save       t0, t1, t2, t3, t4, t5
+c
+      integer    nopx, nbx, nrorth, nitref, nrstrt
+      real       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
+      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

PARPACK/UTIL/BLACS/Makefile.am

@@ -5,5 +5,5 @@ DSRC = pdmout.f pdvout.f
 CSRC = pcmout.f pcvout.f
 ZSRC = pzmout.f pzvout.f
 
-noinst_LTLIBRARIES = libparpackutilblacs.la
-libparpackutilblacs_la_SOURCES = $(SRCS) $(SSRC) $(DSRC) $(CSRC) $(ZSRC)
+noinst_LTLIBRARIES = libparpackutil.la
+libparpackutil_la_SOURCES = $(SRCS) $(SSRC) $(DSRC) $(CSRC) $(ZSRC)

PARPACK/UTIL/MPI/Makefile.am

@@ -7,5 +7,5 @@ DSRC = pdmout.f pdvout.f
 CSRC = pcmout.f pcvout.f
 ZSRC = pzmout.f pzvout.f
 
-noinst_LTLIBRARIES = libparpackutilmpi.la
-libparpackutilmpi_la_SOURCES = $(SRCS) $(SSRC) $(DSRC) $(CSRC) $(ZSRC)
+noinst_LTLIBRARIES = libparpackutil.la
+libparpackutil_la_SOURCES = $(SRCS) $(SSRC) $(DSRC) $(CSRC) $(ZSRC)

PARPACK/UTIL/Makefile.am

@@ -1 +0,0 @@
-SUBDIRS = MPI BLACS

PARPACK_CHANGES

@@ -122,7 +122,7 @@ c  SELECT  Logical array of dimension NCV.  (INPUT/WORKSPACE)
 
 9. 06/04/98:
 
-   p[s,d,c,z]larnv.f  The integer arguement IDIST is not declared.
+   p[s,d,c,z]larnv.f  The integer argument IDIST is not declared.
 
 10. 06/04/98:
    In p[s,d]neupd.f, the declaration of p[s,d]norm2 is missing.
@@ -185,7 +185,7 @@ c  SELECT  Logical array of dimension NCV.  (INPUT/WORKSPACE)
 
 
 18. 08/08/2000
-   In _saitr.f and _naitr.f, add a check for NP = 0 at the very begining of
+   In _saitr.f and _naitr.f, add a check for NP = 0 at the very beginning of
    this routine
 
 19. 09/21/2000
@@ -279,7 +279,7 @@ c  SELECT  Logical array of dimension NCV.  (INPUT/WORKSPACE)
 
 26. 04/10/2001
    In all _seupd codes (both serial and parallel) the call to
-   _sgets had an extra arguement:
+   _sgets had an extra argument:
 
             call dsgets (ishift, which       , nev          ,
      &                np    , workl(irz)  , workl(bounds),
@@ -327,7 +327,11 @@ c  SELECT  Logical array of dimension NCV.  (INPUT/WORKSPACE)
 
 33. 10/24/2003 replace 'slamch' with 'pslamch' in p[c,z]napps.f
 
-34. 10/24/2003 removed an extra arguement to the p[s,d]sgets() call.
+34. 10/24/2003 removed an extra argument to the p[s,d]sgets() call.
 
 35. 10/24/2003 p[c,z]naup2.f contain some lines that were incorrectly
     generated by cpp.  These lines have been removed in the new patch.
+36. 09/18/2016 p*apps.f and p*aitr.f contain condition to fetch machine epsilon.
+    When different p*aupd call use communicator with different number of CPU these
+    conditions cause deadlock. Variables inside conditions are moved to global
+    to be reset each first iteration

README.md

@@ -1,6 +1,9 @@
 ARPACK-NG is a collection of Fortran77 subroutines designed to solve large scale
 eigenvalue problems.
 
+<a href="https://travis-ci.org/opencollab/arpack-ng"><img src="https://travis-ci.org/opencollab/arpack-ng.svg"/></a><br/>
+[![Coverage Status](https://coveralls.io/repos/github/opencollab/arpack-ng/badge.svg?branch=master)](https://coveralls.io/github/opencollab/arpack-ng?branch=master)
+
 Important Features:
 
 * Reverse Communication Interface.
@@ -31,9 +34,9 @@ arpack-ng is replacing arpack almost everywhere.
    The directory SRC contains the top level routines including
    the highest level reverse communication interface routines
 
-      ssaupd, dsaupd - symmetric single and double precision
-      snaupd, dnaupd - non-symmetric single and double precision
-      cnaupd, znaupd - complex non-symmetric single and double precision
+* ssaupd, dsaupd - symmetric single and double precision
+* snaupd, dnaupd - non-symmetric single and double precision
+* cnaupd, znaupd - complex non-symmetric single and double precision
 
    The headers of these routines contain full documentation of calling
    sequence and usage.  Additional information is in the DOCUMENTS directory.
@@ -45,13 +48,13 @@ arpack-ng is replacing arpack almost everywhere.
    data types and precisions may be found in the EXAMPLES directory.
    Upon executing the 'ls EXAMPLES' command you should see
 
-      BAND
-      COMPLEX
-      NONSYM
-      README
-      SIMPLE
-      SVD
-      SYM
+* BAND
+* COMPLEX
+* NONSYM
+* README
+* SIMPLE
+* SVD
+* SYM
 
    Example programs for banded, complex, nonsymmetric, symmetric,
    and singular value decomposition may be found in the directories
@@ -67,30 +70,63 @@ 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.
-   Therefor, the classical:
-   $ sh bootstrap
-   $ ./configure
-   $ make
-   $ make install
-   should work as expected.
+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++).
+
+Therefore, the classical commands should work as expected:
+   
+   
+    $ sh bootstrap
+    $ ./configure
+    $ make
+    $ make check
+    $ make install
+
+Furthermore, ARPACK-NG now provides CMake functionality:
+   
+    $ mkdir build
+    $ cd build
+    $ cmake -D EXAMPLES=ON -D MPI=ON -D BUILD_SHARED_LIBS=ON ..
+    $ make
+    $ make install
+   builds everything including examples and parallel support (with MPI).
+
+To build with code coverage:
+   
+    $ mkdir build
+    $ cd build
+    $ cmake -DCOVERALLS=ON -DCMAKE_BUILD_TYPE=Debug ..
+    $ make all check test coveralls
+
+To get iso_c_binding support:
+   
+    $ ./configure --enable-icb
+    $ cmake -D ICB=ON
+   The install will now provide arpack.h/hpp, parpack.h/hpp and friends.
+   Examples of use can be found in ./TESTS and ./PARPACK/TESTS/MPI.
+   A few related links can be found here:
+
+   * http://fortranwiki.org/fortran/show/ISO_C_BINDING
+   * http://fortranwiki.org/fortran/show/Generating+C+Interfaces
+   * https://www.roguewave.com/sites/rw/files/attachments/StandardizedMixedLanguageProgrammingforCandFortran.pdf
 
 5. Within DOCUMENTS directory there are three files
 
-   ex-sym.doc
-   ex-nonsym.doc and
-   ex-complex.doc
+* ex-sym.doc
+* ex-nonsym.doc and
+* ex-complex.doc
 
    for templates on how to invoke the computational modes of ARPACK.
-   Also look in the README file for explanations concerning the
+   Also look in the README.MD file for explanations concerning the
    other documents.
 
-
-   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
-   Allan Cornet     at  allan.cornet@scilab.org
+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,3 +1,5 @@
+AUTOMAKE_OPTIONS = subdir-objects # Needed as debug_init/icb.f90 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 \
        sgetv0.f ssortc.f ssortr.f ssesrt.f sstqrb.f
@@ -12,7 +14,21 @@ CSRC = cnaitr.f cnapps.f cnaup2.f cnaupd.f cneigh.f cneupd.f cngets.f cstatn.f \
 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
+endif
+
 EXTRA_DIST = debug.h stat.h version.h
 
-noinst_LTLIBRARIES = libarpacksrc.la
-libarpacksrc_la_SOURCES = $(SSRC) $(DSRC) $(CSRC) $(ZSRC)
+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_LIBADD = $(top_builddir)/UTIL/libarpackutil.la $(LAPACK_LIBS) $(BLAS_LIBS)
+libarpack@LIBSUFFIX@_la_LDFLAGS = -no-undefined -version-info 2:0

SRC/cgetv0.f

@@ -239,24 +239,22 @@ c        | the generalized problem when B is possibly (singular).   |
 c        %----------------------------------------------------------%
 c
          call arscnd (t2)
-         if (bmat .eq. 'G') then
-            nopx = nopx + 1
-            ipntr(1) = 1
-            ipntr(2) = n + 1
-            call ccopy (n, resid, 1, workd, 1)
-            ido = -1
-            go to 9000
-         end if
+         nopx = nopx + 1
+         ipntr(1) = 1
+         ipntr(2) = n + 1
+         call ccopy (n, resid, 1, workd, 1)
+         ido = -1
+         go to 9000
       end if
 c 
 c     %----------------------------------------%
-c     | Back from computing B*(initial-vector) |
+c     | Back from computing OP*(initial-vector) |
 c     %----------------------------------------%
 c
       if (first) go to 20
 c
 c     %-----------------------------------------------%
-c     | Back from computing B*(orthogonalized-vector) |
+c     | Back from computing OP*(orthogonalized-vector) |
 c     %-----------------------------------------------%
 c
       if (orth)  go to 40
@@ -271,9 +269,9 @@ c     %------------------------------------------------------%
 c
       call arscnd (t2)
       first = .TRUE.
+      call ccopy (n, workd(n+1), 1, resid, 1)
       if (bmat .eq. 'G') then
          nbx = nbx + 1
-         call ccopy (n, workd(n+1), 1, resid, 1)
          ipntr(1) = n + 1
          ipntr(2) = 1
          ido = 2

SRC/cnapps.f

@@ -209,7 +209,7 @@ c
       cabs1( cdum ) = abs( real( cdum ) ) + abs( aimag( cdum ) )
 c
 c     %----------------%
-c     | Data statments |
+c     | Data statements |
 c     %----------------%
 c
       data       first / .true. /

SRC/cnaup2.f

@@ -26,7 +26,7 @@ c          products (involving the operator OP) per Arnoldi iteration.
 c          The logic for adjusting is contained within the current
 c          subroutine.
 c          If ISHIFT=0, NP is the number of shifts the user needs
-c          to provide via reverse comunication. 0 < NP < NCV-NEV.
+c          to provide via reverse communication. 0 < NP < NCV-NEV.
 c          NP may be less than NCV-NEV since a leading block of the current
 c          upper Hessenberg matrix has split off and contains "unwanted"
 c          Ritz values.

SRC/dgetv0.f

@@ -242,14 +242,12 @@ c        | the generalized problem when B is possibly (singular).   |
 c        %----------------------------------------------------------%
 c
          call arscnd (t2)
-         if (bmat .eq. 'G') then
-            nopx = nopx + 1
-            ipntr(1) = 1
-            ipntr(2) = n + 1
-            call dcopy (n, resid, 1, workd, 1)
-            ido = -1
-            go to 9000
-         end if
+         nopx = nopx + 1
+         ipntr(1) = 1
+         ipntr(2) = n + 1
+         call dcopy (n, resid, 1, workd, 1)
+         ido = -1
+         go to 9000
       end if
 c 
 c     %-----------------------------------------%
@@ -259,7 +257,7 @@ c
       if (first) go to 20
 c
 c     %-----------------------------------------------%
-c     | Back from computing B*(orthogonalized-vector) |
+c     | Back from computing OP*(orthogonalized-vector) |
 c     %-----------------------------------------------%
 c
       if (orth)  go to 40
@@ -276,9 +274,9 @@ c     %------------------------------------------------------%
 c
       call arscnd (t2)
       first = .TRUE.
+      call dcopy (n, workd(n+1), 1, resid, 1)
       if (bmat .eq. 'G') then
          nbx = nbx + 1
-         call dcopy (n, workd(n+1), 1, resid, 1)
          ipntr(1) = n + 1
          ipntr(2) = 1
          ido = 2