README
This folder contains the source code of DCGeig, a collection of dense and sparse
solvers for generalized eigenvalue problems with real symmetric positive
semidefinite matrices.
DCGeig is based on the source code for the master's thesis "Projection Methods
for Generalized Eigenvalue Problems" by Christoph Conrads.
INSTALLATION
You need the following software:
- Python 2.7 or newer
- Cython
- NumPy
- SciPy
- BLAS
- LAPACK 3.6.0 or newer
- GCC C compiler
- GCC C++ compiler 4.7 or newer
- cmake 3.2 or newer
- METIS 5.0 or newer
- git
- bash
BLAS and LAPACK can be provided by OpenBLAS or by Intel MKL 11.2 or newer and
instructions for this build are given below. Note that NumPy and SciPy need to
be linked against Intel MKL, too, if you want to see any speed-up. Instructions
for such a build can be found here:
https://christoph-conrads.name/building-numpy-and-scipy-with-intel-compilers-and-intel-mkl/
All instructions below are to be executed in bash.
Installation with BLAS and LAPACK:
- Locate the source code. In the following, we assume the bash variable
"sourcedir" contains the path to the directory with the source code. Test if
the variable is set correctly by typing
echo "${sourcedir}"
in bash (variables != environment variables).
- Decide where to install and save the path in a variable, e.g., if you want to
install the solvers in your home directory, type
prefixdir=$(echo ~)
- Create a new directory for the build and change into it:
builddir=$(mktemp -d)
cd "${builddir}"
- Configure the build with cmake:
cmake \
-DCMAKE_INSTALL_PREFIX="${prefixdir}" \
"${sourcedir}"
- Build the solvers:
make -j
- Install the solvers:
make install
Installation with Intel MKL
- Proceed as in the installation with BLAS and LAPACK until you have to
configure the build with cmake
- Ensure the Intel MKL environment variable MKLROOT is set:
echo "${MKLROOT}"
- Configure the build with cmake:
cmake \
-DCMAKE_INSTALL_PREFIX="${prefixdir}" \
-DBLAS=mkl -DMKL_ROOT="$MKLROOT" \
"${sourcedir}"
- Build the solvers:
make -j
- Install the solvers:
make install
Updating Search Paths
Some of the software is written in Python 2. The Python interpreter has a
built-in list of paths where it searches for Python packages. Similarly, the
linker tries to find the library containing the dense solvers in certain
standard directories. If you installed the solvers in a non-standard location,
e.g., your home directory, then you need to set the environment variables
LD_LIBRARY_PATH and PYTHONPATH *whenever* you want to run the solvers:
export PYTHONPATH=$PYTHONPATH:${prefixdir}/lib/python2.7/site-packages/
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:${prefixdir}/lib
Furthermore, the variable LD_LIBRARY_PATH needs to contain the path to the Intel
MKL libraries.
You can test the Python code with the following command:
python2 -m unittest discover dcgeig.tests
PROGRAMS
multilevel-gep-solver:
This program solves large, sparse generalized eigenvalue problems with
Hermitian positive definite matrices. Its usage can be queried by calling
multilevel-gep-solver --help
There are two obligatory arguments. The first argument is the path to a file
in Matrix Market format (.mtx) containing the stiffness matrix of a matrix
pencil. The file can be compressed with gzip or bzip2 (.mtx.gz, .mtx.bz2). In
order to find the corresponding mass matrix, the first occurence of the letter
'k' will be replaced by the letter 'm'. If no mass matrix can be found, the
solver assumes it deals with a standard eigenvalue problem. The second
argument is the cutoff value: the solver tries to find all eigenpairs
(\lambda, x) where \lambda is less than or equal to the cutoff.
The '--stats ' causes the solver to print subproblems statistics in the
given file. Below a whole section explains the contents of these files.
benchmark-gep-solver:
Given a file with test problems and cutoff values, this program repeatedly
calls the multilevel-gep-solver saving its (standard) output and its --stats
output. Call
benchmark-gep-solver
to use it. The key is a user-provided sequence of letters and digits used in
the filenames. For example, given the key 'foo123', the standard output will
be saved in the file 'foo123.out' and the --stats output of the problem
bcsstk13.mtx will be saved in 'bcsstk13.mtx.foo123.log'.
solve-gep:
This program solves a dense GEP with all four dense solvers presented in the
thesis:
- the standard solver (reduction to a standard eigenvalue problem);
- deflation of infinite eigenvalues, then calling the standard solver;
- direct GSVD;
- GSVD via QR factorizations and CS decomposition.
A unitary congruence transformation is applied before calling the GEP solvers.
This operation is supposed to avoid diagonal mass and matrices and to force
the GEP solver to decide on the rank of the mass matrix (rank decisions are
trivial with diagonal matrices).
For every solver, the program prints
- the name of the problem,
- the dimension of the problem $n$,
- the solver name,
- the solver time in seconds,
- the maximum backward error divided by n*eps.
To call the program, use
solve-gep [precision]
The stiffness matrix has to be in Matrix Market format (.mtx) and can be
compressed with gzip (.mtx.gz) or with bzip2 (.mtx.bz2). If a mass matrix
cannot be found, the program assumes it has to solve a standard eigenvalue
problem. The optional 'precision' option default to 'double'; you can also set
'single'.
With Netlib BLAS and LAPACK, solving problems with 3600 degrees of freedom by
means of direct GSVD takes several hours on the author's computer. For small
problems (n=1,2,3,4), the largest backward error of the solutions computed by
GSVD-based solvers may be larger than n*eps.
STATS OUTPUT
This section explains the columns in the stats generated by the '--stats'
argument to multilevel-gep-solver. Every row contains data about the subproblems
encountered by the solver. A subproblem is considered solved if every eigenpair
(\lambda, x), \lambda <= \lambda_c, is considered accurate. An eigenpair is
considered accurate if the backward error is less than the single precision
epsilon and if the relative forward error is less than one. An eigenpair
(\lambda, x) is called "desired" if \lambda \leq \lambda_c.
- pid
The id of the problem at hand. For block diagonal matrices, each pair of
blocks on the diagonal is a single problem.
- sid (subproblem id)
A unique id generated by a post-order traversal of the partitioning tree for
the problem at hand. Consequently the lowest id belongs to a subproblem
that will be solved directly while the largest id is the "subproblem" for
finding eigenpairs of the problem at hand.
- lvl
The depth in the recursion tree. The complete GEP is at level 0.
- n
The dimension of the subproblem
- n_c
The number of eigenpairs (\lambda, x) with \lambda <= \lambda_c *after*
solving the subproblem. If there is no such eigenpair, then n_c will be set
to one so that at least one eigenpair has to fulfill the convergence
criterion.
- n_s
The size of the search space *after* solving the subproblem. The size of the
search space *before* solving the subproblem can be calculated as follows:
- If the subproblem was solved directly, then the initial search space was
all of the $n$-dimension space, i.e., n = n_s.
- Otherwise, the problem possesses a nested dissection ordering structure
with two substructure blocks and one coupling block. The initial search
space size is the sum of the search space returned by the two subproblems
in the structure blocks plus the size of the coupling block.
n_s is always at least 32 in order to avoid empty search spaces.
- fill
The fill-in by the Cholesky decomposition of the stiffness matrix computed
for the subspace iteration method.
- norm:K
The Frobenius norm of the stiffness matrix K.
- norm:M
The Frobenius norm of the mass matrix M.
- normK12
The Frobenius norm of the matrix K12, where K is partitioned as 2x2 block
matrix.
- normM12
The Frobenius norm of the matrix M12, where M is partitioned as 2x2 block
matrix.
- min:ev
The minimum eigenvalue in the search space *after* solving a subproblem
divided by lambda_c.
- max:ev
The largest eigenvalue in the search space *after* solving a subproblem
divided by lambda_c.
- min:be
The smallest backward error of all eigenpairs.
- max:be
The largest backward error of all eigenpairs.
- min:fe
The smallest relative forward error of all eigenpairs.
- iter
The number of subspace iterations until all desired eigenpairs fulfilled the
convergence criterion.
- t-sle
The cumulative wall-clock time spent on solving systems of linear equations
Kx=b, during subspace iteration.
- t-rr
The cumulative wall-clock time spent on the Rayleigh-Ritz procedure during
subspace iteration.
- t-wc
If the subproblem is solved directly, then this is the wall-clock time it took
to solve the subproblem directly. Otherwise, divide-and-conquer (divide,
conquer, combine) is applied. In this case, t-wc refers to the wall-clock
time for the combine phase. The divide phase is ignored because it is
comparatively cheap. It holds that t-wc >= t-sle + t-rr.
- t-cpu
Like t-wc but with CPU time instead of wall-clock time.
CREDITS
Christoph Conrads