Shifted inverse iteration example

examples/CMakeLists.txt

@@ -2,6 +2,7 @@ add_subdirectory(3pt_stencil)
 add_subdirectory(asynchronous_stopping_criterion)
 add_subdirectory(custom_matrix_format)
 add_subdirectory(ginkgo_overhead)
+add_subdirectory(inverse_iteration)
 add_subdirectory(minimal_solver_cuda)
 add_subdirectory(preconditioned_solver)
 add_subdirectory(poisson_solver)

examples/inverse_iteration/CMakeLists.txt

@@ -0,0 +1,4 @@
+add_executable(inverse_iteration inverse_iteration.cpp)
+target_link_libraries(inverse_iteration ginkgo)
+target_include_directories(inverse_iteration PRIVATE ${PROJECT_SOURCE_DIR})
+configure_file(data/A.mtx data/A.mtx COPYONLY)

examples/inverse_iteration/build.sh

@@ -0,0 +1,37 @@
+#!/bin/bash
+
+# set up script
+if [ $# -ne 1 ]; then
+    echo -e "Usage: $0 GINKGO_BUILD_DIRECTORY"
+    exit 1
+fi
+BUILD_DIR=$1
+THIS_DIR=$( cd "$( dirname "${BASH_SOURCE[0]}" )" &>/dev/null && pwd )
+
+# copy libraries
+LIBRARY_DIRS="core core/device_hooks reference omp cuda"
+LIBRARY_NAMES="ginkgo ginkgo_reference ginkgo_omp ginkgo_cuda"
+SUFFIXES=".so .dylib .dll d.so d.dylib d.dll"
+for prefix in ${LIBRARY_DIRS}; do
+    for name in ${LIBRARY_NAMES}; do
+        for suffix in ${SUFFIXES}; do
+            cp ${BUILD_DIR}/${prefix}/lib${name}${suffix} \
+                ${THIS_DIR}/lib${name}${suffix} 2>/dev/null
+        done
+    done
+done
+
+# figure out correct compiler flags
+if ls ${THIS_DIR} | grep -F "libginkgo." >/dev/null; then
+    LINK_FLAGS="-lginkgo -lginkgo_omp -lginkgo_cuda -lginkgo_reference"
+else
+    LINK_FLAGS="-lginkgod -lginkgo_ompd -lginkgo_cudad -lginkgo_referenced"
+fi
+if [ -z "${CXX}" ]; then
+    CXX="c++"
+fi
+
+# build
+${CXX} -std=c++11 -o \
+    ${THIS_DIR}/inverse_iteration ${THIS_DIR}/inverse_iteration.cpp \
+    -I${THIS_DIR}/../.. -L${THIS_DIR} ${LINK_FLAGS}

examples/inverse_iteration/data/A.mtx

@@ -0,0 +1,114 @@
+%%MatrixMarket matrix coordinate integer symmetric
+%-------------------------------------------------------------------------------
+% UF Sparse Matrix Collection, Tim Davis
+% http://www.cise.ufl.edu/research/sparse/matrices/JGD_Trefethen/Trefethen_20b
+% name: JGD_Trefethen/Trefethen_20b
+% [Diagonal matrices with primes, Nick Trefethen, Oxford Univ.]
+% id: 2203
+% date: 2008
+% author: N. Trefethen
+% ed: J.-G. Dumas
+% fields: name title A id date author ed kind notes
+% kind: combinatorial problem
+%-------------------------------------------------------------------------------
+% notes:
+% Diagonal matrices with primes, Nick Trefethen, Oxford Univ.          
+% From Jean-Guillaume Dumas' Sparse Integer Matrix Collection,         
+% http://ljk.imag.fr/membres/Jean-Guillaume.Dumas/simc.html            
+%                                                                      
+% Problem 7 of the Hundred-dollar, Hundred-digit Challenge Problems,   
+% SIAM News, vol 35, no. 1.                                            
+%                                                                      
+% 7. Let A be the 20,000 x 20,000 matrix whose entries are zero        
+% everywhere except for the primes 2, 3, 5, 7, . . . , 224737 along the
+% main diagonal and the number 1 in all the positions A(i,j) with      
+% |i-j| = 1,2,4,8, . . . ,16384.  What is the (1,1) entry of inv(A)?   
+%                                                                      
+% http://www.siam.org/news/news.php?id=388                             
+%                                                                      
+% Filename in JGD collection: Trefethen/trefethen_20__19_minor.sms     
+%-------------------------------------------------------------------------------
+19 19 83
+1 1 3
+2 1 1
+3 1 1
+5 1 1
+9 1 1
+17 1 1
+2 2 5
+3 2 1
+4 2 1
+6 2 1
+10 2 1
+18 2 1
+3 3 7
+4 3 1
+5 3 1
+7 3 1
+11 3 1
+19 3 1
+4 4 11
+5 4 1
+6 4 1
+8 4 1
+12 4 1
+5 5 13
+6 5 1
+7 5 1
+9 5 1
+13 5 1
+6 6 17
+7 6 1
+8 6 1
+10 6 1
+14 6 1
+7 7 19
+8 7 1
+9 7 1
+11 7 1
+15 7 1
+8 8 23
+9 8 1
+10 8 1
+12 8 1
+16 8 1
+9 9 29
+10 9 1
+11 9 1
+13 9 1
+17 9 1
+10 10 31
+11 10 1
+12 10 1
+14 10 1
+18 10 1
+11 11 37
+12 11 1
+13 11 1
+15 11 1
+19 11 1
+12 12 41
+13 12 1
+14 12 1
+16 12 1
+13 13 43
+14 13 1
+15 13 1
+17 13 1
+14 14 47
+15 14 1
+16 14 1
+18 14 1
+15 15 53
+16 15 1
+17 15 1
+19 15 1
+16 16 59
+17 16 1
+18 16 1
+17 17 61
+18 17 1
+19 17 1
+18 18 67
+19 18 1
+19 19 71

examples/inverse_iteration/inverse_iteration.cpp

@@ -0,0 +1,217 @@
+/*******************************<GINKGO LICENSE>******************************
+Copyright 2017-2018
+
+Karlsruhe Institute of Technology
+Universitat Jaume I
+University of Tennessee
+
+Redistribution and use in source and binary forms, with or without modification,
+are permitted provided that the following conditions are met:
+
+1. Redistributions of source code must retain the above copyright notice,
+   this list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright notice,
+   this list of conditions and the following disclaimer in the documentation
+   and/or other materials provided with the distribution.
+
+3. Neither the name of the copyright holder nor the names of its contributors
+   may be used to endorse or promote products derived from this software
+   without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
+ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
+ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+******************************<GINKGO LICENSE>*******************************/
+
+/*****************************<COMPILATION>***********************************
+The easiest way to build the example solver is to use the script provided:
+./build.sh <PATH_TO_GINKGO_BUILD_DIR>
+
+Ginkgo should be compiled with `-DBUILD_REFERENCE=on` option.
+
+Alternatively, you can setup the configuration manually:
+
+Go to the <PATH_TO_GINKGO_BUILD_DIR> directory and copy the shared
+libraries located in the following subdirectories:
+
+    + core/
+    + core/device_hooks/
+    + reference/
+    + omp/
+    + cuda/
+
+to this directory.
+
+Then compile the file with the following command line:
+
+c++ -std=c++11 -o inverse_iteration inverse_iteration.cpp -I../.. \
+    -L. -lginkgo -lginkgo_reference -lginkgo_omp -lginkgo_cuda
+
+(if ginkgo was built in debug mode, append 'd' to every library name)
+
+Now you should be able to run the program using:
+
+env LD_LIBRARY_PATH=.:${LD_LIBRARY_PATH} ./inverse_iteration
+
+*****************************<COMPILATION>**********************************/
+
+/*****************************<DECSRIPTION>***********************************
+This example shows how components available in Ginkgo can be used to implement
+higher-level numerical methods. The method used here will be the shifted inverse
+iteration method for eigenvalue computation which find the eigenvalue and
+eigenvector of A closest to z, for some scalar z. The method requires repeatedly
+solving the shifted linear system (A - zI)x = b, as well as performing
+matrix-vector products with the matrix `A`. Here is the complete pseudocode of
+the method:
+
+```
+x_0 = initial guess
+for i = 0 .. max_iterations:
+    solve (A - zI) y_i = x_i for y_i+1
+    x_(i+1) = y_i / || y_i ||      # compute next eigenvector approximation
+    g_(i+1) = x_(i+1)^* A x_(i+1)  # approximate eigenvalue (Rayleigh quotient)
+    if ||A x_(i+1) - g_(i+1)x_(i+1)|| < tol * g_(i+1):  # check convergence
+        break
+```
+ *****************************<DECSRIPTION>**********************************/
+
+#include <include/ginkgo.hpp>
+
+
+#include <cmath>
+#include <complex>
+#include <fstream>
+#include <iomanip>
+#include <iostream>
+#include <string>
+
+
+int main(int argc, char *argv[])
+{
+    // Some shortcuts
+    using precision = std::complex<double>;
+    using real_precision = double;
+    using vec = gko::matrix::Dense<precision>;
+    using mtx = gko::matrix::Csr<precision>;
+    using solver_type = gko::solver::Bicgstab<precision>;
+
+    using std::abs;
+    using std::sqrt;
+
+    // Print version information
+    std::cout << gko::version_info::get() << std::endl;
+
+    std::cout << std::scientific << std::setprecision(8) << std::showpos;
+
+    // Figure out where to run the code
+    std::shared_ptr<gko::Executor> exec;
+    if (argc == 1 || std::string(argv[1]) == "reference") {
+        exec = gko::ReferenceExecutor::create();
+    } else if (argc == 2 && std::string(argv[1]) == "omp") {
+        exec = gko::OmpExecutor::create();
+    } else if (argc == 2 && std::string(argv[1]) == "cuda" &&
+               gko::CudaExecutor::get_num_devices() > 0) {
+        exec = gko::CudaExecutor::create(0, gko::OmpExecutor::create());
+    } else {
+        std::cerr << "Usage: " << argv[0] << " [executor]" << std::endl;
+        std::exit(-1);
+    }
+
+    auto this_exec = exec->get_master();
+
+    // linear system solver parameters
+    auto system_max_iterations = 100u;
+    auto system_residual_goal = real_precision{1e-16};
+
+    // eigensolver parameters
+    auto max_iterations = 20u;
+    auto residual_goal = real_precision{1e-8};
+    auto z = precision{20.0, 2.0};
+
+    // Read data
+    auto A = share(gko::read<mtx>(std::ifstream("data/A.mtx"), exec));
+
+    // Generate shifted matrix  A - zI
+    // - we avoid duplicating memory by not storing both A and A - zI, but
+    //   compute A - zI on the fly by using Ginkgo's utilities for creating
+    //   linear combinations of operators
+    auto one = share(gko::initialize<vec>({precision{1.0}}, exec));
+    auto neg_one = share(gko::initialize<vec>({-precision{1.0}}, exec));
+    auto neg_z = gko::initialize<vec>({-z}, exec);
+
+    auto system_matrix = share(gko::Combination<precision>::create(
+        one, A, gko::initialize<vec>({-z}, exec),
+        gko::matrix::Identity<precision>::create(exec, A->get_size()[0])));
+
+    // Generate solver operator  (A - zI)^-1
+    auto solver =
+        solver_type::Factory::create()
+            .with_criterion(
+                gko::stop::Combined::Factory::create()
+                    .with_criteria(
+                        gko::stop::Iteration::Factory::create()
+                            .with_max_iters(system_max_iterations)
+                            .on_executor(exec),
+                        gko::stop::ResidualNormReduction<
+                            precision>::Factory::create()
+                            .with_reduction_factor(system_residual_goal)
+                            .on_executor(exec))
+                    .on_executor(exec))
+            .on_executor(exec)
+            ->generate(system_matrix);
+
+    // inverse iterations
+
+    // start with guess [1, 1, ..., 1]
+    auto x = [&] {
+        auto work = vec::create(this_exec, gko::dim<2>{A->get_size()[0], 1});
+        const auto n = work->get_size()[0];
+        for (int i = 0; i < n; ++i) {
+            work->get_values()[i] = precision{1.0} / sqrt(n);
+        }
+        return clone(exec, work);
+    }();
+    auto y = clone(x);
+    auto tmp = clone(x);
+    auto norm = clone(one);
+    auto inv_norm = clone(this_exec, one);
+    auto g = clone(one);
+
+    for (auto i = 0u; i < max_iterations; ++i) {
+        std::cout << "{ ";
+        // (A - zI)y = x
+        solver->apply(lend(x), lend(y));
+        system_matrix->apply(lend(one), lend(y), lend(neg_one), lend(x));
+        x->compute_dot(lend(x), lend(norm));
+        std::cout << "\"system_residual\": "
+                  << sqrt(clone(this_exec, norm)->get_values()[0]) << ", ";
+        x->copy_from(lend(y));
+        // x = y / || y ||
+        x->compute_dot(lend(x), lend(norm));
+        inv_norm->get_values()[0] =
+            precision{1.0} / sqrt(clone(this_exec, norm)->get_values()[0]);
+        x->scale(lend(clone(exec, inv_norm)));
+        // g = x^* A x
+        A->apply(lend(x), lend(tmp));
+        x->compute_dot(lend(tmp), lend(g));
+        auto g_val = clone(this_exec, g)->get_values()[0];
+        std::cout << "\"eigenvalue\": " << g_val << ", ";
+        // ||Ax - gx|| < tol * g
+        auto v = gko::initialize<vec>({-g_val}, exec);
+        tmp->add_scaled(lend(v), lend(x));
+        tmp->compute_dot(lend(tmp), lend(norm));
+        auto res_val = sqrt(clone(exec->get_master(), norm)->get_values()[0]);
+        std::cout << "\"residual\": " << res_val / g_val << " }," << std::endl;
+        if (abs(res_val) < residual_goal * abs(g_val)) {
+            break;
+        }
+    }
+}

include/ginkgo.hpp

@@ -69,6 +69,7 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 #include "core/solver/cg.hpp"
 #include "core/solver/cgs.hpp"
 #include "core/solver/fcg.hpp"
+#include "core/solver/gmres.hpp"
 
 #include "core/stop/combined.hpp"
 #include "core/stop/iteration.hpp"