phynexis-deps/arpack-ng
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

pyarpack: python binding based on Boost.Python.Numpy exposing C++ API. (#238)

Browse Source 
* pyarpack: python binding based on Boost.Python.Numpy exposing C++ API.

* py3.
This commit is contained in:
Franck HOUSSEN
2020-02-26 12:32:59 +01:00
committed by GitHub
parent ad82dcbc0b
commit cfdb5f680f
24 changed files with 2634 additions and 1 deletions
Download Patch File Download Diff File Expand all files Collapse all files
.travis.yml
+11 -1
View File
@@ -178,13 +178,23 @@ jobs:
apt-get -y install build-essential && \
apt-get -y install git gfortran gcc g++ openmpi-bin libopenmpi-dev cmake && \
apt-get -y install libblas-dev liblapack-dev && \
apt-get -y install libeigen3-dev && \
apt-get -y install python3-minimal python3-pip python3-numpy && \
pip3 install numpy && \
apt-get -y install wget && \
wget https://sourceforge.net/projects/boost/files/boost/1.67.0/boost_1_67_0.tar.gz && \
tar -xf boost_1_67_0.tar.gz && cd boost_1_67_0 && \
./bootstrap.sh --with-libraries=python --with-python=/usr/bin/python3 && \
./b2 install && \
apt-get install locate && \
updatedb && \
cd /tmp && \
cd arpack-ng && \
git status && \
git log -2 && \
sed -e 's/mpirun /mpirun --allow-run-as-root --oversubscribe /' -i CMakeLists.txt && \
mkdir -p build && cd build && \
cmake -DEXAMPLES=ON -DMPI=ON -DICB=ON .. && \
cmake -DEXAMPLES=ON -DMPI=ON -DPYTHON3=ON -DBOOST_PYTHON_LIBSUFFIX='36' .. && \
export VERBOSE=1 && \
make all && \
make test && \
CHANGES
+1
View File
@@ -4,6 +4,7 @@ arpack-ng - 3.8.0
* [BUG FIX]: fix 'Unknown CMake command "check_symbol_exists".' when ICB=ON.
[ Franck Houssen ]
* pyarpack: python binding based on Boost.Python.Numpy exposing C++ API.
* [CLEAN] arpackSolver API: more convenient, suppress template parameters when possible.
* [BUG FIX] ICB using rvec/select: rvec/select turned to integer
bool should be, but, is not always supported (depend on compiler, options).
CMakeLists.txt
+72
View File
@@ -18,6 +18,8 @@ endif ()
option(MPI "Enable parallel support" OFF)
option(ICB "Enable support for *[ae]upd_c with ISO_C_BINDING" OFF)
option(ICBEXMM "Enable support for matrix market example based on ICB" OFF)
option(PYTHON3 "Enable python3 support" OFF)
set(BOOST_PYTHON_LIBSUFFIX "" CACHE STRING "suffix to add to custom boost python libs")
#option(SYSTEM_BLAS "Use system BLAS" ON)
#option(SYSTEM_LAPACK "Use system LAPACK" ON)
option(EXAMPLES "Compile ARPACK examples" OFF)
@@ -64,6 +66,14 @@ function(pexamples list_name)
endforeach()
endfunction(pexamples)
if (PYTHON3)
find_package(PythonInterp 3 REQUIRED)
find_package(PythonLibs 3 REQUIRED)
find_package(Boost COMPONENTS python${BOOST_PYTHON_LIBSUFFIX} numpy${BOOST_PYTHON_LIBSUFFIX} REQUIRED)
set(ICBEXMM "ON")
endif ()
if (ICBEXMM)
find_package(Eigen3)
if (NOT EIGEN3_FOUND) # If not found, piggy-back pkg-config files.
@@ -312,7 +322,9 @@ if (MPI)
set_target_properties(parpack PROPERTIES SOVERSION 2)
endif ()
set(PYINT "int32" CACHE STRING "int type to be used in python scripts") # PYINT : int used with python.
if (INTERFACE64)
set(PYINT "int64" CACHE STRING "int type to be used in python scripts") # PYINT : int used with python.
set(INTERFACE64 1)
if ("${CMAKE_Fortran_COMPILER_ID}" MATCHES "GNU")
set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -fdefault-integer-8")
@@ -592,6 +604,61 @@ if(ICB)
add_test(NAME arpackmm_tst WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY} COMMAND ${BASH_PROGRAM} arpackmm.sh)
endif()
if (PYTHON3)
python_add_module(pyarpack ${arpackutil_STAT_SRCS} ${arpacksrc_STAT_SRCS} ${arpacksrc_ICB} ${PROJECT_SOURCE_DIR}/EXAMPLES/PYARPACK/pyarpack.cpp)
set(pyarpack_HDR ${PROJECT_SOURCE_DIR}/ICB ${PROJECT_SOURCE_DIR}/EXAMPLES/MATRIX_MARKET ${PROJECT_SOURCE_DIR}/EXAMPLES/PYARPACK)
target_include_directories(pyarpack PUBLIC ${pyarpack_HDR} ${EIGEN3_INCLUDE_DIR} ${Boost_INCLUDE_DIRS} ${PYTHON_INCLUDE_DIRS})
target_link_libraries(pyarpack ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES} ${Boost_LIBRARIES} ${PYTHON_LIBRARIES})
install(TARGETS pyarpack
ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR}/pyarpack
LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}/pyarpack)
configure_file("${PROJECT_SOURCE_DIR}/EXAMPLES/PYARPACK/pyarpackSparseBiCGDiag.py.in" "${CMAKE_BINARY_DIR}/pyarpackSparseBiCGDiag.py" @ONLY)
add_test(NAME pyarpackSparseBiCGDiag_tst COMMAND ${PYTHON_EXECUTABLE} pyarpackSparseBiCGDiag.py)
set_tests_properties(pyarpackSparseBiCGDiag_tst PROPERTIES ENVIRONMENT PYTHONPATH=${CMAKE_BINARY_DIR}/lib:$ENV{PYTHONPATH})
configure_file("${PROJECT_SOURCE_DIR}/EXAMPLES/PYARPACK/pyarpackSparseBiCGILU.py.in" "${CMAKE_BINARY_DIR}/pyarpackSparseBiCGILU.py" @ONLY)
add_test(NAME pyarpackSparseBiCGILU_tst COMMAND ${PYTHON_EXECUTABLE} pyarpackSparseBiCGILU.py)
set_tests_properties(pyarpackSparseBiCGILU_tst PROPERTIES ENVIRONMENT PYTHONPATH=${CMAKE_BINARY_DIR}/lib:$ENV{PYTHONPATH})
configure_file("${PROJECT_SOURCE_DIR}/EXAMPLES/PYARPACK/pyarpackSparseCGDiag.py.in" "${CMAKE_BINARY_DIR}/pyarpackSparseCGDiag.py" @ONLY)
add_test(NAME pyarpackSparseCGDiag_tst COMMAND ${PYTHON_EXECUTABLE} pyarpackSparseCGDiag.py)
set_tests_properties(pyarpackSparseCGDiag_tst PROPERTIES ENVIRONMENT PYTHONPATH=${CMAKE_BINARY_DIR}/lib:$ENV{PYTHONPATH})
configure_file("${PROJECT_SOURCE_DIR}/EXAMPLES/PYARPACK/pyarpackSparseCGILU.py.in" "${CMAKE_BINARY_DIR}/pyarpackSparseCGILU.py" @ONLY)
add_test(NAME pyarpackSparseCGILU_tst COMMAND ${PYTHON_EXECUTABLE} pyarpackSparseCGILU.py)
set_tests_properties(pyarpackSparseCGILU_tst PROPERTIES ENVIRONMENT PYTHONPATH=${CMAKE_BINARY_DIR}/lib:$ENV{PYTHONPATH})
configure_file("${PROJECT_SOURCE_DIR}/EXAMPLES/PYARPACK/pyarpackSparseLLT.py.in" "${CMAKE_BINARY_DIR}/pyarpackSparseLLT.py" @ONLY)
add_test(NAME pyarpackSparseLLT_tst COMMAND ${PYTHON_EXECUTABLE} pyarpackSparseLLT.py)
set_tests_properties(pyarpackSparseLLT_tst PROPERTIES ENVIRONMENT PYTHONPATH=${CMAKE_BINARY_DIR}/lib:$ENV{PYTHONPATH})
configure_file("${PROJECT_SOURCE_DIR}/EXAMPLES/PYARPACK/pyarpackSparseLDLT.py.in" "${CMAKE_BINARY_DIR}/pyarpackSparseLDLT.py" @ONLY)
add_test(NAME pyarpackSparseLDLT_tst COMMAND ${PYTHON_EXECUTABLE} pyarpackSparseLDLT.py)
set_tests_properties(pyarpackSparseLDLT_tst PROPERTIES ENVIRONMENT PYTHONPATH=${CMAKE_BINARY_DIR}/lib:$ENV{PYTHONPATH})
configure_file("${PROJECT_SOURCE_DIR}/EXAMPLES/PYARPACK/pyarpackSparseLU.py.in" "${CMAKE_BINARY_DIR}/pyarpackSparseLU.py" @ONLY)
add_test(NAME pyarpackSparseLU_tst COMMAND ${PYTHON_EXECUTABLE} pyarpackSparseLU.py)
set_tests_properties(pyarpackSparseLU_tst PROPERTIES ENVIRONMENT PYTHONPATH=${CMAKE_BINARY_DIR}/lib:$ENV{PYTHONPATH})
configure_file("${PROJECT_SOURCE_DIR}/EXAMPLES/PYARPACK/pyarpackSparseQR.py.in" "${CMAKE_BINARY_DIR}/pyarpackSparseQR.py" @ONLY)
add_test(NAME pyarpackSparseQR_tst COMMAND ${PYTHON_EXECUTABLE} pyarpackSparseQR.py)
set_tests_properties(pyarpackSparseQR_tst PROPERTIES ENVIRONMENT PYTHONPATH=${CMAKE_BINARY_DIR}/lib:$ENV{PYTHONPATH})
configure_file("${PROJECT_SOURCE_DIR}/EXAMPLES/PYARPACK/pyarpackDenseLLT.py.in" "${CMAKE_BINARY_DIR}/pyarpackDenseLLT.py" @ONLY)
add_test(NAME pyarpackDenseLLT_tst COMMAND ${PYTHON_EXECUTABLE} pyarpackDenseLLT.py)
set_tests_properties(pyarpackDenseLLT_tst PROPERTIES ENVIRONMENT PYTHONPATH=${CMAKE_BINARY_DIR}/lib:$ENV{PYTHONPATH})
configure_file("${PROJECT_SOURCE_DIR}/EXAMPLES/PYARPACK/pyarpackDenseLDLT.py.in" "${CMAKE_BINARY_DIR}/pyarpackDenseLDLT.py" @ONLY)
add_test(NAME pyarpackDenseLDLT_tst COMMAND ${PYTHON_EXECUTABLE} pyarpackDenseLDLT.py)
set_tests_properties(pyarpackDenseLDLT_tst PROPERTIES ENVIRONMENT PYTHONPATH=${CMAKE_BINARY_DIR}/lib:$ENV{PYTHONPATH})
configure_file("${PROJECT_SOURCE_DIR}/EXAMPLES/PYARPACK/pyarpackDenseLURR.py.in" "${CMAKE_BINARY_DIR}/pyarpackDenseLURR.py" @ONLY)
add_test(NAME pyarpackDenseLURR_tst COMMAND ${PYTHON_EXECUTABLE} pyarpackDenseLURR.py)
set_tests_properties(pyarpackDenseLURR_tst PROPERTIES ENVIRONMENT PYTHONPATH=${CMAKE_BINARY_DIR}/lib:$ENV{PYTHONPATH})
configure_file("${PROJECT_SOURCE_DIR}/EXAMPLES/PYARPACK/pyarpackDenseQRRR.py.in" "${CMAKE_BINARY_DIR}/pyarpackDenseQRRR.py" @ONLY)
add_test(NAME pyarpackDenseQRRR_tst COMMAND ${PYTHON_EXECUTABLE} pyarpackDenseQRRR.py)
set_tests_properties(pyarpackDenseQRRR_tst PROPERTIES ENVIRONMENT PYTHONPATH=${CMAKE_BINARY_DIR}/lib:$ENV{PYTHONPATH})
configure_file("${PROJECT_SOURCE_DIR}/EXAMPLES/PYARPACK/pyarpackDenseLUPP.py.in" "${CMAKE_BINARY_DIR}/pyarpackDenseLUPP.py" @ONLY)
add_test(NAME pyarpackDenseLUPP_tst COMMAND ${PYTHON_EXECUTABLE} pyarpackDenseLUPP.py)
set_tests_properties(pyarpackDenseLUPP_tst PROPERTIES ENVIRONMENT PYTHONPATH=${CMAKE_BINARY_DIR}/lib:$ENV{PYTHONPATH})
configure_file("${PROJECT_SOURCE_DIR}/EXAMPLES/PYARPACK/pyarpackDenseQRPP.py.in" "${CMAKE_BINARY_DIR}/pyarpackDenseQRPP.py" @ONLY)
add_test(NAME pyarpackDenseQRPP_tst COMMAND ${PYTHON_EXECUTABLE} pyarpackDenseQRPP.py)
set_tests_properties(pyarpackDenseQRPP_tst PROPERTIES ENVIRONMENT PYTHONPATH=${CMAKE_BINARY_DIR}/lib:$ENV{PYTHONPATH})
configure_file("${PROJECT_SOURCE_DIR}/EXAMPLES/PYARPACK/pyarpackRestart.py.in" "${CMAKE_BINARY_DIR}/pyarpackRestart.py" @ONLY)
add_test(NAME pyarpackRestart_tst COMMAND ${PYTHON_EXECUTABLE} pyarpackRestart.py)
set_tests_properties(pyarpackRestart_tst PROPERTIES ENVIRONMENT PYTHONPATH=${CMAKE_BINARY_DIR}/lib:$ENV{PYTHONPATH})
endif()
if (MPI)
set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/PARPACK/TESTS/MPI)
@@ -757,3 +824,8 @@ libsummary("LAPACK" "" "${LAPACK_LIBRARIES}")
if (ICBEXMM)
libsummary("EIGEN3" "${EIGEN3_INCLUDE_DIR}" "")
endif()
if (PYTHON3)
libsummary("PYTHON" "${PYTHON_INCLUDE_DIRS}" "${PYTHON_LIBRARIES}")
message(" -- exe: ${PYTHON_EXECUTABLE}")
libsummary("BOOST" "${Boost_INCLUDE_DIRS}" "${Boost_LIBRARIES}")
endif()
EXAMPLES/PYARPACK/README
+55
View File
@@ -0,0 +1,55 @@
pyarpack: python binding based on Boost.Python.Numpy used to expose arpack C++ API
Installation:
-------------
Python3: ~/arpack-ng/build> cmake -DCMAKE_INSTALL_PREFIX=/tmp/local -DPYTHON3=ON -DBOOST_PYTHON_LIBSUFFIX="3" ..
~/arpack-ng/build> make all test
Note: Boost must have been compiled for Python3.
Usage:
------
>> export PYTHONPATH="/tmp/local/lib/pyarpack:${PYTHONPATH}"
>> python
>> import pyarpack
>> help(pyarpack)
You can use sparse or dense matrices, and, play with iterative or direct mode solvers (CG, LU, ...):
1. choose arpack solver with a given mode solver
1.1. if you need to handle sparse matrices
>> from pyarpack import sparseBiCG as pyarpackSlv
1.2. if you need to handle dense matrices
>> from pyarpack import denseBiCG as pyarpackSlv
2. choose arpack data type (float, double, ...)
>> arpackSlv = pyarpackSlv.double()
3. solve the eigen problem
>> arpackSlv.solve(A, B)
4. get eigen values and vectors
>> print(arpackSlv.vec)
>> print(arpackSlv.val)
You can also:
1. restart a solve from the workspace of a previous solve: check out pyarpackRestart.py.in.
2. compute eigen and / or schur vectors.
Note:
1. arpack data type (float, double, ...) must be consistent with A/B numpy dtypes (float32, float64, ...).
2. sparse matrices must be provided in coo format (n, i, j, Mij), that is, as a tuple where:
2.1. n is an integer.
2.2. i, j, Mij are 1 x nnz numpy arrays.
3. dense matrices must be provided in raw format (Mij, rowOrdered), that is, as a tuple where:
3.1. Mij is an n x n numpy array.
3.2. rowOrdered is a boolean (column ordered if False).
4. arpack mode solver are provided by eigen:
4.1. when solver is iterative, A and B can be sparse only.
4.2. when solver is direct, A and B can be sparse or dense.
Examples:
---------
~/arpack-ng> find . -name *.py.in (template files from which python scripts will result)
EXAMPLES/PYARPACK/pyarpack.cpp
+291
View File
@@ -0,0 +1,291 @@
#include <string>
#include <vector>
#include <sstream> // ostringstream.
#include <complex>
#include <boost/python.hpp>
#include <boost/python/numpy.hpp>
#include <boost/python/suite/indexing/vector_indexing_suite.hpp>
#include <Python.h> // PyErr_SetString.
#include <pyarpackItrSolver.hpp>
#include <pyarpackDrtSolver.hpp>
namespace bp = boost::python;
namespace bn = boost::python::numpy;
template<typename RC, typename FD, typename EM, typename SLV>
void exportArpackSparseItr(bp::scope & pySlv, std::string const & dtype) {
// Created nested namespace in module.
pySlv.attr(dtype.c_str()) = bp::class_<pyarpackSparseItrSolver<RC, FD, EM, SLV>>(dtype.c_str(),
"arpack data type (must be consistent with numpy dtype)")
.def ("solve", &pyarpackSparseItrSolver<RC, FD, EM, SLV>::solve,
(bp::arg("A"), bp::arg("B") = bp::tuple()),
"solve standard or generalised eigen problem where A and B must be sparse and provided in coo format: (dimension, row-indice array, column-indice array, matrice-value array) tuple")
.def ("checkEigVec", &pyarpackSparseItrSolver<RC, FD, EM, SLV>::checkEigVec,
(bp::arg("A"), bp::arg("B") = bp::tuple(), bp::arg("diffTol") = 1.e-3),
"check eigen vectors accuracy where A and B must be sparse and provided in coo format: (dimension, row-indice array, column-indice array, matrice-value array) tuple")
ARPACKSOLVERMEMBER(pyarpackSparseItrSolver)
.def_readwrite("slvTol", &pyarpackSparseItrSolver<RC, FD, EM, SLV>::slvTol,
"tolerance of the iterative mode solver - default: 1.e-6")
.def_readwrite("slvMaxIt", &pyarpackSparseItrSolver<RC, FD, EM, SLV>::slvMaxIt,
"maximum number of iterations of the iterative mode solver - default: 100")
.def_readwrite("slvILUDropTol", &pyarpackSparseItrSolver<RC, FD, EM, SLV>::slvILUDropTol,
"drop tolerance of the ILU preconditioner (if any) of the iterative mode solver - default: 1")
.def_readwrite("slvILUFillFactor", &pyarpackSparseItrSolver<RC, FD, EM, SLV>::slvILUFillFactor,
"fill factor of the ILU preconditioner (if any) of the iterative mode solver - default: 2")
;
};
template<typename RC, typename FD, typename EM, typename SLV>
void exportArpackSparseDrt(bp::scope & pySlv, std::string const & dtype) {
// Created nested namespace in module.
pySlv.attr(dtype.c_str()) = bp::class_<pyarpackSparseDrtSolver<RC, FD, EM, SLV>>(dtype.c_str(),
"arpack data type (must be consistent with numpy dtype)")
.def ("solve", &pyarpackSparseDrtSolver<RC, FD, EM, SLV>::solve,
(bp::arg("A"), bp::arg("B") = bp::tuple()),
"solve standard or generalised eigen problem where A and B must be sparse and provided in coo format: (dimension, row-indice array, column-indice array, matrice-value array) tuple")
.def ("checkEigVec", &pyarpackSparseDrtSolver<RC, FD, EM, SLV>::checkEigVec,
(bp::arg("A"), bp::arg("B") = bp::tuple(), bp::arg("diffTol") = 1.e-3),
"check eigen vectors accuracy where A and B must be sparse and provided in coo format: (dimension, row-indice array, column-indice array, matrice-value array) tuple")
ARPACKSOLVERMEMBER(pyarpackSparseDrtSolver)
.def_readwrite("slvPvtThd", &pyarpackSparseDrtSolver<RC, FD, EM, SLV>::slvPvtThd,
"pivoting tolerance of the direct mode solver - default: 1.e-6")
.def_readwrite("slvOffset", &pyarpackSparseDrtSolver<RC, FD, EM, SLV>::slvOffset,
"cholesky offset (LLT, LDLT) of the direct mode solver - default: 0.")
.def_readwrite("slvScale", &pyarpackSparseDrtSolver<RC, FD, EM, SLV>::slvScale,
"cholesky scale (LLT, LDLT) of the direct mode solver - default: 1.")
;
};
template<typename RC, typename FD, typename EM, typename SLV>
void exportArpackDenseDrt(bp::scope & pySlv, std::string const & dtype) {
// Created nested namespace in module.
pySlv.attr(dtype.c_str()) = bp::class_<pyarpackDenseDrtSolver<RC, FD, EM, SLV>>(dtype.c_str(),
"arpack data type (must be consistent with numpy dtype)")
.def ("solve", &pyarpackDenseDrtSolver<RC, FD, EM, SLV>::solve,
(bp::arg("A"), bp::arg("B") = bp::tuple()),
"solve standard or generalised eigen problem where A and B must be dense and provided in raw format: (n-squared matrice-value array, row or column ordered boolean)")
.def ("checkEigVec", &pyarpackDenseDrtSolver<RC, FD, EM, SLV>::checkEigVec,
(bp::arg("A"), bp::arg("B") = bp::tuple(), bp::arg("diffTol") = 1.e-3),
"check eigen vectors accuracy where A and B must be dense and provided in raw format: (n-squared matrice-value array, row or column ordered boolean)")
ARPACKSOLVERMEMBER(pyarpackDenseDrtSolver)
.def_readwrite("slvPvtThd", &pyarpackDenseDrtSolver<RC, FD, EM, SLV>::slvPvtThd,
"pivoting tolerance of the direct mode solver - default: 1.e-6")
.def_readwrite("slvOffset", &pyarpackDenseDrtSolver<RC, FD, EM, SLV>::slvOffset,
"cholesky offset (LLT, LDLT) of the direct mode solver - default: 0.")
.def_readwrite("slvScale", &pyarpackDenseDrtSolver<RC, FD, EM, SLV>::slvScale,
"cholesky scale (LLT, LDLT) of the direct mode solver - default: 1.")
;
};
class sparseBiCGDiag {};
class sparseBiCGILU {};
class sparseCGDiag {};
class sparseCGILU {};
class sparseLLT {};
class sparseLDLT {};
class sparseLU {};
class sparseQR {};
class denseLLT {};
class denseLDLT {};
class denseLURR {};
class denseQRRR {};
class denseLUPP {};
class denseQRPP {};
std::complex<double> EigVecZGetItem(Eigen::Matrix<std::complex<double>, Eigen::Dynamic, 1> & M, int idx) {
if (idx < 0 || idx >= M.size()) {pyarpackThrowError("index out of range"); return std::complex<double>();}
return M[idx];
};
std::string EigVecZToString(EigVecZ const & vec) {
std::ostringstream s;
s << vec;
return s.str();
};
BOOST_PYTHON_MODULE(pyarpack)
{
// Initialize.
bn::initialize();
bp::class_<std::vector<std::complex<double>>>("StdVecZ")
.def(bp::vector_indexing_suite<std::vector<std::complex<double>>>())
;
bp::class_<Eigen::Matrix<std::complex<double>, Eigen::Dynamic, 1>>("EigVecZ")
.def("__getitem__", &EigVecZGetItem)
.def("__str__", &EigVecZToString)
;
bp::class_<std::vector<EigVecZ>>("StdVecEVZ")
.def("__iter__", bp::iterator<std::vector<EigVecZ>>())
.def(bp::vector_indexing_suite<std::vector<EigVecZ>>())
;
// Specify that this module is actually a package.
bp::object package = bp::scope();
package.attr("__path__") = "pyarpack";
// Create python module.
std::string module = "pyarpack";
bp::object pyModule(bp::handle<>(bp::borrowed(PyImport_AddModule(module.c_str()))));
// Create modules.
{
std::string slv = "sparseBiCGDiag";
std::string slvHelp = "arpack internal mode solver (mode > 1): BiCG with diagonal (Jacobi) preconditioner";
bp::scope pySlvBiCGDiag = bp::class_<sparseBiCGDiag>(slv.c_str(), slvHelp.c_str());
exportArpackSparseItr< float , float, EigSMxS, EigSBiCGS>(pySlvBiCGDiag, "float" );
exportArpackSparseItr< double , double, EigSMxD, EigSBiCGD>(pySlvBiCGDiag, "double");
exportArpackSparseItr<std::complex< float>, float, EigSMxC, EigSBiCGC>(pySlvBiCGDiag, "complexFloat" );
exportArpackSparseItr<std::complex<double>, double, EigSMxZ, EigSBiCGZ>(pySlvBiCGDiag, "complexDouble");
}
{
std::string slv = "sparseBiCGILU";
std::string slvHelp = "arpack internal mode solver (mode > 1): BiCG with ILU preconditioner";
bp::scope pySlvBiCGILU = bp::class_<sparseBiCGILU>(slv.c_str(), slvHelp.c_str());
exportArpackSparseItr< float , float, EigSMxS, EigSBiCGILUS>(pySlvBiCGILU, "float" );
exportArpackSparseItr< double , double, EigSMxD, EigSBiCGILUD>(pySlvBiCGILU, "double");
exportArpackSparseItr<std::complex< float>, float, EigSMxC, EigSBiCGILUC>(pySlvBiCGILU, "complexFloat" );
exportArpackSparseItr<std::complex<double>, double, EigSMxZ, EigSBiCGILUZ>(pySlvBiCGILU, "complexDouble");
}
{
std::string slv = "sparseCGDiag";
std::string slvHelp = "arpack internal mode solver (mode > 1): CG with diagonal (Jacobi) preconditioner";
bp::scope pySlvCGDiag = bp::class_<sparseCGDiag>(slv.c_str(), slvHelp.c_str());
exportArpackSparseItr< float , float, EigSMxS, EigSCGS>(pySlvCGDiag, "float" );
exportArpackSparseItr< double , double, EigSMxD, EigSCGD>(pySlvCGDiag, "double");
exportArpackSparseItr<std::complex< float>, float, EigSMxC, EigSCGC>(pySlvCGDiag, "complexFloat" );
exportArpackSparseItr<std::complex<double>, double, EigSMxZ, EigSCGZ>(pySlvCGDiag, "complexDouble");
}
{
std::string slv = "sparseCGILU";
std::string slvHelp = "arpack internal mode solver (mode > 1): CG with ILU preconditioner";
bp::scope pySlvCGILU = bp::class_<sparseCGILU>(slv.c_str(), slvHelp.c_str());
exportArpackSparseItr< float , float, EigSMxS, EigSCGILUS>(pySlvCGILU, "float" );
exportArpackSparseItr< double , double, EigSMxD, EigSCGILUD>(pySlvCGILU, "double");
exportArpackSparseItr<std::complex< float>, float, EigSMxC, EigSCGILUC>(pySlvCGILU, "complexFloat" );
exportArpackSparseItr<std::complex<double>, double, EigSMxZ, EigSCGILUZ>(pySlvCGILU, "complexDouble");
}
{
std::string slv = "sparseLLT";
std::string slvHelp = "arpack internal mode solver (mode > 1): LLT";
bp::scope pySlvLLT = bp::class_<sparseLLT>(slv.c_str(), slvHelp.c_str());
exportArpackSparseDrt< float , float, EigSMxS, EigSLLTS>(pySlvLLT, "float" );
exportArpackSparseDrt< double , double, EigSMxD, EigSLLTD>(pySlvLLT, "double");
exportArpackSparseDrt<std::complex< float>, float, EigSMxC, EigSLLTC>(pySlvLLT, "complexFloat" );
exportArpackSparseDrt<std::complex<double>, double, EigSMxZ, EigSLLTZ>(pySlvLLT, "complexDouble");
}
{
std::string slv = "sparseLDLT";
std::string slvHelp = "arpack internal mode solver (mode > 1): LDLT";
bp::scope pySlvLDLT = bp::class_<sparseLDLT>(slv.c_str(), slvHelp.c_str());
exportArpackSparseDrt< float , float, EigSMxS, EigSLDLTS>(pySlvLDLT, "float" );
exportArpackSparseDrt< double , double, EigSMxD, EigSLDLTD>(pySlvLDLT, "double");
exportArpackSparseDrt<std::complex< float>, float, EigSMxC, EigSLDLTC>(pySlvLDLT, "complexFloat" );
exportArpackSparseDrt<std::complex<double>, double, EigSMxZ, EigSLDLTZ>(pySlvLDLT, "complexDouble");
}
{
std::string slv = "sparseLU";
std::string slvHelp = "arpack internal mode solver (mode > 1): LU";
bp::scope pySlvLU = bp::class_<sparseLU>(slv.c_str(), slvHelp.c_str());
exportArpackSparseDrt< float , float, EigSMxS, EigSLUS>(pySlvLU, "float" );
exportArpackSparseDrt< double , double, EigSMxD, EigSLUD>(pySlvLU, "double");
exportArpackSparseDrt<std::complex< float>, float, EigSMxC, EigSLUC>(pySlvLU, "complexFloat" );
exportArpackSparseDrt<std::complex<double>, double, EigSMxZ, EigSLUZ>(pySlvLU, "complexDouble");
}
{
std::string slv = "sparseQR";
std::string slvHelp = "arpack internal mode solver (mode > 1): QR";
bp::scope pySlvQR = bp::class_<sparseQR>(slv.c_str(), slvHelp.c_str());
exportArpackSparseDrt< float , float, EigSMxS, EigSQRS>(pySlvQR, "float" );
exportArpackSparseDrt< double , double, EigSMxD, EigSQRD>(pySlvQR, "double");
exportArpackSparseDrt<std::complex< float>, float, EigSMxC, EigSQRC>(pySlvQR, "complexFloat" );
exportArpackSparseDrt<std::complex<double>, double, EigSMxZ, EigSQRZ>(pySlvQR, "complexDouble");
}
{
std::string slv = "denseLLT";
std::string slvHelp = "arpack internal mode solver (mode > 1): LLT";
bp::scope pySlvLLT = bp::class_<denseLLT>(slv.c_str(), slvHelp.c_str());
exportArpackDenseDrt< float , float, EigDMxS, EigDLLTS>(pySlvLLT, "float" );
exportArpackDenseDrt< double , double, EigDMxD, EigDLLTD>(pySlvLLT, "double");
exportArpackDenseDrt<std::complex< float>, float, EigDMxC, EigDLLTC>(pySlvLLT, "complexFloat" );
exportArpackDenseDrt<std::complex<double>, double, EigDMxZ, EigDLLTZ>(pySlvLLT, "complexDouble");
}
{
std::string slv = "denseLDLT";
std::string slvHelp = "arpack internal mode solver (mode > 1): LDLT";
bp::scope pySlvLDLT = bp::class_<denseLDLT>(slv.c_str(), slvHelp.c_str());
exportArpackDenseDrt< float , float, EigDMxS, EigDLDLTS>(pySlvLDLT, "float" );
exportArpackDenseDrt< double , double, EigDMxD, EigDLDLTD>(pySlvLDLT, "double");
exportArpackDenseDrt<std::complex< float>, float, EigDMxC, EigDLDLTC>(pySlvLDLT, "complexFloat" );
exportArpackDenseDrt<std::complex<double>, double, EigDMxZ, EigDLDLTZ>(pySlvLDLT, "complexDouble");
}
{
std::string slv = "denseLURR";
std::string slvHelp = "arpack internal mode solver (mode > 1): LU Rank Revealing (slower, more stable)";
bp::scope pySlvLURR = bp::class_<denseLURR>(slv.c_str(), slvHelp.c_str());
exportArpackDenseDrt< float , float, EigDMxS, EigDFLUS>(pySlvLURR, "float" );
exportArpackDenseDrt< double , double, EigDMxD, EigDFLUD>(pySlvLURR, "double");
exportArpackDenseDrt<std::complex< float>, float, EigDMxC, EigDFLUC>(pySlvLURR, "complexFloat" );
exportArpackDenseDrt<std::complex<double>, double, EigDMxZ, EigDFLUZ>(pySlvLURR, "complexDouble");
}
{
std::string slv = "denseQRRR";
std::string slvHelp = "arpack internal mode solver (mode > 1): QR Rank Revealing (slower, more stable)";
bp::scope pySlvQRRR = bp::class_<denseQRRR>(slv.c_str(), slvHelp.c_str());
exportArpackDenseDrt< float , float, EigDMxS, EigDFQRS>(pySlvQRRR, "float" );
exportArpackDenseDrt< double , double, EigDMxD, EigDFQRD>(pySlvQRRR, "double");
exportArpackDenseDrt<std::complex< float>, float, EigDMxC, EigDFQRC>(pySlvQRRR, "complexFloat" );
exportArpackDenseDrt<std::complex<double>, double, EigDMxZ, EigDFQRZ>(pySlvQRRR, "complexDouble");
}
{
std::string slv = "denseLUPP";
std::string slvHelp = "arpack internal mode solver (mode > 1): LU Partial Pivoting (faster, less stable)";
bp::scope pySlvLUPP = bp::class_<denseLUPP>(slv.c_str(), slvHelp.c_str());
exportArpackDenseDrt< float , float, EigDMxS, EigDPLUS>(pySlvLUPP, "float" );
exportArpackDenseDrt< double , double, EigDMxD, EigDPLUD>(pySlvLUPP, "double");
exportArpackDenseDrt<std::complex< float>, float, EigDMxC, EigDPLUC>(pySlvLUPP, "complexFloat" );
exportArpackDenseDrt<std::complex<double>, double, EigDMxZ, EigDPLUZ>(pySlvLUPP, "complexDouble");
}
{
std::string slv = "denseQRPP";
std::string slvHelp = "arpack internal mode solver (mode > 1): QR Partial Pivoting (faster, less stable)";
bp::scope pySlvQPPR = bp::class_<denseQRPP>(slv.c_str(), slvHelp.c_str());
exportArpackDenseDrt< float , float, EigDMxS, EigDPQRS>(pySlvQPPR, "float" );
exportArpackDenseDrt< double , double, EigDMxD, EigDPQRD>(pySlvQPPR, "double");
exportArpackDenseDrt<std::complex< float>, float, EigDMxC, EigDPQRC>(pySlvQPPR, "complexFloat" );
exportArpackDenseDrt<std::complex<double>, double, EigDMxZ, EigDPQRZ>(pySlvQPPR, "complexDouble");
}
}
// Local Variables:
// mode: c++
// c-file-style:"stroustrup"
// show-trailing-whitespace: t
// End:
/* vim: set sw=2 ts=2 et smartindent :*/
EXAMPLES/PYARPACK/pyarpackDenseLDLT.py.in
+119
View File
@@ -0,0 +1,119 @@
#!/usr/bin/env python
from __future__ import print_function
import numpy as np
from pyarpack import denseLDLT as pyarpackSlv
# Build laplacian.
n = 4
Aij = np.array([], dtype='complex128')
for k in range(n):
for l in range(n):
if l == k:
Aij = np.append(Aij, np.complex128(np.complex( 200., 200.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
elif l == k-1:
Aij = np.append(Aij, np.complex128(np.complex(-101., -101.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
elif l == k+1:
Aij = np.append(Aij, np.complex128(np.complex( -99., -99.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
else:
Aij = np.append(Aij, np.complex128(np.complex( 0., 0.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
for idx, val in enumerate(Aij):
print("A[", idx, "] =", val)
A = (Aij, False) # raw format: Aij values, row ordered (or not).
# Get and tune arpack solver.
arpackSlv = pyarpackSlv.complexDouble() # Caution: complexDouble <=> np.array(..., dtype='complex128')
arpackSlv.verbose = 3 # Set to 0 to get a quiet solve.
arpackSlv.debug = 1 # Set to 0 to get a quiet solve.
arpackSlv.nbEV = 1
arpackSlv.nbCV = 2*arpackSlv.nbEV + 1
arpackSlv.mag = 'LM'
arpackSlv.maxIt = 200
arpackSlv.slvOffset = 0.
arpackSlv.slvScale = 1.
arpackSlv.symPb = False
# Solve eigen problem.
rc = arpackSlv.solve(A)
assert rc == 0, "bad solve"
rc = arpackSlv.checkEigVec(A)
assert rc == 0, "bad checkEigVec"
# Print out results (mode selected, eigen vectors, eigen values, ...).
assert arpackSlv.nbEV == len(arpackSlv.val), "bad result"
print("\nresults:\n")
print("mode selected:", arpackSlv.mode)
print("nb iterations:", arpackSlv.nbIt)
print("Reverse Communication Interface time:", arpackSlv.rciTime, "s")
for val, vec in zip(arpackSlv.val, arpackSlv.vec):
print("eigen value:", val)
print("eigen vector:")
print(vec)
#######################################################################################
print("\n##########################################################################\n")
#######################################################################################
# Build laplacian.
n = 8
Aij = np.array([], dtype='complex64')
Bij = np.array([], dtype='complex64')
for k in range(n):
for l in range(n):
if l == k:
Aij = np.append(Aij, np.complex64(np.complex( 200., 200.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.complex64(np.complex( 33.3, 33.3))) # Casting value on append is MANDATORY or C++ won't get the expected type.
elif l == k-1:
Aij = np.append(Aij, np.complex64(np.complex(-101., -101.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.complex64(np.complex( 16.6, 16.6))) # Casting value on append is MANDATORY or C++ won't get the expected type.
elif l == k+1:
Aij = np.append(Aij, np.complex64(np.complex( -99., -99.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.complex64(np.complex( 16.6, 16.6))) # Casting value on append is MANDATORY or C++ won't get the expected type.
else:
Aij = np.append(Aij, np.complex64(np.complex( 0., 0.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.complex64(np.complex( 0., 0.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
for idx, val in enumerate(Aij):
print("A[", idx, "] =", val)
for idx, val in enumerate(Bij):
print("B[", idx, "] =", val)
A = (Aij, False) # raw format: Aij values, row ordered (or not).
B = (Bij, True) # raw format: Bij values, row ordered (or not).
# Get and tune arpack solver.
arpackSlv = pyarpackSlv.complexFloat() # Caution: complexFloat <=> np.array(..., dtype='complex64')
arpackSlv.verbose = 3 # Set to 0 to get a quiet solve.
arpackSlv.debug = 1 # Set to 0 to get a quiet solve.
arpackSlv.nbEV = 2
arpackSlv.nbCV = 2*arpackSlv.nbEV + 1
arpackSlv.mag = 'LM'
arpackSlv.maxIt = 200
arpackSlv.slvOffset = 0.
arpackSlv.slvScale = 1.
arpackSlv.sigmaReal = 1
arpackSlv.symPb = False
# Solve eigen problem.
rc = arpackSlv.solve(A, B)
assert rc == 0, "bad solve"
rc = arpackSlv.checkEigVec(A, B, 1.e-2)
assert rc == 0, "bad checkEigVec"
# Print out results (mode selected, eigen vectors, eigen values, ...).
assert arpackSlv.nbEV == len(arpackSlv.val), "bad result"
print("\nresults:\n")
print("mode selected:", arpackSlv.mode)
print("nb iterations:", arpackSlv.nbIt)
print("Reverse Communication Interface time:", arpackSlv.rciTime, "s")
for val, vec in zip(arpackSlv.val, arpackSlv.vec):
print("eigen value:", val)
print("eigen vector:")
for v in range(n):
print(vec[v])
EXAMPLES/PYARPACK/pyarpackDenseLLT.py.in
+111
View File
@@ -0,0 +1,111 @@
#!/usr/bin/env python
from __future__ import print_function
import numpy as np
from pyarpack import denseLLT as pyarpackSlv
# Build laplacian.
n = 4
Aij = np.array([], dtype='float64')
for k in range(n):
for l in range(n):
if l == k:
Aij = np.append(Aij, np.float64( 200.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
elif l == k-1 or l == k+1:
Aij = np.append(Aij, np.float64(-100.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
else:
Aij = np.append(Aij, np.float64( 0.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
for idx, val in enumerate(Aij):
print("A[", idx, "] =", val)
A = (Aij, False) # raw format: Aij values, row ordered (or not).
# Get and tune arpack solver.
arpackSlv = pyarpackSlv.double() # Caution: double <=> np.array(..., dtype='float64')
arpackSlv.verbose = 3 # Set to 0 to get a quiet solve.
arpackSlv.debug = 1 # Set to 0 to get a quiet solve.
arpackSlv.nbEV = 1
arpackSlv.nbCV = 2*arpackSlv.nbEV + 1
arpackSlv.mag = 'LM'
arpackSlv.maxIt = 200
arpackSlv.slvOffset = 0.
arpackSlv.slvScale = 1.
# Solve eigen problem.
rc = arpackSlv.solve(A)
assert rc == 0, "bad solve"
rc = arpackSlv.checkEigVec(A)
assert rc == 0, "bad checkEigVec"
# Print out results (mode selected, eigen vectors, eigen values, ...).
assert arpackSlv.nbEV == len(arpackSlv.val), "bad result"
print("\nresults:\n")
print("mode selected:", arpackSlv.mode)
print("nb iterations:", arpackSlv.nbIt)
print("Reverse Communication Interface time:", arpackSlv.rciTime, "s")
for val, vec in zip(arpackSlv.val, arpackSlv.vec):
print("eigen value:", val)
print("eigen vector:")
print(vec)
#######################################################################################
print("\n##########################################################################\n")
#######################################################################################
# Build laplacian.
n = 8
Aij = np.array([], dtype='float32')
Bij = np.array([], dtype='float32')
for k in range(n):
for l in range(n):
if l == k:
Aij = np.append(Aij, np.float32( 200.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.float32( 33.3)) # Casting value on append is MANDATORY or C++ won't get the expected type.
elif l == k-1 or l == k+1:
Aij = np.append(Aij, np.float32(-100.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.float32( 16.6)) # Casting value on append is MANDATORY or C++ won't get the expected type.
else:
Aij = np.append(Aij, np.float32( 0.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.float32( 0.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
for idx, val in enumerate(Aij):
print("A[", idx, "] =", val)
for idx, val in enumerate(Bij):
print("B[", idx, "] =", val)
A = (Aij, False) # raw format: Aij values, row ordered (or not).
B = (Bij, True) # raw format: Bij values, row ordered (or not).
# Get and tune arpack solver.
arpackSlv = pyarpackSlv.float() # Caution: float <=> np.array(..., dtype='float32')
arpackSlv.verbose = 3 # Set to 0 to get a quiet solve.
arpackSlv.debug = 1 # Set to 0 to get a quiet solve.
arpackSlv.nbEV = 2
arpackSlv.nbCV = 2*arpackSlv.nbEV + 1
arpackSlv.mag = 'LM'
arpackSlv.maxIt = 200
arpackSlv.slvOffset = 0.
arpackSlv.slvScale = 1.
# Solve eigen problem.
rc = arpackSlv.solve(A, B)
assert rc == 0, "bad solve"
rc = arpackSlv.checkEigVec(A, B, 1.e-2)
assert rc == 0, "bad checkEigVec"
# Print out results (mode selected, eigen vectors, eigen values, ...).
assert arpackSlv.nbEV == len(arpackSlv.val), "bad result"
print("\nresults:\n")
print("mode selected:", arpackSlv.mode)
print("nb iterations:", arpackSlv.nbIt)
print("Reverse Communication Interface time:", arpackSlv.rciTime, "s")
for val, vec in zip(arpackSlv.val, arpackSlv.vec):
print("eigen value:", val)
print("eigen vector:")
for v in range(n):
print(vec[v])
EXAMPLES/PYARPACK/pyarpackDenseLUPP.py.in
+110
View File
@@ -0,0 +1,110 @@
#!/usr/bin/env python
from __future__ import print_function
import numpy as np
from pyarpack import denseLUPP as pyarpackSlv
# Build laplacian.
n = 4
Aij = np.array([], dtype='float64')
for k in range(n):
for l in range(n):
if l == k:
Aij = np.append(Aij, np.float64( 200.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
elif l == k-1 or l == k+1:
Aij = np.append(Aij, np.float64(-100.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
else:
Aij = np.append(Aij, np.float64( 0.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
for idx, val in enumerate(Aij):
print("A[", idx, "] =", val)
A = (Aij, False) # raw format: Aij values, row ordered (or not).
# Get and tune arpack solver.
arpackSlv = pyarpackSlv.double() # Caution: double <=> np.array(..., dtype='float64')
arpackSlv.verbose = 3 # Set to 0 to get a quiet solve.
arpackSlv.debug = 1 # Set to 0 to get a quiet solve.
arpackSlv.nbEV = 1
arpackSlv.nbCV = 2*arpackSlv.nbEV + 1
arpackSlv.mag = 'LM'
arpackSlv.maxIt = 200
arpackSlv.slvPvtThd = 1.e-6
# Solve eigen problem.
rc = arpackSlv.solve(A)
assert rc == 0, "bad solve"
rc = arpackSlv.checkEigVec(A)
assert rc == 0, "bad checkEigVec"
# Print out results (mode selected, eigen vectors, eigen values, ...).
assert arpackSlv.nbEV == len(arpackSlv.val), "bad result"
print("\nresults:\n")
print("mode selected:", arpackSlv.mode)
print("nb iterations:", arpackSlv.nbIt)
print("Reverse Communication Interface time:", arpackSlv.rciTime, "s")
for val, vec in zip(arpackSlv.val, arpackSlv.vec):
print("eigen value:", val)
print("eigen vector:")
print(vec)
#######################################################################################
print("\n##########################################################################\n")
#######################################################################################
# Build laplacian.
n = 8
Aij = np.array([], dtype='float32')
Bij = np.array([], dtype='float32')
for k in range(n):
for l in range(n):
if l == k:
Aij = np.append(Aij, np.float32( 200.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.float32( 33.3)) # Casting value on append is MANDATORY or C++ won't get the expected type.
elif l == k-1 or l == k+1:
Aij = np.append(Aij, np.float32(-100.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.float32( 16.6)) # Casting value on append is MANDATORY or C++ won't get the expected type.
else:
Aij = np.append(Aij, np.float32( 0.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.float32( 0.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
for idx, val in enumerate(Aij):
print("A[", idx, "] =", val)
for idx, val in enumerate(Bij):
print("B[", idx, "] =", val)
A = (Aij, False) # raw format: Aij values, row ordered (or not).
B = (Bij, True) # raw format: Bij values, row ordered (or not).
# Get and tune arpack solver.
arpackSlv = pyarpackSlv.float() # Caution: float <=> np.array(..., dtype='float32')
arpackSlv.verbose = 3 # Set to 0 to get a quiet solve.
arpackSlv.debug = 1 # Set to 0 to get a quiet solve.
arpackSlv.nbEV = 2
arpackSlv.nbCV = 2*arpackSlv.nbEV + 1
arpackSlv.mag = 'LM'
arpackSlv.maxIt = 200
arpackSlv.slvPvtThd = 1.e-6
arpackSlv.sigmaReal = 1
# Solve eigen problem.
rc = arpackSlv.solve(A, B)
assert rc == 0, "bad solve"
rc = arpackSlv.checkEigVec(A, B, 1.e-2)
assert rc == 0, "bad checkEigVec"
# Print out results (mode selected, eigen vectors, eigen values, ...).
assert arpackSlv.nbEV == len(arpackSlv.val), "bad result"
print("\nresults:\n")
print("mode selected:", arpackSlv.mode)
print("nb iterations:", arpackSlv.nbIt)
print("Reverse Communication Interface time:", arpackSlv.rciTime, "s")
for val, vec in zip(arpackSlv.val, arpackSlv.vec):
print("eigen value:", val)
print("eigen vector:")
for v in range(n):
print(vec[v])
EXAMPLES/PYARPACK/pyarpackDenseLURR.py.in
+110
View File
@@ -0,0 +1,110 @@
#!/usr/bin/env python
from __future__ import print_function
import numpy as np
from pyarpack import denseLURR as pyarpackSlv
# Build laplacian.
n = 4
Aij = np.array([], dtype='float64')
for k in range(n):
for l in range(n):
if l == k:
Aij = np.append(Aij, np.float64( 200.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
elif l == k-1 or l == k+1:
Aij = np.append(Aij, np.float64(-100.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
else:
Aij = np.append(Aij, np.float64( 0.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
for idx, val in enumerate(Aij):
print("A[", idx, "] =", val)
A = (Aij, False) # raw format: Aij values, row ordered (or not).
# Get and tune arpack solver.
arpackSlv = pyarpackSlv.double() # Caution: double <=> np.array(..., dtype='float64')
arpackSlv.verbose = 3 # Set to 0 to get a quiet solve.
arpackSlv.debug = 1 # Set to 0 to get a quiet solve.
arpackSlv.nbEV = 1
arpackSlv.nbCV = 2*arpackSlv.nbEV + 1
arpackSlv.mag = 'LM'
arpackSlv.maxIt = 200
arpackSlv.slvPvtThd = 1.e-6
# Solve eigen problem.
rc = arpackSlv.solve(A)
assert rc == 0, "bad solve"
rc = arpackSlv.checkEigVec(A)
assert rc == 0, "bad checkEigVec"
# Print out results (mode selected, eigen vectors, eigen values, ...).
assert arpackSlv.nbEV == len(arpackSlv.val), "bad result"
print("\nresults:\n")
print("mode selected:", arpackSlv.mode)
print("nb iterations:", arpackSlv.nbIt)
print("Reverse Communication Interface time:", arpackSlv.rciTime, "s")
for val, vec in zip(arpackSlv.val, arpackSlv.vec):
print("eigen value:", val)
print("eigen vector:")
print(vec)
#######################################################################################
print("\n##########################################################################\n")
#######################################################################################
# Build laplacian.
n = 8
Aij = np.array([], dtype='float32')
Bij = np.array([], dtype='float32')
for k in range(n):
for l in range(n):
if l == k:
Aij = np.append(Aij, np.float32( 200.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.float32( 33.3)) # Casting value on append is MANDATORY or C++ won't get the expected type.
elif l == k-1 or l == k+1:
Aij = np.append(Aij, np.float32(-100.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.float32( 16.6)) # Casting value on append is MANDATORY or C++ won't get the expected type.
else:
Aij = np.append(Aij, np.float32( 0.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.float32( 0.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
for idx, val in enumerate(Aij):
print("A[", idx, "] =", val)
for idx, val in enumerate(Bij):
print("B[", idx, "] =", val)
A = (Aij, False) # raw format: Aij values, row ordered (or not).
B = (Bij, True) # raw format: Bij values, row ordered (or not).
# Get and tune arpack solver.
arpackSlv = pyarpackSlv.float() # Caution: float <=> np.array(..., dtype='float32')
arpackSlv.verbose = 3 # Set to 0 to get a quiet solve.
arpackSlv.debug = 1 # Set to 0 to get a quiet solve.
arpackSlv.nbEV = 2
arpackSlv.nbCV = 2*arpackSlv.nbEV + 1
arpackSlv.mag = 'LM'
arpackSlv.maxIt = 200
arpackSlv.slvPvtThd = 1.e-6
arpackSlv.sigmaReal = 1
# Solve eigen problem.
rc = arpackSlv.solve(A, B)
assert rc == 0, "bad solve"
rc = arpackSlv.checkEigVec(A, B, 1.e-2)
assert rc == 0, "bad checkEigVec"
# Print out results (mode selected, eigen vectors, eigen values, ...).
assert arpackSlv.nbEV == len(arpackSlv.val), "bad result"
print("\nresults:\n")
print("mode selected:", arpackSlv.mode)
print("nb iterations:", arpackSlv.nbIt)
print("Reverse Communication Interface time:", arpackSlv.rciTime, "s")
for val, vec in zip(arpackSlv.val, arpackSlv.vec):
print("eigen value:", val)
print("eigen vector:")
for v in range(n):
print(vec[v])
EXAMPLES/PYARPACK/pyarpackDenseQRPP.py.in
+117
View File
@@ -0,0 +1,117 @@
#!/usr/bin/env python
from __future__ import print_function
import numpy as np
from pyarpack import denseQRPP as pyarpackSlv
# Build laplacian.
n = 4
Aij = np.array([], dtype='complex128')
for k in range(n):
for l in range(n):
if l == k:
Aij = np.append(Aij, np.complex128(np.complex( 200., 200.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
elif l == k-1:
Aij = np.append(Aij, np.complex128(np.complex(-101., -101.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
elif l == k+1:
Aij = np.append(Aij, np.complex128(np.complex( -99., -99.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
else:
Aij = np.append(Aij, np.complex128(np.complex( 0., 0.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
for idx, val in enumerate(Aij):
print("A[", idx, "] =", val)
A = (Aij, False) # raw format: Aij values, row ordered (or not).
# Get and tune arpack solver.
arpackSlv = pyarpackSlv.complexDouble() # Caution: complexDouble <=> np.array(..., dtype='complex128')
arpackSlv.verbose = 3 # Set to 0 to get a quiet solve.
arpackSlv.debug = 1 # Set to 0 to get a quiet solve.
arpackSlv.nbEV = 1
arpackSlv.nbCV = 2*arpackSlv.nbEV + 1
arpackSlv.mag = 'LM'
arpackSlv.maxIt = 200
arpackSlv.slvPvtThd = 1.e-6
arpackSlv.symPb = False
# Solve eigen problem.
rc = arpackSlv.solve(A)
assert rc == 0, "bad solve"
rc = arpackSlv.checkEigVec(A)
assert rc == 0, "bad checkEigVec"
# Print out results (mode selected, eigen vectors, eigen values, ...).
assert arpackSlv.nbEV == len(arpackSlv.val), "bad result"
print("\nresults:\n")
print("mode selected:", arpackSlv.mode)
print("nb iterations:", arpackSlv.nbIt)
print("Reverse Communication Interface time:", arpackSlv.rciTime, "s")
for val, vec in zip(arpackSlv.val, arpackSlv.vec):
print("eigen value:", val)
print("eigen vector:")
print(vec)
#######################################################################################
print("\n##########################################################################\n")
#######################################################################################
# Build laplacian.
n = 8
Aij = np.array([], dtype='complex64')
Bij = np.array([], dtype='complex64')
for k in range(n):
for l in range(n):
if l == k:
Aij = np.append(Aij, np.complex64(np.complex( 200., 200.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.complex64(np.complex( 33.3, 33.3))) # Casting value on append is MANDATORY or C++ won't get the expected type.
elif l == k-1:
Aij = np.append(Aij, np.complex64(np.complex(-101., -101.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.complex64(np.complex( 16.6, 16.6))) # Casting value on append is MANDATORY or C++ won't get the expected type.
elif l == k+1:
Aij = np.append(Aij, np.complex64(np.complex( -99., -99.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.complex64(np.complex( 16.6, 16.6))) # Casting value on append is MANDATORY or C++ won't get the expected type.
else:
Aij = np.append(Aij, np.complex64(np.complex( 0., 0.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.complex64(np.complex( 0., 0.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
for idx, val in enumerate(Aij):
print("A[", idx, "] =", val)
for idx, val in enumerate(Bij):
print("B[", idx, "] =", val)
A = (Aij, False) # raw format: Aij values, row ordered (or not).
B = (Bij, True) # raw format: Bij values, row ordered (or not).
# Get and tune arpack solver.
arpackSlv = pyarpackSlv.complexFloat() # Caution: complexFloat <=> np.array(..., dtype='complex64')
arpackSlv.verbose = 3 # Set to 0 to get a quiet solve.
arpackSlv.debug = 1 # Set to 0 to get a quiet solve.
arpackSlv.nbEV = 2
arpackSlv.nbCV = 2*arpackSlv.nbEV + 1
arpackSlv.mag = 'LM'
arpackSlv.maxIt = 200
arpackSlv.slvPvtThd = 1.e-6
arpackSlv.sigmaReal = 1
arpackSlv.symPb = False
# Solve eigen problem.
rc = arpackSlv.solve(A, B)
assert rc == 0, "bad solve"
rc = arpackSlv.checkEigVec(A, B, 1.e-2)
assert rc == 0, "bad checkEigVec"
# Print out results (mode selected, eigen vectors, eigen values, ...).
assert arpackSlv.nbEV == len(arpackSlv.val), "bad result"
print("\nresults:\n")
print("mode selected:", arpackSlv.mode)
print("nb iterations:", arpackSlv.nbIt)
print("Reverse Communication Interface time:", arpackSlv.rciTime, "s")
for val, vec in zip(arpackSlv.val, arpackSlv.vec):
print("eigen value:", val)
print("eigen vector:")
for v in range(n):
print(vec[v])
EXAMPLES/PYARPACK/pyarpackDenseQRRR.py.in
+117
View File
@@ -0,0 +1,117 @@
#!/usr/bin/env python
from __future__ import print_function
import numpy as np
from pyarpack import denseQRRR as pyarpackSlv
# Build laplacian.
n = 4
Aij = np.array([], dtype='complex128')
for k in range(n):
for l in range(n):
if l == k:
Aij = np.append(Aij, np.complex128(np.complex( 200., 200.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
elif l == k-1:
Aij = np.append(Aij, np.complex128(np.complex(-101., -101.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
elif l == k+1:
Aij = np.append(Aij, np.complex128(np.complex( -99., -99.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
else:
Aij = np.append(Aij, np.complex128(np.complex( 0., 0.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
for idx, val in enumerate(Aij):
print("A[", idx, "] =", val)
A = (Aij, False) # raw format: Aij values, row ordered (or not).
# Get and tune arpack solver.
arpackSlv = pyarpackSlv.complexDouble() # Caution: complexDouble <=> np.array(..., dtype='complex128')
arpackSlv.verbose = 3 # Set to 0 to get a quiet solve.
arpackSlv.debug = 1 # Set to 0 to get a quiet solve.
arpackSlv.nbEV = 1
arpackSlv.nbCV = 2*arpackSlv.nbEV + 1
arpackSlv.mag = 'LM'
arpackSlv.maxIt = 200
arpackSlv.slvPvtThd = 1.e-6
arpackSlv.symPb = False
# Solve eigen problem.
rc = arpackSlv.solve(A)
assert rc == 0, "bad solve"
rc = arpackSlv.checkEigVec(A)
assert rc == 0, "bad checkEigVec"
# Print out results (mode selected, eigen vectors, eigen values, ...).
assert arpackSlv.nbEV == len(arpackSlv.val), "bad result"
print("\nresults:\n")
print("mode selected:", arpackSlv.mode)
print("nb iterations:", arpackSlv.nbIt)
print("Reverse Communication Interface time:", arpackSlv.rciTime, "s")
for val, vec in zip(arpackSlv.val, arpackSlv.vec):
print("eigen value:", val)
print("eigen vector:")
print(vec)
#######################################################################################
print("\n##########################################################################\n")
#######################################################################################
# Build laplacian.
n = 8
Aij = np.array([], dtype='complex64')
Bij = np.array([], dtype='complex64')
for k in range(n):
for l in range(n):
if l == k:
Aij = np.append(Aij, np.complex64(np.complex( 200., 200.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.complex64(np.complex( 33.3, 33.3))) # Casting value on append is MANDATORY or C++ won't get the expected type.
elif l == k-1:
Aij = np.append(Aij, np.complex64(np.complex(-101., -101.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.complex64(np.complex( 16.6, 16.6))) # Casting value on append is MANDATORY or C++ won't get the expected type.
elif l == k+1:
Aij = np.append(Aij, np.complex64(np.complex( -99., -99.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.complex64(np.complex( 16.6, 16.6))) # Casting value on append is MANDATORY or C++ won't get the expected type.
else:
Aij = np.append(Aij, np.complex64(np.complex( 0., 0.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.complex64(np.complex( 0., 0.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
for idx, val in enumerate(Aij):
print("A[", idx, "] =", val)
for idx, val in enumerate(Bij):
print("B[", idx, "] =", val)
A = (Aij, False) # raw format: Aij values, row ordered (or not).
B = (Bij, True) # raw format: Bij values, row ordered (or not).
# Get and tune arpack solver.
arpackSlv = pyarpackSlv.complexFloat() # Caution: complexFloat <=> np.array(..., dtype='complex64')
arpackSlv.verbose = 3 # Set to 0 to get a quiet solve.
arpackSlv.debug = 1 # Set to 0 to get a quiet solve.
arpackSlv.nbEV = 2
arpackSlv.nbCV = 2*arpackSlv.nbEV + 1
arpackSlv.mag = 'LM'
arpackSlv.maxIt = 200
arpackSlv.slvPvtThd = 1.e-6
arpackSlv.sigmaReal = 1
arpackSlv.symPb = False
# Solve eigen problem.
rc = arpackSlv.solve(A, B)
assert rc == 0, "bad solve"
rc = arpackSlv.checkEigVec(A, B, 1.e-2)
assert rc == 0, "bad checkEigVec"
# Print out results (mode selected, eigen vectors, eigen values, ...).
assert arpackSlv.nbEV == len(arpackSlv.val), "bad result"
print("\nresults:\n")
print("mode selected:", arpackSlv.mode)
print("nb iterations:", arpackSlv.nbIt)
print("Reverse Communication Interface time:", arpackSlv.rciTime, "s")
for val, vec in zip(arpackSlv.val, arpackSlv.vec):
print("eigen value:", val)
print("eigen vector:")
for v in range(n):
print(vec[v])
EXAMPLES/PYARPACK/pyarpackDrtSolver.hpp
+138
View File
@@ -0,0 +1,138 @@
#ifndef __PYARPACKDRTSOLVER_HPP__
#define __PYARPACKDRTSOLVER_HPP__
#include <string>
#include <arpackSolver.hpp>
#include <pyarpackServices.hpp>
#include "debug_c.hpp"
#include "stat_c.hpp"
#include <boost/python.hpp>
#include <boost/python/numpy.hpp>
namespace bp = boost::python;
namespace bn = boost::python::numpy;
template<typename RC, typename FD, typename EM, typename SLV>
class pyarpackSparseDrtSolver: public arpackDrtSolver<RC, FD, EM, SLV> {
// Public methods.
public:
pyarpackSparseDrtSolver(): arpackDrtSolver<RC, FD, EM, SLV>() {
debug = 0;
nopx = 0, nbx = 0, nrorth = 0, nitref = 0, nrstrt = 0;
tsaupd = 0., tsaup2 = 0., tsaitr = 0., tseigt = 0., tsgets = 0., tsapps = 0., tsconv = 0.;
tnaupd = 0., tnaup2 = 0., tnaitr = 0., tneigt = 0., tngets = 0., tnapps = 0., tnconv = 0.;
tcaupd = 0., tcaup2 = 0., tcaitr = 0., tceigt = 0., tcgets = 0., tcapps = 0., tcconv = 0.;
tmvopx = 0., tmvbx = 0., tgetv0 = 0., titref = 0., trvec = 0.;
};
int solve(bp::tuple & A, bp::tuple B = bp::tuple()) {
ARPACKSOLVERDEBUGSTAT();
EM M;
int rc = pyarpackServices<RC, EM>::buildSparseMatrice(A, M, debug, "A");
if (rc != 0) {pyarpackThrowError("build matrice from A KO"); return rc;}
bool stdPb = (bp::len(B) > 0) ? false : true;
EM N;
if (!stdPb) {
rc = pyarpackServices<RC, EM>::buildSparseMatrice(B, N, debug, "B");
if (rc != 0) {pyarpackThrowError("build matrice from B KO"); return rc;}
}
return arpackDrtSolver<RC, FD, EM, SLV>::solve(M, (stdPb ? NULL : &N));
};
int checkEigVec(bp::tuple const & A, bp::tuple const B = bp::tuple(), double const diffTol = 1.e-3) {
ARPACKSOLVERDEBUGSTAT();
EM M;
int rc = pyarpackServices<RC, EM>::buildSparseMatrice(A, M, debug, "A");
if (rc != 0) {pyarpackThrowError("build matrice from A KO"); return rc;}
bool stdPb = (bp::len(B) > 0) ? false : true;
EM N;
if (!stdPb) {
rc = pyarpackServices<RC, EM>::buildSparseMatrice(B, N, debug, "B");
if (rc != 0) {pyarpackThrowError("build matrice from B KO"); return rc;}
}
return arpackDrtSolver<RC, FD, EM, SLV>::checkEigVec(M, (stdPb ? NULL : &N), &diffTol);
};
// Public members.
public:
a_int debug;
a_int nopx, nbx, nrorth, nitref, nrstrt;
float tsaupd, tsaup2, tsaitr, tseigt, tsgets, tsapps, tsconv;
float tnaupd, tnaup2, tnaitr, tneigt, tngets, tnapps, tnconv;
float tcaupd, tcaup2, tcaitr, tceigt, tcgets, tcapps, tcconv;
float tmvopx, tmvbx, tgetv0, titref, trvec;
};
template<typename RC, typename FD, typename EM, typename SLV>
class pyarpackDenseDrtSolver: public arpackDrtSolver<RC, FD, EM, SLV> {
// Public methods.
public:
pyarpackDenseDrtSolver(): arpackDrtSolver<RC, FD, EM, SLV>() {
debug = 0;
nopx = 0, nbx = 0, nrorth = 0, nitref = 0, nrstrt = 0;
tsaupd = 0., tsaup2 = 0., tsaitr = 0., tseigt = 0., tsgets = 0., tsapps = 0., tsconv = 0.;
tnaupd = 0., tnaup2 = 0., tnaitr = 0., tneigt = 0., tngets = 0., tnapps = 0., tnconv = 0.;
tcaupd = 0., tcaup2 = 0., tcaitr = 0., tceigt = 0., tcgets = 0., tcapps = 0., tcconv = 0.;
tmvopx = 0., tmvbx = 0., tgetv0 = 0., titref = 0., trvec = 0.;
};
int solve(bp::tuple & A, bp::tuple B = bp::tuple()) {
ARPACKSOLVERDEBUGSTAT();
EM M;
int rc = pyarpackServices<RC, EM>::buildDenseMatrice(A, M, debug, "A");
if (rc != 0) {pyarpackThrowError("build matrice from A KO"); return rc;}
bool stdPb = (bp::len(B) > 0) ? false : true;
EM N;
if (!stdPb) {
rc = pyarpackServices<RC, EM>::buildDenseMatrice(B, N, debug, "B");
if (rc != 0) {pyarpackThrowError("build matrice from B KO"); return rc;}
}
return arpackDrtSolver<RC, FD, EM, SLV>::solve(M, (stdPb ? NULL : &N));
};
int checkEigVec(bp::tuple const & A, bp::tuple const B = bp::tuple(), double const diffTol = 1.e-3) {
ARPACKSOLVERDEBUGSTAT();
EM M;
int rc = pyarpackServices<RC, EM>::buildDenseMatrice(A, M, debug, "A");
if (rc != 0) {pyarpackThrowError("build matrice from A KO"); return rc;}
bool stdPb = (bp::len(B) > 0) ? false : true;
EM N;
if (!stdPb) {
rc = pyarpackServices<RC, EM>::buildDenseMatrice(B, N, debug, "B");
if (rc != 0) {pyarpackThrowError("build matrice from B KO"); return rc;}
}
return arpackDrtSolver<RC, FD, EM, SLV>::checkEigVec(M, (stdPb ? NULL : &N), &diffTol);
};
// Public members.
public:
a_int debug;
a_int nopx, nbx, nrorth, nitref, nrstrt;
float tsaupd, tsaup2, tsaitr, tseigt, tsgets, tsapps, tsconv;
float tnaupd, tnaup2, tnaitr, tneigt, tngets, tnapps, tnconv;
float tcaupd, tcaup2, tcaitr, tceigt, tcgets, tcapps, tcconv;
float tmvopx, tmvbx, tgetv0, titref, trvec;
};
#endif
// Local Variables:
// mode: c++
// c-file-style:"stroustrup"
// show-trailing-whitespace: t
// End:
/* vim: set sw=2 ts=2 et smartindent :*/
EXAMPLES/PYARPACK/pyarpackItrSolver.hpp
+80
View File
@@ -0,0 +1,80 @@
#ifndef __PYARPACKITRSOLVER_HPP__
#define __PYARPACKITRSOLVER_HPP__
#include <string>
#include <arpackSolver.hpp>
#include <pyarpackServices.hpp>
#include "debug_c.hpp"
#include "stat_c.hpp"
#include <boost/python.hpp>
#include <boost/python/numpy.hpp>
namespace bp = boost::python;
template<typename RC, typename FD, typename EM, typename SLV>
class pyarpackSparseItrSolver: public arpackItrSolver<RC, FD, EM, SLV> {
// Public methods.
public:
pyarpackSparseItrSolver(): arpackItrSolver<RC, FD, EM, SLV>() {
debug = 0;
nopx = 0, nbx = 0, nrorth = 0, nitref = 0, nrstrt = 0;
tsaupd = 0., tsaup2 = 0., tsaitr = 0., tseigt = 0., tsgets = 0., tsapps = 0., tsconv = 0.;
tnaupd = 0., tnaup2 = 0., tnaitr = 0., tneigt = 0., tngets = 0., tnapps = 0., tnconv = 0.;
tcaupd = 0., tcaup2 = 0., tcaitr = 0., tceigt = 0., tcgets = 0., tcapps = 0., tcconv = 0.;
tmvopx = 0., tmvbx = 0., tgetv0 = 0., titref = 0., trvec = 0.;
};
int solve(bp::tuple & A, bp::tuple B = bp::tuple()) {
ARPACKSOLVERDEBUGSTAT();
EM M;
int rc = pyarpackServices<RC, EM>::buildSparseMatrice(A, M, debug, "A");
if (rc != 0) {pyarpackThrowError("build matrice from A KO"); return rc;}
bool stdPb = (bp::len(B) > 0) ? false : true;
EM N;
if (!stdPb) {
rc = pyarpackServices<RC, EM>::buildSparseMatrice(B, N, debug, "B");
if (rc != 0) {pyarpackThrowError("build matrice from B KO"); return rc;}
}
return arpackItrSolver<RC, FD, EM, SLV>::solve(M, (stdPb ? NULL : &N));
};
int checkEigVec(bp::tuple const & A, bp::tuple const B = bp::tuple(), double const diffTol = 1.e-3) {
ARPACKSOLVERDEBUGSTAT();
EM M;
int rc = pyarpackServices<RC, EM>::buildSparseMatrice(A, M, debug, "A");
if (rc != 0) {pyarpackThrowError("build matrice from A KO"); return rc;}
bool stdPb = (bp::len(B) > 0) ? false : true;
EM N;
if (!stdPb) {
rc = pyarpackServices<RC, EM>::buildSparseMatrice(B, N, debug, "B");
if (rc != 0) {pyarpackThrowError("build matrice from B KO"); return rc;}
}
return arpackItrSolver<RC, FD, EM, SLV>::checkEigVec(M, (stdPb ? NULL : &N), &diffTol);
};
// Public members.
public:
a_int debug;
a_int nopx, nbx, nrorth, nitref, nrstrt;
float tsaupd, tsaup2, tsaitr, tseigt, tsgets, tsapps, tsconv;
float tnaupd, tnaup2, tnaitr, tneigt, tngets, tnapps, tnconv;
float tcaupd, tcaup2, tcaitr, tceigt, tcgets, tcapps, tcconv;
float tmvopx, tmvbx, tgetv0, titref, trvec;
};
#endif
// Local Variables:
// mode: c++
// c-file-style:"stroustrup"
// show-trailing-whitespace: t
// End:
/* vim: set sw=2 ts=2 et smartindent :*/
EXAMPLES/PYARPACK/pyarpackRestart.py.in
+120
View File
@@ -0,0 +1,120 @@
#!/usr/bin/env python
from __future__ import print_function
import numpy as np
from pyarpack import sparseBiCGDiag as pyarpackSlv
# Build laplacian.
n = 8
i = np.array([], dtype='@PYINT@')
j = np.array([], dtype='@PYINT@')
Aij = np.array([], dtype='float64')
for k in range(n):
for l in [k-1, k, k+1]:
if l < 0 or l > n-1:
continue
i = np.append(i, np.@PYINT@(k)) # Casting value on append is MANDATORY or C++ won't get the expected type.
j = np.append(j, np.@PYINT@(l)) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k:
Aij = np.append(Aij, np.float64( 200.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k-1 or l == k+1:
Aij = np.append(Aij, np.float64(-100.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
for k, l, Akl in zip(i, j, Aij):
print("A[", k, ",", l, "] =", Akl)
A = (n, i, j, Aij) # coo format: dimension, i 0-based indices, j 0-based indices, Aij values.
# Get and tune arpack solver.
arpackSlv = pyarpackSlv.double() # Caution: double <=> np.array(..., dtype='float64')
arpackSlv.verbose = 3 # Set to 0 to get a quiet solve.
arpackSlv.debug = 1 # Set to 0 to get a quiet solve.
arpackSlv.nbEV = 1
arpackSlv.nbCV = 2*arpackSlv.nbEV + 1
arpackSlv.mag = 'LM'
arpackSlv.maxIt = 200
arpackSlv.slvTol = 1.e-6
arpackSlv.slvMaxIt = 100
arpackSlv.dumpToFile = True # Dump eigen vectors to arpackSolver.*.out files.
arpackSlv.schur = True # Schur vectors and eigenvectors of A are the same if A is a normal matrix.
# Solve eigen problem.
rc = arpackSlv.solve(A)
assert rc == 0, "bad solve"
rc = arpackSlv.checkEigVec(A)
assert rc == 0, "bad checkEigVec"
nbIt1 = arpackSlv.nbIt
# Print out results (mode selected, eigen vectors, eigen values, ...).
assert arpackSlv.nbEV == len(arpackSlv.val), "bad result"
print("\nresults:\n")
print("mode selected:", arpackSlv.mode)
print("nb iterations:", arpackSlv.nbIt)
print("Reverse Communication Interface time:", arpackSlv.rciTime, "s")
for val, vec in zip(arpackSlv.val, arpackSlv.vec):
print("eigen value:", val)
print("eigen vector:")
print(vec)
#######################################################################################
print("\n##########################################################################\n")
#######################################################################################
# Build laplacian (similar-but-different from the previous one).
n = 8
i = np.array([], dtype='@PYINT@')
j = np.array([], dtype='@PYINT@')
Aij = np.array([], dtype='float64')
for k in range(n):
for l in [k-1, k, k+1]:
if l < 0 or l > n-1:
continue
i = np.append(i, np.@PYINT@(k)) # Casting value on append is MANDATORY or C++ won't get the expected type.
j = np.append(j, np.@PYINT@(l)) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k:
Aij = np.append(Aij, np.float64( 210.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k-1 or l == k+1:
Aij = np.append(Aij, np.float64( -90.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
for k, l, Akl in zip(i, j, Aij):
print("A[", k, ",", l, "] =", Akl)
A = (n, i, j, Aij) # coo format: dimension, i 0-based indices, j 0-based indices, Aij values.
# Get and tune arpack solver.
arpackSlv = pyarpackSlv.double() # Caution: double <=> np.array(..., dtype='float64')
arpackSlv.verbose = 3 # Set to 0 to get a quiet solve.
arpackSlv.debug = 1 # Set to 0 to get a quiet solve.
arpackSlv.nbEV = 1
arpackSlv.nbCV = 2*arpackSlv.nbEV + 1
arpackSlv.mag = 'LM'
arpackSlv.maxIt = 200
arpackSlv.slvTol = 1.e-6
arpackSlv.slvMaxIt = 100
arpackSlv.restartFromFile = True # Restart from eigen vectors found in arpackSolver.*.out files.
arpackSlv.schur = True # Schur vectors and eigenvectors of A are the same if A is a normal matrix.
# Solve eigen problem.
rc = arpackSlv.solve(A)
assert rc == 0, "bad solve"
rc = arpackSlv.checkEigVec(A)
assert rc == 0, "bad checkEigVec"
nbIt2 = arpackSlv.nbIt
assert nbIt2 < nbIt1, "bad restart" # Restart from the first solve to run the second solve for a similar-but-different A.
# Print out results (mode selected, eigen vectors, eigen values, ...).
assert arpackSlv.nbEV == len(arpackSlv.val), "bad result"
print("\nresults:\n")
print("mode selected:", arpackSlv.mode)
print("nb iterations:", arpackSlv.nbIt)
print("Reverse Communication Interface time:", arpackSlv.rciTime, "s")
for val, vec in zip(arpackSlv.val, arpackSlv.vec):
print("eigen value:", val)
print("eigen vector:")
print(vec)
EXAMPLES/PYARPACK/pyarpackServices.hpp
+194
View File
@@ -0,0 +1,194 @@
#ifndef __PYARPACKSERVICES_HPP__
#define __PYARPACKSERVICES_HPP__
#include <vector>
#include <string>
#include <complex>
#include <iostream>
#include <cmath> // sqrt.
#include <Eigen/Sparse>
#include <boost/python.hpp>
#include <boost/python/numpy.hpp>
namespace bp = boost::python;
namespace bn = boost::python::numpy;
#define ARPACKSOLVERMEMBER(pyarpackSolver) \
.def_readwrite("symPb", &pyarpackSolver<RC, FD, EM, SLV>::symPb, \
"symmetric problem - default: true") \
.def_readwrite("nbEV", &pyarpackSolver<RC, FD, EM, SLV>::nbEV, \
"number of eigen vectors to find - default: 1") \
.def_readwrite("nbCV", &pyarpackSolver<RC, FD, EM, SLV>::nbCV, \
"number of dimensions of the workspace - default: 3") \
.def_readwrite("tol", &pyarpackSolver<RC, FD, EM, SLV>::tol, \
"tolerance - default: 1.e-6") \
.def_readwrite("sigmaReal", &pyarpackSolver<RC, FD, EM, SLV>::sigmaReal, \
"shift over real axis - default: 0.") \
.def_readwrite("sigmaImag", &pyarpackSolver<RC, FD, EM, SLV>::sigmaImag, \
"shift over imaginary axis - default: 0.") \
.def_readwrite("dumpToFile", &pyarpackSolver<RC, FD, EM, SLV>::dumpToFile, \
"dump eigen vectors to arpackSolver.*.out files - default: false") \
.def_readwrite("restartFromFile", &pyarpackSolver<RC, FD, EM, SLV>::restartFromFile, \
"restart from eigen vectors found in arpackSolver.*.out files - default: false") \
.def_readwrite("mag", &pyarpackSolver<RC, FD, EM, SLV>::mag, \
"magnitude - default: LM") \
.def_readwrite("maxIt", &pyarpackSolver<RC, FD, EM, SLV>::maxIt, \
"maximum number of arpack iterations - default: 100") \
.def_readwrite("schur", &pyarpackSolver<RC, FD, EM, SLV>::schur, \
"compute schur vectors - default: false") \
.def_readwrite("verbose", &pyarpackSolver<RC, FD, EM, SLV>::verbose, \
"verbosity level - default: 0") \
.def_readonly ("stdPb", &pyarpackSolver<RC, FD, EM, SLV>::stdPb, \
"standard or generalised problem - default: true") \
.def_readonly ("val", &pyarpackSolver<RC, FD, EM, SLV>::val, \
"eigen values found") \
.def_readonly ("vec", &pyarpackSolver<RC, FD, EM, SLV>::vec, \
"eigen vectors found") \
.def_readonly ("mode", &pyarpackSolver<RC, FD, EM, SLV>::mode, \
"selected arpack mode (according to input options: std/gen, shift, ...)") \
.def_readonly ("nbIt", &pyarpackSolver<RC, FD, EM, SLV>::nbIt, \
"number of arpack iterations") \
.def_readonly ("imsTime", &pyarpackSolver<RC, FD, EM, SLV>::imsTime, \
"time spent to initialize the mode solver if needed") \
.def_readonly ("rciTime", &pyarpackSolver<RC, FD, EM, SLV>::rciTime, \
"time spent in Reverse Communication Interface") \
.def_readwrite("debug", &pyarpackSolver<RC, FD, EM, SLV>::debug, \
"debug traces (up to 3) - default: 0") \
#define ARPACKSOLVERDEBUGSTAT() \
if (debug > 3) debug = 3; \
debug_c(6, -6, debug, debug, debug, debug, debug, debug, debug, debug, debug, debug, debug, \
debug, debug, debug, debug, debug, debug, debug, debug, debug, debug, debug); \
stat_c(nopx, nbx, nrorth, nitref, nrstrt, tsaupd, tsaup2, \
tsaitr, tseigt, tsgets, tsapps, tsconv, tnaupd, tnaup2, \
tnaitr, tneigt, tngets, tnapps, tnconv, tcaupd, tcaup2, \
tcaitr, tceigt, tcgets, tcapps, tcconv, tmvopx, tmvbx, \
tgetv0, titref, trvec); \
void pyarpackThrowError(std::string const & msg) {
std::string const info = "Error: " + msg;
std::cerr << info << std::endl;
PyErr_SetString(PyExc_IndexError, info.c_str());
bp::throw_error_already_set();
};
template<typename RC, typename EM>
class pyarpackServices {
// Public methods.
public:
static int buildSparseMatrice(bp::tuple const & T, Eigen::SparseMatrix<RC> & M,
a_int const & debug, std::string const & msg) {
// Get boost data as C++ data.
if (bp::len(T) != 4) {pyarpackThrowError(msg + " must be a 3-tuple"); return 1;}
bp::extract<int> nExt(T[0]);
bp::extract<bn::ndarray> iExt(T[1]);
bp::extract<bn::ndarray> jExt(T[2]);
bp::extract<bn::ndarray> mijExt(T[3]);
if (! nExt.check()) {pyarpackThrowError(msg + "[0] must be an integer" ); return 1;}
if (! iExt.check()) {pyarpackThrowError(msg + "[1] must be numpy.array"); return 1;}
if (! jExt.check()) {pyarpackThrowError(msg + "[2] must be numpy.array"); return 1;}
if (!mijExt.check()) {pyarpackThrowError(msg + "[3] must be numpy.array"); return 1;}
bn::ndarray iArray = iExt();
bn::ndarray jArray = jExt();
bn::ndarray mijArray = mijExt();
if (iArray.get_dtype() != bn::dtype::get_builtin<a_int>()) {pyarpackThrowError(msg + "[1] type is not consistent"); return 1;}
if (jArray.get_dtype() != bn::dtype::get_builtin<a_int>()) {pyarpackThrowError(msg + "[2] type is not consistent"); return 1;}
if (mijArray.get_dtype() != bn::dtype::get_builtin<RC>() ) {pyarpackThrowError(msg + "[3] type is not consistent with arpack type"); return 1;}
a_int iSz = iArray.shape(0);
a_int * iPtr = reinterpret_cast<a_int*>(iArray.get_data());
a_int jSz = jArray.shape(0);
a_int * jPtr = reinterpret_cast<a_int*>(jArray.get_data());
a_int mSz = mijArray.shape(0);
RC * mPtr = reinterpret_cast<RC*>(mijArray.get_data());
if (iSz != jSz) {pyarpackThrowError(msg + "[1] and " + msg + "[2] must have same lenght"); return 1;}
if (iSz != mSz) {pyarpackThrowError(msg + "[1] and " + msg + "[3] must have same lenght"); return 1;}
// Debug on demand: casting value on numpy.append is MANDATORY or C++ won't get the expected type..
for (auto k = 0; debug && k < mSz; k++) {
std::cout << "pyarpackServices::buildSparseMatrice - " << msg << "[" << iPtr[k] << ", " << jPtr[k] << "] = " << mPtr[k] << std::endl;
};
// Build sparse matrice.
a_uint n = nExt();
a_uint iMin = n+1, jMin = n+1;
for (auto k = 0; k < mSz; k++) {
if (iPtr[k] < iMin) iMin = iPtr[k];
if (jPtr[k] < jMin) jMin = jPtr[k];
};
if (iMin != 0 && iMin != 1) {pyarpackThrowError(msg + ": smallest row indice must be 0 or 1"); return 1;}
if (jMin != 0 && jMin != 1) {pyarpackThrowError(msg + ": smallest column indice must be 0 or 1"); return 1;}
a_int iBased = 0, jBased = 0;
if (iMin == 1) iBased = 1;
if (jMin == 1) jBased = 1;
M = Eigen::SparseMatrix<RC>(n, n); // Set matrice dimensions.
std::vector<Eigen::Triplet<RC>> triplets;
a_uint nnz = mSz;
triplets.reserve(nnz);
for (auto k = 0; k < nnz; k++) triplets.emplace_back(iPtr[k] - iBased, jPtr[k] - jBased, mPtr[k]);
M.setFromTriplets(triplets.begin(), triplets.end()); // Set all (i, j, Mij).
return 0;
};
static int buildDenseMatrice(bp::tuple const & T, Eigen::Matrix<RC, Eigen::Dynamic, Eigen::Dynamic> & M,
a_int const & debug, std::string const & msg) {
// Get boost data as C++ data.
if (bp::len(T) != 2) {pyarpackThrowError(msg + " must be a 2-tuple"); return 1;}
bp::extract<bn::ndarray> mijExt(T[0]);
bp::extract<bool> oExt(T[1]);
if (!mijExt.check()) {pyarpackThrowError(msg + " must be numpy.array"); return 1;}
if ( !oExt.check()) {pyarpackThrowError(msg + " must be a boolean"); return 1;}
bn::ndarray mijArray = mijExt();
bool rowOrdered = oExt();
if (mijArray.get_dtype() != bn::dtype::get_builtin<RC>()) {pyarpackThrowError(msg + " type is not consistent with arpack type"); return 1;}
a_int mSz = mijArray.shape(0);
RC * mPtr = reinterpret_cast<RC*>(mijArray.get_data());
a_uint n = std::sqrt(mSz);
if (n*n != mSz) {pyarpackThrowError(msg + " must be a squared matrice"); return 1;}
// Debug on demand: casting value on numpy.append is MANDATORY or C++ won't get the expected type..
for (auto k = 0; debug && k < mSz; k++) {
std::cout << "pyarpackServices::buildDenseMatrice - " << msg << "[" << k << "] = " << mPtr[k] << std::endl;
};
// Build dense matrice.
M = Eigen::Matrix<RC, Eigen::Dynamic, Eigen::Dynamic>(n, n); // Set matrice dimensions.
M.setZero(n, n); // Avoid spurious/random values which may break solves (LU, QR, ...).
if (rowOrdered) {
for (size_t k = 0; k < n; k++) {
for (size_t l = 0; l < n; l++) M(k, l) = mPtr[l+k*n];
}
}
else {
for (size_t l = 0; l < n; l++) {
for (size_t k = 0; k < n; k++) M(k, l) = mPtr[k+l*n];
}
}
return 0;
};
};
#endif
// Local Variables:
// mode: c++
// c-file-style:"stroustrup"
// show-trailing-whitespace: t
// End:
/* vim: set sw=2 ts=2 et smartindent :*/
EXAMPLES/PYARPACK/pyarpackSparseBiCGDiag.py.in
+119
View File
@@ -0,0 +1,119 @@
#!/usr/bin/env python
from __future__ import print_function
import numpy as np
from pyarpack import sparseBiCGDiag as pyarpackSlv
# Build laplacian.
n = 4
i = np.array([], dtype='@PYINT@')
j = np.array([], dtype='@PYINT@')
Aij = np.array([], dtype='float64')
for k in range(n):
for l in [k-1, k, k+1]:
if l < 0 or l > n-1:
continue
i = np.append(i, np.@PYINT@(k)) # Casting value on append is MANDATORY or C++ won't get the expected type.
j = np.append(j, np.@PYINT@(l)) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k:
Aij = np.append(Aij, np.float64( 200.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k-1 or l == k+1:
Aij = np.append(Aij, np.float64(-100.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
for k, l, Akl in zip(i, j, Aij):
print("A[", k, ",", l, "] =", Akl)
A = (n, i, j, Aij) # coo format: dimension, i 0-based indices, j 0-based indices, Aij values.
# Get and tune arpack solver.
arpackSlv = pyarpackSlv.double() # Caution: double <=> np.array(..., dtype='float64')
arpackSlv.verbose = 3 # Set to 0 to get a quiet solve.
arpackSlv.debug = 1 # Set to 0 to get a quiet solve.
arpackSlv.nbEV = 1
arpackSlv.nbCV = 2*arpackSlv.nbEV + 1
arpackSlv.mag = 'LM'
arpackSlv.maxIt = 200
arpackSlv.slvTol = 1.e-6
arpackSlv.slvMaxIt = 100
# Solve eigen problem.
rc = arpackSlv.solve(A)
assert rc == 0, "bad solve"
rc = arpackSlv.checkEigVec(A)
assert rc == 0, "bad checkEigVec"
# Print out results (mode selected, eigen vectors, eigen values, ...).
assert arpackSlv.nbEV == len(arpackSlv.val), "bad result"
print("\nresults:\n")
print("mode selected:", arpackSlv.mode)
print("nb iterations:", arpackSlv.nbIt)
print("Reverse Communication Interface time:", arpackSlv.rciTime, "s")
for val, vec in zip(arpackSlv.val, arpackSlv.vec):
print("eigen value:", val)
print("eigen vector:")
print(vec)
#######################################################################################
print("\n##########################################################################\n")
#######################################################################################
# Build laplacian.
n = 8
i = np.array([], dtype='@PYINT@')
j = np.array([], dtype='@PYINT@')
Aij = np.array([], dtype='float32')
Bij = np.array([], dtype='float32')
for k in range(n):
for l in [k-1, k, k+1]:
if l < 0 or l > n-1:
continue
i = np.append(i, np.@PYINT@(k+1)) # Casting value on append is MANDATORY or C++ won't get the expected type.
j = np.append(j, np.@PYINT@(l+1)) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k:
Aij = np.append(Aij, np.float32( 200.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.float32( 33.3)) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k-1 or l == k+1:
Aij = np.append(Aij, np.float32(-100.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.float32( 16.6)) # Casting value on append is MANDATORY or C++ won't get the expected type.
for k, l, Akl in zip(i, j, Aij):
print("A[", k, ",", l, "] =", Akl)
for k, l, Bkl in zip(i, j, Bij):
print("B[", k, ",", l, "] =", Bkl)
A = (n, i, j, Aij) # coo format: dimension, i 1-based indices, j 1-based indices, Aij values.
B = (n, i, j, Bij) # coo format: dimension, i 1-based indices, j 1-based indices, Bij values.
# Get and tune arpack solver.
arpackSlv = pyarpackSlv.float() # Caution: float <=> np.array(..., dtype='float32')
arpackSlv.verbose = 3 # Set to 0 to get a quiet solve.
arpackSlv.debug = 1 # Set to 0 to get a quiet solve.
arpackSlv.nbEV = 2
arpackSlv.nbCV = 2*arpackSlv.nbEV + 1
arpackSlv.mag = 'LM'
arpackSlv.maxIt = 200
arpackSlv.slvTol = 1.e-6
arpackSlv.slvMaxIt = 100
arpackSlv.sigmaReal = 1
# Solve eigen problem.
rc = arpackSlv.solve(A, B)
assert rc == 0, "bad solve"
rc = arpackSlv.checkEigVec(A, B, 1.e-2)
assert rc == 0, "bad checkEigVec"
# Print out results (mode selected, eigen vectors, eigen values, ...).
assert arpackSlv.nbEV == len(arpackSlv.val), "bad result"
print("\nresults:\n")
print("mode selected:", arpackSlv.mode)
print("nb iterations:", arpackSlv.nbIt)
print("Reverse Communication Interface time:", arpackSlv.rciTime, "s")
for val, vec in zip(arpackSlv.val, arpackSlv.vec):
print("eigen value:", val)
print("eigen vector:")
for v in range(n):
print(vec[v])
EXAMPLES/PYARPACK/pyarpackSparseBiCGILU.py.in
+131
View File
@@ -0,0 +1,131 @@
#!/usr/bin/env python
from __future__ import print_function
import numpy as np
from pyarpack import sparseBiCGILU as pyarpackSlv
# Build laplacian.
n = 4
i = np.array([], dtype='@PYINT@')
j = np.array([], dtype='@PYINT@')
Aij = np.array([], dtype='complex128')
for k in range(n):
for l in [k-1, k, k+1]:
if l < 0 or l > n-1:
continue
i = np.append(i, np.@PYINT@(k)) # Casting value on append is MANDATORY or C++ won't get the expected type.
j = np.append(j, np.@PYINT@(l)) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k:
Aij = np.append(Aij, np.complex128(np.complex( 200., 200.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k-1:
Aij = np.append(Aij, np.complex128(np.complex(-101., -101.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k+1:
Aij = np.append(Aij, np.complex128(np.complex( -99., -99.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
for k, l, Akl in zip(i, j, Aij):
print("A[", k, ",", l, "] =", Akl)
A = (n, i, j, Aij) # coo format: dimension, i 0-based indices, j 0-based indices, Aij values.
# Get and tune arpack solver.
arpackSlv = pyarpackSlv.complexDouble() # Caution: complexDouble <=> np.array(..., dtype='complex128')
arpackSlv.verbose = 3 # Set to 0 to get a quiet solve.
arpackSlv.debug = 1 # Set to 0 to get a quiet solve.
arpackSlv.nbEV = 1
arpackSlv.nbCV = 2*arpackSlv.nbEV + 1
arpackSlv.mag = 'LM'
arpackSlv.maxIt = 200
arpackSlv.slvTol = 1.e-6
arpackSlv.slvMaxIt = 200
arpackSlv.slvILUDropTol = 1.
arpackSlv.slvILUFillFactor = 2
arpackSlv.symPb = False
# Solve eigen problem.
rc = arpackSlv.solve(A)
assert rc == 0, "bad solve"
rc = arpackSlv.checkEigVec(A)
assert rc == 0, "bad checkEigVec"
# Print out results (mode selected, eigen vectors, eigen values, ...).
assert arpackSlv.nbEV == len(arpackSlv.val), "bad result"
print("\nresults:\n")
print("mode selected:", arpackSlv.mode)
print("nb iterations:", arpackSlv.nbIt)
print("Reverse Communication Interface time:", arpackSlv.rciTime, "s")
for val, vec in zip(arpackSlv.val, arpackSlv.vec):
print("eigen value:", val)
print("eigen vector:")
print(vec)
#######################################################################################
print("\n##########################################################################\n")
#######################################################################################
# Build laplacian.
n = 8
i = np.array([], dtype='@PYINT@')
j = np.array([], dtype='@PYINT@')
Aij = np.array([], dtype='complex64')
Bij = np.array([], dtype='complex64')
for k in range(n):
for l in [k-1, k, k+1]:
if l < 0 or l > n-1:
continue
i = np.append(i, np.@PYINT@(k+1)) # Casting value on append is MANDATORY or C++ won't get the expected type.
j = np.append(j, np.@PYINT@(l+1)) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k:
Aij = np.append(Aij, np.complex64(np.complex( 200., 200.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.complex64(np.complex( 33.3, 33.3))) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k-1:
Aij = np.append(Aij, np.complex64(np.complex(-101., -101.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.complex64(np.complex( 16.6, 16.6))) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k+1:
Aij = np.append(Aij, np.complex64(np.complex( -99., -99.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.complex64(np.complex( 16.6, 16.6))) # Casting value on append is MANDATORY or C++ won't get the expected type.
for k, l, Akl in zip(i, j, Aij):
print("A[", k, ",", l, "] =", Akl)
for k, l, Bkl in zip(i, j, Bij):
print("B[", k, ",", l, "] =", Bkl)
A = (n, i, j, Aij) # coo format: dimension, i 1-based indices, j 1-based indices, Aij values.
B = (n, i, j, Bij) # coo format: dimension, i 1-based indices, j 1-based indices, Bij values.
# Get and tune arpack solver.
arpackSlv = pyarpackSlv.complexFloat() # Caution: complexFloat <=> np.array(..., dtype='complex64')
arpackSlv.verbose = 3 # Set to 0 to get a quiet solve.
arpackSlv.debug = 1 # Set to 0 to get a quiet solve.
arpackSlv.nbEV = 2
arpackSlv.nbCV = 2*arpackSlv.nbEV + 1
arpackSlv.mag = 'LM'
arpackSlv.maxIt = 200
arpackSlv.slvTol = 1.e-6
arpackSlv.slvMaxIt = 200
arpackSlv.slvILUDropTol = 1.
arpackSlv.slvILUFillFactor = 2
arpackSlv.sigmaReal = 1
arpackSlv.sigmaImag = 1
arpackSlv.symPb = False
# Solve eigen problem.
rc = arpackSlv.solve(A, B)
assert rc == 0, "bad solve"
rc = arpackSlv.checkEigVec(A, B, 1.e-2)
assert rc == 0, "bad checkEigVec"
# Print out results (mode selected, eigen vectors, eigen values, ...).
assert arpackSlv.nbEV == len(arpackSlv.val), "bad result"
print("\nresults:\n")
print("mode selected:", arpackSlv.mode)
print("nb iterations:", arpackSlv.nbIt)
print("Reverse Communication Interface time:", arpackSlv.rciTime, "s")
for val, vec in zip(arpackSlv.val, arpackSlv.vec):
print("eigen value:", val)
print("eigen vector:")
for v in range(n):
print(vec[v])
EXAMPLES/PYARPACK/pyarpackSparseCGDiag.py.in
+119
View File
@@ -0,0 +1,119 @@
#!/usr/bin/env python
from __future__ import print_function
import numpy as np
from pyarpack import sparseCGDiag as pyarpackSlv
# Build laplacian.
n = 4
i = np.array([], dtype='@PYINT@')
j = np.array([], dtype='@PYINT@')
Aij = np.array([], dtype='float64')
for k in range(n):
for l in [k-1, k, k+1]:
if l < 0 or l > n-1:
continue
i = np.append(i, np.@PYINT@(k)) # Casting value on append is MANDATORY or C++ won't get the expected type.
j = np.append(j, np.@PYINT@(l)) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k:
Aij = np.append(Aij, np.float64( 200.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k-1 or l == k+1:
Aij = np.append(Aij, np.float64(-100.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
for k, l, Akl in zip(i, j, Aij):
print("A[", k, ",", l, "] =", Akl)
A = (n, i, j, Aij) # coo format: dimension, i 0-based indices, j 0-based indices, Aij values.
# Get and tune arpack solver.
arpackSlv = pyarpackSlv.double() # Caution: double <=> np.array(..., dtype='float64')
arpackSlv.verbose = 3 # Set to 0 to get a quiet solve.
arpackSlv.debug = 1 # Set to 0 to get a quiet solve.
arpackSlv.nbEV = 1
arpackSlv.nbCV = 2*arpackSlv.nbEV + 1
arpackSlv.mag = 'LM'
arpackSlv.maxIt = 200
arpackSlv.slvTol = 1.e-6
arpackSlv.slvMaxIt = 100
# Solve eigen problem.
rc = arpackSlv.solve(A)
assert rc == 0, "bad solve"
rc = arpackSlv.checkEigVec(A)
assert rc == 0, "bad checkEigVec"
# Print out results (mode selected, eigen vectors, eigen values, ...).
assert arpackSlv.nbEV == len(arpackSlv.val), "bad result"
print("\nresults:\n")
print("mode selected:", arpackSlv.mode)
print("nb iterations:", arpackSlv.nbIt)
print("Reverse Communication Interface time:", arpackSlv.rciTime, "s")
for val, vec in zip(arpackSlv.val, arpackSlv.vec):
print("eigen value:", val)
print("eigen vector:")
print(vec)
#######################################################################################
print("\n##########################################################################\n")
#######################################################################################
# Build laplacian.
n = 8
i = np.array([], dtype='@PYINT@')
j = np.array([], dtype='@PYINT@')
Aij = np.array([], dtype='float32')
Bij = np.array([], dtype='float32')
for k in range(n):
for l in [k-1, k, k+1]:
if l < 0 or l > n-1:
continue
i = np.append(i, np.@PYINT@(k+1)) # Casting value on append is MANDATORY or C++ won't get the expected type.
j = np.append(j, np.@PYINT@(l+1)) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k:
Aij = np.append(Aij, np.float32( 200.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.float32( 33.3)) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k-1 or l == k+1:
Aij = np.append(Aij, np.float32(-100.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.float32( 16.6)) # Casting value on append is MANDATORY or C++ won't get the expected type.
for k, l, Akl in zip(i, j, Aij):
print("A[", k, ",", l, "] =", Akl)
for k, l, Bkl in zip(i, j, Bij):
print("B[", k, ",", l, "] =", Bkl)
A = (n, i, j, Aij) # coo format: dimension, i 1-based indices, j 1-based indices, Aij values.
B = (n, i, j, Bij) # coo format: dimension, i 1-based indices, j 1-based indices, Bij values.
# Get and tune arpack solver.
arpackSlv = pyarpackSlv.float() # Caution: float <=> np.array(..., dtype='float32')
arpackSlv.verbose = 3 # Set to 0 to get a quiet solve.
arpackSlv.debug = 1 # Set to 0 to get a quiet solve.
arpackSlv.nbEV = 2
arpackSlv.nbCV = 2*arpackSlv.nbEV + 1
arpackSlv.mag = 'LM'
arpackSlv.maxIt = 200
arpackSlv.slvTol = 1.e-6
arpackSlv.slvMaxIt = 100
arpackSlv.sigmaReal = 1
# Solve eigen problem.
rc = arpackSlv.solve(A, B)
assert rc == 0, "bad solve"
rc = arpackSlv.checkEigVec(A, B, 1.e-2)
assert rc == 0, "bad checkEigVec"
# Print out results (mode selected, eigen vectors, eigen values, ...).
assert arpackSlv.nbEV == len(arpackSlv.val), "bad result"
print("\nresults:\n")
print("mode selected:", arpackSlv.mode)
print("nb iterations:", arpackSlv.nbIt)
print("Reverse Communication Interface time:", arpackSlv.rciTime, "s")
for val, vec in zip(arpackSlv.val, arpackSlv.vec):
print("eigen value:", val)
print("eigen vector:")
for v in range(n):
print(vec[v])
EXAMPLES/PYARPACK/pyarpackSparseCGILU.py.in
+131
View File
@@ -0,0 +1,131 @@
#!/usr/bin/env python
from __future__ import print_function
import numpy as np
from pyarpack import sparseCGILU as pyarpackSlv
# Build laplacian.
n = 4
i = np.array([], dtype='@PYINT@')
j = np.array([], dtype='@PYINT@')
Aij = np.array([], dtype='complex128')
for k in range(n):
for l in [k-1, k, k+1]:
if l < 0 or l > n-1:
continue
i = np.append(i, np.@PYINT@(k)) # Casting value on append is MANDATORY or C++ won't get the expected type.
j = np.append(j, np.@PYINT@(l)) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k:
Aij = np.append(Aij, np.complex128(np.complex( 200., 200.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k-1:
Aij = np.append(Aij, np.complex128(np.complex(-101., -101.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k+1:
Aij = np.append(Aij, np.complex128(np.complex( -99., -99.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
for k, l, Akl in zip(i, j, Aij):
print("A[", k, ",", l, "] =", Akl)
A = (n, i, j, Aij) # coo format: dimension, i 0-based indices, j 0-based indices, Aij values.
# Get and tune arpack solver.
arpackSlv = pyarpackSlv.complexDouble() # Caution: complexDouble <=> np.array(..., dtype='complex128')
arpackSlv.verbose = 3 # Set to 0 to get a quiet solve.
arpackSlv.debug = 1 # Set to 0 to get a quiet solve.
arpackSlv.nbEV = 1
arpackSlv.nbCV = 2*arpackSlv.nbEV + 1
arpackSlv.mag = 'LM'
arpackSlv.maxIt = 200
arpackSlv.slvTol = 1.e-6
arpackSlv.slvMaxIt = 200
arpackSlv.slvILUDropTol = 1.
arpackSlv.slvILUFillFactor = 2
arpackSlv.symPb = False
# Solve eigen problem.
rc = arpackSlv.solve(A)
assert rc == 0, "bad solve"
rc = arpackSlv.checkEigVec(A)
assert rc == 0, "bad checkEigVec"
# Print out results (mode selected, eigen vectors, eigen values, ...).
assert arpackSlv.nbEV == len(arpackSlv.val), "bad result"
print("\nresults:\n")
print("mode selected:", arpackSlv.mode)
print("nb iterations:", arpackSlv.nbIt)
print("Reverse Communication Interface time:", arpackSlv.rciTime, "s")
for val, vec in zip(arpackSlv.val, arpackSlv.vec):
print("eigen value:", val)
print("eigen vector:")
print(vec)
#######################################################################################
print("\n##########################################################################\n")
#######################################################################################
# Build laplacian.
n = 8
i = np.array([], dtype='@PYINT@')
j = np.array([], dtype='@PYINT@')
Aij = np.array([], dtype='complex64')
Bij = np.array([], dtype='complex64')
for k in range(n):
for l in [k-1, k, k+1]:
if l < 0 or l > n-1:
continue
i = np.append(i, np.@PYINT@(k+1)) # Casting value on append is MANDATORY or C++ won't get the expected type.
j = np.append(j, np.@PYINT@(l+1)) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k:
Aij = np.append(Aij, np.complex64(np.complex( 200., 200.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.complex64(np.complex( 33.3, 33.3))) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k-1:
Aij = np.append(Aij, np.complex64(np.complex(-101., -101.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.complex64(np.complex( 16.6, 16.6))) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k+1:
Aij = np.append(Aij, np.complex64(np.complex( -99., -99.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.complex64(np.complex( 16.6, 16.6))) # Casting value on append is MANDATORY or C++ won't get the expected type.
for k, l, Akl in zip(i, j, Aij):
print("A[", k, ",", l, "] =", Akl)
for k, l, Bkl in zip(i, j, Bij):
print("B[", k, ",", l, "] =", Bkl)
A = (n, i, j, Aij) # coo format: dimension, i 1-based indices, j 1-based indices, Aij values.
B = (n, i, j, Bij) # coo format: dimension, i 1-based indices, j 1-based indices, Bij values.
# Get and tune arpack solver.
arpackSlv = pyarpackSlv.complexFloat() # Caution: complexFloat <=> np.array(..., dtype='complex64')
arpackSlv.verbose = 3 # Set to 0 to get a quiet solve.
arpackSlv.debug = 1 # Set to 0 to get a quiet solve.
arpackSlv.nbEV = 2
arpackSlv.nbCV = 2*arpackSlv.nbEV + 1
arpackSlv.mag = 'LM'
arpackSlv.maxIt = 200
arpackSlv.slvTol = 1.e-6
arpackSlv.slvMaxIt = 200
arpackSlv.slvILUDropTol = 1.
arpackSlv.slvILUFillFactor = 2
arpackSlv.sigmaReal = 1
arpackSlv.sigmaImag = 1
arpackSlv.symPb = False
# Solve eigen problem.
rc = arpackSlv.solve(A, B)
assert rc == 0, "bad solve"
rc = arpackSlv.checkEigVec(A, B, 1.)
assert rc == 0, "bad checkEigVec"
# Print out results (mode selected, eigen vectors, eigen values, ...).
assert arpackSlv.nbEV == len(arpackSlv.val), "bad result"
print("\nresults:\n")
print("mode selected:", arpackSlv.mode)
print("nb iterations:", arpackSlv.nbIt)
print("Reverse Communication Interface time:", arpackSlv.rciTime, "s")
for val, vec in zip(arpackSlv.val, arpackSlv.vec):
print("eigen value:", val)
print("eigen vector:")
for v in range(n):
print(vec[v])
EXAMPLES/PYARPACK/pyarpackSparseLDLT.py.in
+127
View File
@@ -0,0 +1,127 @@
#!/usr/bin/env python
from __future__ import print_function
import numpy as np
from pyarpack import sparseLDLT as pyarpackSlv
# Build laplacian.
n = 4
i = np.array([], dtype='@PYINT@')
j = np.array([], dtype='@PYINT@')
Aij = np.array([], dtype='complex128')
for k in range(n):
for l in [k-1, k, k+1]:
if l < 0 or l > n-1:
continue
i = np.append(i, np.@PYINT@(k)) # Casting value on append is MANDATORY or C++ won't get the expected type.
j = np.append(j, np.@PYINT@(l)) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k:
Aij = np.append(Aij, np.complex128(np.complex( 200., 200.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k-1:
Aij = np.append(Aij, np.complex128(np.complex(-101., -101.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k+1:
Aij = np.append(Aij, np.complex128(np.complex( -99., -99.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
for k, l, Akl in zip(i, j, Aij):
print("A[", k, ",", l, "] =", Akl)
A = (n, i, j, Aij) # coo format: dimension, i 0-based indices, j 0-based indices, Aij values.
# Get and tune arpack solver.
arpackSlv = pyarpackSlv.complexDouble() # Caution: complexDouble <=> np.array(..., dtype='complex128')
arpackSlv.verbose = 3 # Set to 0 to get a quiet solve.
arpackSlv.debug = 1 # Set to 0 to get a quiet solve.
arpackSlv.nbEV = 1
arpackSlv.nbCV = 2*arpackSlv.nbEV + 1
arpackSlv.mag = 'LM'
arpackSlv.maxIt = 200
arpackSlv.slvOffset = 0.
arpackSlv.slvScale = 1.
arpackSlv.symPb = False
# Solve eigen problem.
rc = arpackSlv.solve(A)
assert rc == 0, "bad solve"
rc = arpackSlv.checkEigVec(A)
assert rc == 0, "bad checkEigVec"
# Print out results (mode selected, eigen vectors, eigen values, ...).
assert arpackSlv.nbEV == len(arpackSlv.val), "bad result"
print("\nresults:\n")
print("mode selected:", arpackSlv.mode)
print("nb iterations:", arpackSlv.nbIt)
print("Reverse Communication Interface time:", arpackSlv.rciTime, "s")
for val, vec in zip(arpackSlv.val, arpackSlv.vec):
print("eigen value:", val)
print("eigen vector:")
print(vec)
#######################################################################################
print("\n##########################################################################\n")
#######################################################################################
# Build laplacian.
n = 8
i = np.array([], dtype='@PYINT@')
j = np.array([], dtype='@PYINT@')
Aij = np.array([], dtype='complex64')
Bij = np.array([], dtype='complex64')
for k in range(n):
for l in [k-1, k, k+1]:
if l < 0 or l > n-1:
continue
i = np.append(i, np.@PYINT@(k+1)) # Casting value on append is MANDATORY or C++ won't get the expected type.
j = np.append(j, np.@PYINT@(l+1)) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k:
Aij = np.append(Aij, np.complex64(np.complex( 200., 200.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.complex64(np.complex( 33.3, 33.3))) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k-1:
Aij = np.append(Aij, np.complex64(np.complex(-101., -101.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.complex64(np.complex( 16.6, 16.6))) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k+1:
Aij = np.append(Aij, np.complex64(np.complex( -99., -99.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.complex64(np.complex( 16.6, 16.6))) # Casting value on append is MANDATORY or C++ won't get the expected type.
for k, l, Akl in zip(i, j, Aij):
print("A[", k, ",", l, "] =", Akl)
for k, l, Bkl in zip(i, j, Bij):
print("B[", k, ",", l, "] =", Bkl)
A = (n, i, j, Aij) # coo format: dimension, i 1-based indices, j 1-based indices, Aij values.
B = (n, i, j, Bij) # coo format: dimension, i 1-based indices, j 1-based indices, Bij values.
# Get and tune arpack solver.
arpackSlv = pyarpackSlv.complexFloat() # Caution: complexFloat <=> np.array(..., dtype='complex64')
arpackSlv.verbose = 3 # Set to 0 to get a quiet solve.
arpackSlv.debug = 1 # Set to 0 to get a quiet solve.
arpackSlv.nbEV = 2
arpackSlv.nbCV = 2*arpackSlv.nbEV + 1
arpackSlv.mag = 'LM'
arpackSlv.maxIt = 200
arpackSlv.slvOffset = 0.
arpackSlv.slvScale = 1.
arpackSlv.sigmaReal = 1
arpackSlv.sigmaImag = 1
arpackSlv.symPb = False
# Solve eigen problem.
rc = arpackSlv.solve(A, B)
assert rc == 0, "bad solve"
rc = arpackSlv.checkEigVec(A, B, 1.e-2)
assert rc == 0, "bad checkEigVec"
# Print out results (mode selected, eigen vectors, eigen values, ...).
assert arpackSlv.nbEV == len(arpackSlv.val), "bad result"
print("\nresults:\n")
print("mode selected:", arpackSlv.mode)
print("nb iterations:", arpackSlv.nbIt)
print("Reverse Communication Interface time:", arpackSlv.rciTime, "s")
for val, vec in zip(arpackSlv.val, arpackSlv.vec):
print("eigen value:", val)
print("eigen vector:")
for v in range(n):
print(vec[v])
EXAMPLES/PYARPACK/pyarpackSparseLLT.py.in
+118
View File
@@ -0,0 +1,118 @@
#!/usr/bin/env python
from __future__ import print_function
import numpy as np
from pyarpack import sparseLLT as pyarpackSlv
# Build laplacian.
n = 4
i = np.array([], dtype='@PYINT@')
j = np.array([], dtype='@PYINT@')
Aij = np.array([], dtype='float64')
for k in range(n):
for l in [k-1, k, k+1]:
if l < 0 or l > n-1:
continue
i = np.append(i, np.@PYINT@(k)) # Casting value on append is MANDATORY or C++ won't get the expected type.
j = np.append(j, np.@PYINT@(l)) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k:
Aij = np.append(Aij, np.float64( 200.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k-1 or l == k+1:
Aij = np.append(Aij, np.float64(-100.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
for k, l, Akl in zip(i, j, Aij):
print("A[", k, ",", l, "] =", Akl)
A = (n, i, j, Aij) # coo format: dimension, i 0-based indices, j 0-based indices, Aij values.
# Get and tune arpack solver.
arpackSlv = pyarpackSlv.double() # Caution: double <=> np.array(..., dtype='float64')
arpackSlv.verbose = 3 # Set to 0 to get a quiet solve.
arpackSlv.debug = 1 # Set to 0 to get a quiet solve.
arpackSlv.nbEV = 1
arpackSlv.nbCV = 2*arpackSlv.nbEV + 1
arpackSlv.mag = 'LM'
arpackSlv.maxIt = 200
arpackSlv.slvOffset = 0.
arpackSlv.slvScale = 1.
# Solve eigen problem.
rc = arpackSlv.solve(A)
assert rc == 0, "bad solve"
rc = arpackSlv.checkEigVec(A)
assert rc == 0, "bad checkEigVec"
# Print out results (mode selected, eigen vectors, eigen values, ...).
assert arpackSlv.nbEV == len(arpackSlv.val), "bad result"
print("\nresults:\n")
print("mode selected:", arpackSlv.mode)
print("nb iterations:", arpackSlv.nbIt)
print("Reverse Communication Interface time:", arpackSlv.rciTime, "s")
for val, vec in zip(arpackSlv.val, arpackSlv.vec):
print("eigen value:", val)
print("eigen vector:")
print(vec)
#######################################################################################
print("\n##########################################################################\n")
#######################################################################################
# Build laplacian.
n = 8
i = np.array([], dtype='@PYINT@')
j = np.array([], dtype='@PYINT@')
Aij = np.array([], dtype='float32')
Bij = np.array([], dtype='float32')
for k in range(n):
for l in [k-1, k, k+1]:
if l < 0 or l > n-1:
continue
i = np.append(i, np.@PYINT@(k+1)) # Casting value on append is MANDATORY or C++ won't get the expected type.
j = np.append(j, np.@PYINT@(l+1)) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k:
Aij = np.append(Aij, np.float32( 200.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.float32( 33.3)) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k-1 or l == k+1:
Aij = np.append(Aij, np.float32(-100.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.float32( 16.6)) # Casting value on append is MANDATORY or C++ won't get the expected type.
for k, l, Akl in zip(i, j, Aij):
print("A[", k, ",", l, "] =", Akl)
for k, l, Bkl in zip(i, j, Bij):
print("B[", k, ",", l, "] =", Bkl)
A = (n, i, j, Aij) # coo format: dimension, i 1-based indices, j 1-based indices, Aij values.
B = (n, i, j, Bij) # coo format: dimension, i 1-based indices, j 1-based indices, Bij values.
# Get and tune arpack solver.
arpackSlv = pyarpackSlv.float() # Caution: float <=> np.array(..., dtype='float32')
arpackSlv.verbose = 3 # Set to 0 to get a quiet solve.
arpackSlv.debug = 1 # Set to 0 to get a quiet solve.
arpackSlv.nbEV = 2
arpackSlv.nbCV = 2*arpackSlv.nbEV + 1
arpackSlv.mag = 'LM'
arpackSlv.maxIt = 200
arpackSlv.slvOffset = 0.
arpackSlv.slvScale = 1.
# Solve eigen problem.
rc = arpackSlv.solve(A, B)
assert rc == 0, "bad solve"
rc = arpackSlv.checkEigVec(A, B, 1.e-2)
assert rc == 0, "bad checkEigVec"
# Print out results (mode selected, eigen vectors, eigen values, ...).
assert arpackSlv.nbEV == len(arpackSlv.val), "bad result"
print("\nresults:\n")
print("mode selected:", arpackSlv.mode)
print("nb iterations:", arpackSlv.nbIt)
print("Reverse Communication Interface time:", arpackSlv.rciTime, "s")
for val, vec in zip(arpackSlv.val, arpackSlv.vec):
print("eigen value:", val)
print("eigen vector:")
for v in range(n):
print(vec[v])
EXAMPLES/PYARPACK/pyarpackSparseLU.py.in
+117
View File
@@ -0,0 +1,117 @@
#!/usr/bin/env python
from __future__ import print_function
import numpy as np
from pyarpack import sparseLU as pyarpackSlv
# Build laplacian.
n = 4
i = np.array([], dtype='@PYINT@')
j = np.array([], dtype='@PYINT@')
Aij = np.array([], dtype='float64')
for k in range(n):
for l in [k-1, k, k+1]:
if l < 0 or l > n-1:
continue
i = np.append(i, np.@PYINT@(k)) # Casting value on append is MANDATORY or C++ won't get the expected type.
j = np.append(j, np.@PYINT@(l)) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k:
Aij = np.append(Aij, np.float64( 200.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k-1 or l == k+1:
Aij = np.append(Aij, np.float64(-100.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
for k, l, Akl in zip(i, j, Aij):
print("A[", k, ",", l, "] =", Akl)
A = (n, i, j, Aij) # coo format: dimension, i 0-based indices, j 0-based indices, Aij values.
# Get and tune arpack solver.
arpackSlv = pyarpackSlv.double() # Caution: double <=> np.array(..., dtype='float64')
arpackSlv.verbose = 3 # Set to 0 to get a quiet solve.
arpackSlv.debug = 1 # Set to 0 to get a quiet solve.
arpackSlv.nbEV = 1
arpackSlv.nbCV = 2*arpackSlv.nbEV + 1
arpackSlv.mag = 'LM'
arpackSlv.maxIt = 200
arpackSlv.slvPvtThd = 1.e-6
# Solve eigen problem.
rc = arpackSlv.solve(A)
assert rc == 0, "bad solve"
rc = arpackSlv.checkEigVec(A)
assert rc == 0, "bad checkEigVec"
# Print out results (mode selected, eigen vectors, eigen values, ...).
assert arpackSlv.nbEV == len(arpackSlv.val), "bad result"
print("\nresults:\n")
print("mode selected:", arpackSlv.mode)
print("nb iterations:", arpackSlv.nbIt)
print("Reverse Communication Interface time:", arpackSlv.rciTime, "s")
for val, vec in zip(arpackSlv.val, arpackSlv.vec):
print("eigen value:", val)
print("eigen vector:")
print(vec)
#######################################################################################
print("\n##########################################################################\n")
#######################################################################################
# Build laplacian.
n = 8
i = np.array([], dtype='@PYINT@')
j = np.array([], dtype='@PYINT@')
Aij = np.array([], dtype='float32')
Bij = np.array([], dtype='float32')
for k in range(n):
for l in [k-1, k, k+1]:
if l < 0 or l > n-1:
continue
i = np.append(i, np.@PYINT@(k+1)) # Casting value on append is MANDATORY or C++ won't get the expected type.
j = np.append(j, np.@PYINT@(l+1)) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k:
Aij = np.append(Aij, np.float32( 200.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.float32( 33.3)) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k-1 or l == k+1:
Aij = np.append(Aij, np.float32(-100.)) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.float32( 16.6)) # Casting value on append is MANDATORY or C++ won't get the expected type.
for k, l, Akl in zip(i, j, Aij):
print("A[", k, ",", l, "] =", Akl)
for k, l, Bkl in zip(i, j, Bij):
print("B[", k, ",", l, "] =", Bkl)
A = (n, i, j, Aij) # coo format: dimension, i 1-based indices, j 1-based indices, Aij values.
B = (n, i, j, Bij) # coo format: dimension, i 1-based indices, j 1-based indices, Bij values.
# Get and tune arpack solver.
arpackSlv = pyarpackSlv.float() # Caution: float <=> np.array(..., dtype='float32')
arpackSlv.verbose = 3 # Set to 0 to get a quiet solve.
arpackSlv.debug = 1 # Set to 0 to get a quiet solve.
arpackSlv.nbEV = 2
arpackSlv.nbCV = 2*arpackSlv.nbEV + 1
arpackSlv.mag = 'LM'
arpackSlv.maxIt = 200
arpackSlv.slvPvtThd = 1.e-6
arpackSlv.sigmaReal = 1
# Solve eigen problem.
rc = arpackSlv.solve(A, B)
assert rc == 0, "bad solve"
rc = arpackSlv.checkEigVec(A, B, 1.e-2)
assert rc == 0, "bad checkEigVec"
# Print out results (mode selected, eigen vectors, eigen values, ...).
assert arpackSlv.nbEV == len(arpackSlv.val), "bad result"
print("\nresults:\n")
print("mode selected:", arpackSlv.mode)
print("nb iterations:", arpackSlv.nbIt)
print("Reverse Communication Interface time:", arpackSlv.rciTime, "s")
for val, vec in zip(arpackSlv.val, arpackSlv.vec):
print("eigen value:", val)
print("eigen vector:")
for v in range(n):
print(vec[v])
EXAMPLES/PYARPACK/pyarpackSparseQR.py.in
+125
View File
@@ -0,0 +1,125 @@
#!/usr/bin/env python
from __future__ import print_function
import numpy as np
from pyarpack import sparseQR as pyarpackSlv
# Build laplacian.
n = 4
i = np.array([], dtype='@PYINT@')
j = np.array([], dtype='@PYINT@')
Aij = np.array([], dtype='complex128')
for k in range(n):
for l in [k-1, k, k+1]:
if l < 0 or l > n-1:
continue
i = np.append(i, np.@PYINT@(k)) # Casting value on append is MANDATORY or C++ won't get the expected type.
j = np.append(j, np.@PYINT@(l)) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k:
Aij = np.append(Aij, np.complex128(np.complex( 200., 200.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k-1:
Aij = np.append(Aij, np.complex128(np.complex(-101., -101.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k+1:
Aij = np.append(Aij, np.complex128(np.complex( -99., -99.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
for k, l, Akl in zip(i, j, Aij):
print("A[", k, ",", l, "] =", Akl)
A = (n, i, j, Aij) # coo format: dimension, i 0-based indices, j 0-based indices, Aij values.
# Get and tune arpack solver.
arpackSlv = pyarpackSlv.complexDouble() # Caution: complexDouble <=> np.array(..., dtype='complex128')
arpackSlv.verbose = 3 # Set to 0 to get a quiet solve.
arpackSlv.debug = 1 # Set to 0 to get a quiet solve.
arpackSlv.nbEV = 1
arpackSlv.nbCV = 2*arpackSlv.nbEV + 1
arpackSlv.mag = 'LM'
arpackSlv.maxIt = 200
arpackSlv.slvPvtThd = 1.e-6
arpackSlv.symPb = False
# Solve eigen problem.
rc = arpackSlv.solve(A)
assert rc == 0, "bad solve"
rc = arpackSlv.checkEigVec(A)
assert rc == 0, "bad checkEigVec"
# Print out results (mode selected, eigen vectors, eigen values, ...).
assert arpackSlv.nbEV == len(arpackSlv.val), "bad result"
print("\nresults:\n")
print("mode selected:", arpackSlv.mode)
print("nb iterations:", arpackSlv.nbIt)
print("Reverse Communication Interface time:", arpackSlv.rciTime, "s")
for val, vec in zip(arpackSlv.val, arpackSlv.vec):
print("eigen value:", val)
print("eigen vector:")
print(vec)
#######################################################################################
print("\n##########################################################################\n")
#######################################################################################
# Build laplacian.
n = 8
i = np.array([], dtype='@PYINT@')
j = np.array([], dtype='@PYINT@')
Aij = np.array([], dtype='complex64')
Bij = np.array([], dtype='complex64')
for k in range(n):
for l in [k-1, k, k+1]:
if l < 0 or l > n-1:
continue
i = np.append(i, np.@PYINT@(k+1)) # Casting value on append is MANDATORY or C++ won't get the expected type.
j = np.append(j, np.@PYINT@(l+1)) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k:
Aij = np.append(Aij, np.complex64(np.complex( 200., 200.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.complex64(np.complex( 33.3, 33.3))) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k-1:
Aij = np.append(Aij, np.complex64(np.complex(-101., -101.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.complex64(np.complex( 16.6, 16.6))) # Casting value on append is MANDATORY or C++ won't get the expected type.
if l == k+1:
Aij = np.append(Aij, np.complex64(np.complex( -99., -99.))) # Casting value on append is MANDATORY or C++ won't get the expected type.
Bij = np.append(Bij, np.complex64(np.complex( 16.6, 16.6))) # Casting value on append is MANDATORY or C++ won't get the expected type.
for k, l, Akl in zip(i, j, Aij):
print("A[", k, ",", l, "] =", Akl)
for k, l, Bkl in zip(i, j, Bij):
print("B[", k, ",", l, "] =", Bkl)
A = (n, i, j, Aij) # coo format: dimension, i 1-based indices, j 1-based indices, Aij values.
B = (n, i, j, Bij) # coo format: dimension, i 1-based indices, j 1-based indices, Bij values.
# Get and tune arpack solver.
arpackSlv = pyarpackSlv.complexFloat() # Caution: complexFloat <=> np.array(..., dtype='complex64')
arpackSlv.verbose = 3 # Set to 0 to get a quiet solve.
arpackSlv.debug = 1 # Set to 0 to get a quiet solve.
arpackSlv.nbEV = 2
arpackSlv.nbCV = 2*arpackSlv.nbEV + 1
arpackSlv.mag = 'LM'
arpackSlv.maxIt = 200
arpackSlv.slvPvtThd = 1.e-6
arpackSlv.sigmaReal = 1
arpackSlv.sigmaImag = 1
arpackSlv.symPb = False
# Solve eigen problem.
rc = arpackSlv.solve(A, B)
assert rc == 0, "bad solve"
rc = arpackSlv.checkEigVec(A, B, 1.e-2)
assert rc == 0, "bad checkEigVec"
# Print out results (mode selected, eigen vectors, eigen values, ...).
assert arpackSlv.nbEV == len(arpackSlv.val), "bad result"
print("\nresults:\n")
print("mode selected:", arpackSlv.mode)
print("nb iterations:", arpackSlv.nbIt)
print("Reverse Communication Interface time:", arpackSlv.rciTime, "s")
for val, vec in zip(arpackSlv.val, arpackSlv.vec):
print("eigen value:", val)
print("eigen vector:")
for v in range(n):
print(vec[v])
README.md
+1
View File
@@ -35,6 +35,7 @@ Important Features:
$ export INTERFACE64=1
$ ./configure --with-blas=mkl_gf_ilp64 --with-lapack=mkl_gf_ilp64
$ make all check```
* pyarpack: python support based on Boost.Python.Numpy exposing C++ API.
This project started as a joint project between Debian, Octave and Scilab in order to
provide a common and maintained version of arpack.