AMDiS/lib/mtl4/MTLConfig.cmake

+set(MTL_COMMON_CONFIG "${MTL_DIR}/tools/cmake/MTLCommon.cmake")
+if( EXISTS ${MTL_COMMON_CONFIG})
+	include(${MTL_COMMON_CONFIG})
+else()
+	message(ERROR "could not find the common mtl configuration. this is possibly a wrong installation of mtl.")
+endif()

AMDiS/lib/mtl4/README

+-*- text -*-
+
+Software License for MTL
+
+Copyright (c) 2007-2008 The Trustees of Indiana University. All rights reserved.
+Authors: Peter Gottschling and Andrew Lumsdaine
+
+This file is part of the Matrix Template Library
+
+See also license.mtl.txt in the distribution.
+---------------------------------------------------------------
+
+              The Matrix Template Library
+                      Version 4
+
+I. Introduction
+===============
+
+The Matrix Template Library is a C++ class library for basic linear
+algebra. The MTL is designed for high-performance while at the same
+time taking advantage of the generic programming paradigm (ala the
+STL) to allow much greater flexibility and breadth of
+functionality. Many new and advanced programming techniques were used
+in the construction of this library.
+
+The MTL is a low level library in the sense that the user must be
+conscious of the matrix type being used, and that all computationally
+expensive operations are explicit. The MTL is not a C++ Matlab or
+Mathematica. Nevertheless, the interface is designed to be simple and easy
+to use.
+
+The matrix types provided include compressed sparse row/column,
+traditional row- and column-major dense, Morton-order, and block-recursive
+matrices. All matrix types share a common and easy to use interface.
+
+The algorithms consist of the traditional basic linear algebra
+routines (from the BLAS level-1 to 3) which includes matrix and vector
+arithmetic.
+
+
+II. Compilers
+=============
+
+The Matrix Template Library is written in compliance with the C++ standard
+and should be compilable with every compiler compliant with the standard.
+It has been tested (and passed) with the following compilers and architectures:
+
+Linux & MacOS
+-------------
+    clang++ 3.0 (32 bit)
+    g++ 4.1.3 (32 and 64 bit)
+    g++ 4.2.4 (64bit)
+    g++ 4.3.4 (32 and 64 bit)
+    g++ 4.4.3 (32 and 64 bit)
+    g++ 4.5.1 (64 bit)
+    g++ 4.6.0 (64 bit)
+    g++ 4.6.1 (64 bit)
+    icc 9.1 (32 and 64 bit)
+    icc 10.0 (32 and 64 bit)
+    icc 10.1 (32 and 64 bit)
+    icc 11.0 (32 and 64 bit)
+    icc 11.1 (32 and 64 bit)
+
+Windows
+-------
+    VC 8.0 from Visual Studio 2005 (32 bit)
+    VC 9.0 from Visual Studio 2008 (32 bit)
+    VC 10.0 from Visual Studio 2010 (32 bit)
+
+Compilers that are not standard-compliant (e.g. VC 6.0 from VS 2003) are not subject to support.
+
+
+III. Documentation
+==================
+
+The documentation can be generated with doxygen in the main directory and
+is found in the /libs/numeric/mtl/doc. You
+can also view the documentation at our web site:
+http://www.simunova.com/node/184
+At the moment, the documentation illustrates the basic operations and
+partially provides a reference manual of some classes and functions.
+A detailed description of more features, especially tuning opportunities,
+is in the works.
+
+There are also a large number of examples in the /libs/numeric/mtl/examples
+directory.
+
+
+IV. Installation Procedure
+==========================
+
+See the file INSTALL in this directory.
+
+
+
+V. Contact Information
+=======================
+
+The Matrix Template Library is available at the main distribution site:
+
+    http://www.mtl4.org
+
+This distribution includes:
+
+    - Source code for the MTL (all header files /boost/numeric/mtl)
+    - Some completing sources (header files /boost/numeric/meta_math)
+    - A test suite (/libs/numeric/mtl/test)
+    - Some example code (/libs/numeric/mtl/examples)
+
+Bug reports should be sent to mtl4@osl.iu.edu.
+
+Questions, comments, suggestions, and requests for additional
+information should also be directed to mtl4@osl.iu.edu.
+
+Subscribe at:
+
+     http://www.osl.iu.edu/mailman/listinfo.cgi/mtl4

AMDiS/lib/mtl4/SConstruct

+import getopt, sys, os, os.path, re, string
+
+######################
+# Helper functions
+######################
+
+# platform-dependent preprocessor definitions
+def macro_flag(env, flag):
+    if sys.platform == 'win32' and conf.env['CC'] == 'cl':
+        return '/D' + flag
+    else:
+        return '-D' + flag
+
+# add debugging flags to environment
+def add_debug(env):
+    if sys.platform == 'win32' and conf.env['CC'] == 'cl':
+        env.Append(CCFLAGS = '/Z7')
+        env.Append(CCFLAGS = '/MDd')
+    else:
+        env.Append(CCFLAGS = '-g')
+    #print macro_flag(env, 'MTL_ASSERT_FOR_THROW')
+    env.Append(CCFLAGS = macro_flag(env, 'MTL_ASSERT_FOR_THROW'))
+
+
+def check_g5():
+    print "platform ", sys.platform
+    if sys.platform == "darwin" :
+
+	system("/usr/bin/system_profiler SPHardwareDataType")
+	profile = os.popen("system_profiler SPHardwareDataType").read()
+	print profile
+	if profile.find("Mac G5") != -1:
+	    print "is G5"
+	    return 1
+    print "isn't G5"
+    return None
+
+
+def add_platform_flags(env):
+    if env['CC'] == 'cl':
+        # Turn off some warnings: 1. for problems with CRTP 2. with VC's own headers
+        #                         3. alleged ambiguos assignment
+        # I hope we find a cleaner way to get rid of this warning
+        wd = ' /wd4355 /wd4996 /wd4522'
+        env.Append(CCFLAGS = '/EHa /D_CONSOLE /D_CRT_SECURE_NO_DEPRECATE /DNOMINMAX' +
+			     ' /D_SECURE_SCL=0' + wd)
+
+        # The following two if shouldn't be necessary
+        # Note that INCLUDE must be split and LIB must not :-P
+        if os.environ.has_key('INCLUDE') :
+            for p in os.environ['INCLUDE'].split(';'):
+                env.Append(CPPPATH = p)
+        if os.environ.has_key('LIB') :
+            env.Append(LIBPATH = os.environ['LIB'])
+
+
+# add optimization flags to environment
+# normal optimization for opt=1, for opt=2 high optimization
+# at the moment additional flags are added in this function
+# -> this needs to be cleaned up!
+def add_opt(env, opt):
+    if int(opt):
+        env.Append(CCFLAGS = macro_flag(env, 'NDEBUG'))
+    if sys.platform == 'win32' and conf.env['CC'] == 'cl':
+        if int(opt) == 1:
+            env.Append(CCFLAGS = '/O2 /Ot /Ob2')
+        elif int(opt) == 2:
+            env.Append(CCFLAGS = '/Ox')
+            env.Append(CCFLAGS = '/MD')
+    else:
+        if int(opt) == 2:
+    	    env.Append(CCFLAGS = '-O3')
+    	    if conf.env['CC'] == 'gcc':
+    	        env.Append(CCFLAGS = '-ffast-math')
+    	        #check_g5()
+            # env.Append(CCFLAGS = '-pg') # for profiling
+    	    # env.Append(LINKFLAGS = '-pg') # for profiling
+        elif int(opt) == 1:
+    	    env.Append(CCFLAGS = '-O2')
+            # env.Append(CCFLAGS = '-pg') # for profiling
+    	    # env.Append(LINKFLAGS = '-pg') # for profiling
+
+    if env['CC'] != 'cl':
+        env.Append(CXXFLAGS = conf.env['CCFLAGS'])
+
+
+def check_no_long_double(conf):
+    cc = conf.env['CC']
+    if cc == 'gcc':
+        # conf.env.Append(CXXFLAGS = '-Werror -Wreorder -pedantic -Wpointer-arith -Wcast-align -Wcast-qual -Wwrite-strings')
+	output = os.popen(cc + " --version").read()
+	st = output.split('\n')[0]
+	tmatch = re.search(".* \(GCC\) (\d\.\d)", st)
+	if tmatch and string.atof(tmatch.group(1)) == 4.0:
+	    conf.env.Append(CCFLAGS = '-Wno-long-double')
+	    # conf.env.Append(CXXFLAGS = '-Wno-long-double')
+
+
+def add_user_opt(env, add_optflag):
+     env.Append(CCFLAGS = add_optflag)
+     env.Append(CXXFLAGS = add_optflag)
+
+###############################
+# check for installed blas ...
+# search different libraries
+###############################
+def detected_blas(env): # blas library found
+   env.Append(CPPDEFINES = 'MTL_HAS_BLAS')
+
+def check_for_blas(env):
+ # symbolname
+ symname = 'dgemm_'
+ found_blas = 0
+
+ # get command line blas path
+ blas_path = ARGUMENTS.get('blas_path', '')
+ if blas_path:
+   print 'adding ' + blas_path + ' to LIBPATH.'
+   env.Append(LIBPATH = [ blas_path, blas_path + '/lib/' ])
+
+ # extra linker flags for blas
+ blas_ldflags = ARGUMENTS.get('blas_ldflags', '')
+ if blas_ldflags:
+   env.Append(_LIBFLAGS=blas_ldflags.split())
+
+ # get command line blas lib
+ blas_lib = ARGUMENTS.get('blas_lib', '')
+ if blas_lib:
+   myenv = env.Copy()
+   conf = Configure(myenv)
+   # check supplied lib
+   print "Checking for lib " + blas_lib + "..."
+   if conf.CheckLib(blas_lib, symname):
+     myenv = conf.Finish()
+     detected_blas(myenv)
+     myenv.Append(LIBS=[blas_lib])
+     return myenv
+ else:
+   print "Autodetecting BLAS support. See config.log for details!"
+
+   ########################
+   # check for acml
+   myenv = env.Copy()
+   # new configure object
+   conf = Configure(myenv)
+   # additional libs needed for acml
+   myenv.Append(LIBS=['m', 'g2c'])
+   if conf.CheckLib('acml', symname):
+     detected_blas(myenv)
+     found_blas = 1
+   myenv = conf.Finish()
+   if(found_blas == 1):
+     return myenv
+
+   ########################
+   # check for goto
+   myenv = env.Copy()
+   # new configure object
+   conf = Configure(myenv)
+   # additional libs needed for goto
+   myenv.Append(LIBS=['pthread'])
+   # extra linker flags for libgoto
+   # myenv.Append(_LIBFLAGS=['libs/numeric/mtl/build/xerbla.c'])
+   myenv.Append(_LIBFLAGS=['libs/numeric/mtl/build/xerbla.o']) # not portable !!!
+   # myenv.Library('build/xerbla', 'build/xerbla.c')
+   # myenv.Append(LIBS=['xerbla'])
+   # myenv.Append(LIBPATH=['build'])
+
+   if conf.CheckLib('goto', symname) or conf.CheckLib('goto_opteron-64', symname) or \
+      conf.CheckLib('goto_coppermine-32', symname) or conf.CheckLib('goto_opteron-32', symname) or \
+      conf.CheckLib('goto_itanium2-64', symname) or conf.CheckLib('goto_katmai-32', symname) or \
+      conf.CheckLib('goto_northwood-32', symname) or conf.CheckLib('goto_prescott-32', symname) or \
+      conf.CheckLib('goto_prescott-64', symname):
+     detected_blas(myenv)
+     found_blas = 1
+   myenv = conf.Finish()
+   if(found_blas == 1):
+     return myenv
+
+   ########################
+   # check for ATLAS
+   myenv = env.Copy()
+   # new configure object
+   conf = Configure(myenv)
+   # additional libs needed for goto
+   myenv.Append(LIBS=['f77blas', 'g2c'])
+   if conf.CheckLib('atlas', symname):
+     detected_blas(myenv)
+     found_blas = 1
+   myenv = conf.Finish()
+   if(found_blas == 1):
+     return myenv
+
+ return env
+
+
+
+###################################
+# Add UMFPACK (extremely simplistic)
+###################################
+
+def detected_umfpack(env): # umfpack library found
+   env.Append(CPPDEFINES = 'MTL_HAS_UMFPACK')
+
+def check_for_umfpack(env):
+ ufconfig_path = ARGUMENTS.get('ufconfig_path', '')
+ if ufconfig_path:
+   env.Prepend(CPPPATH = [ ufconfig_path ])
+ else:
+   print "UMFPACK only works with the UFconfig library."
+   print "Provide ufconfig_path."
+   print "Building continued without UMFPACK."
+   return env
+
+ amd_path = ARGUMENTS.get('amd_path', '')
+ if amd_path:
+   # print 'adding ' + amd_path + '/Lib to LIBPATH.'
+   env.Append(LIBPATH = [ amd_path + '/Lib' ])
+   #print 'adding ' + amd_path + '/Include to CPPPATH.'
+   env.Prepend(CPPPATH = [ amd_path + '/Include' ])
+   env.Prepend(LIBS=['amd'])
+   # print "nachher: ", env['CPPPATH']
+ else:
+   print "UMFPACK only works with the AMD library."
+   print "Provide amd_path."
+   print "Building continued without UMFPACK."
+   return env
+
+ # get command line umfpack path
+ umfpack_path = ARGUMENTS.get('umfpack_path', '')
+ if umfpack_path:
+   # print 'adding ' + umfpack_path + '/Lib to LIBPATH.'
+   env.Prepend(LIBPATH = [ umfpack_path + '/Lib' ])
+   env.Prepend(LIBS=['umfpack'])
+   # print 'adding ' + umfpack_path + '/Include to CPPPATH.'
+   env.Append(CPPPATH = [ umfpack_path + '/Include' ])
+   detected_umfpack(env)
+   return env
+
+
+
+
+######################
+# Start setting up env
+######################
+
+
+SourceSignatures('timestamp')
+
+env = Environment()
+# env.Decider('timestamp-newer')
+
+conf = Configure(env)
+check_no_long_double(conf)
+env = conf.Finish()
+
+# Search only for BLAS when explicitly asked for
+with_blas = ARGUMENTS.get('with_blas', 0)
+if int(with_blas):
+    env = check_for_blas(env)
+
+# Search only for UMFPACK when explicitly asked for
+with_umfpack = ARGUMENTS.get('with_umfpack', 0)
+if int(with_umfpack):
+    env = check_for_umfpack(env)
+
+
+
+
+######################
+# Include paths
+######################
+
+
+my_cpppath = []
+if os.environ.has_key('MTL_BOOST_ROOT') :
+    mtlp = os.environ['MTL_BOOST_ROOT']
+    my_cpppath.append(mtlp)
+if os.environ.has_key('BOOST_ROOT') :
+    my_cpppath.append(os.environ['BOOST_ROOT'])
+
+# In branches, the current directory is prepended
+pwd = os.getcwd()
+if not os.environ.has_key('MTL_BOOST_ROOT') \
+or not hasattr(os.path, 'samefile') or not os.path.samefile(pwd, mtlp):
+     my_cpppath = [pwd] + my_cpppath
+
+# Add directory of extern libraries if exist
+# Not needed currently
+#if os.path.exists('extern'):
+#    my_cpppath.append(os.path.abspath('extern'))
+
+env.Append(CPPPATH = my_cpppath)
+
+# Hack to select compiler explicitly
+# env['CXX'] = 'g++-4.1'
+
+
+######################
+# Opt. and debug flags
+######################
+
+
+opts = Options()
+opts.Add('opt', 'Set to 1 for normal optimization, 2 for high optimization', 0)
+opts.Add('debug', 'Set to 1 for debugging', 0)
+opts.Add('check', 'Set to 1 if test programs shall be ran', 0)
+opts.Add('with_concepts', 'Set to 1 if test programs shall be compiled with conceptg++ (must be in ~/bin)', 0)
+opts.Add('full_warnings', 'Set to 1 to turn on all possible warnings and transform into errors (currently only gcc)', 0)
+
+
+opts.Add('add_ccflag', 'Add extra CC compiler flags that aren\'t automatically set', '')
+opts.Add('add_cxxflag', 'Add extra CXX compiler flags that aren\'t automatically set', '')
+opts.Add('add_optflag', 'Add extra optimization flags that aren\'t automatically set', '')
+
+############
+# BLAS flags
+############
+
+opts.Add('with_blas', 'Link with BLAS (have to given each time :-( until we can save configurations)', 0)
+opts.Add('blas_path', 'Add path of blas library', '')
+opts.Add('blas_lib', 'Add name of blas library', '')
+opts.Add('blas_ldflags', 'Add LDFLAGS for blas library', '')
+
+############
+# UMFPACK flags
+############
+opts.Add('with_umfpack', 'Use UMFPACK. Needs umfpack_path and amd_path.', 0)
+opts.Add('umfpack_path', 'Add path of UMFPACK library', '')
+opts.Add('amd_path', 'Add path of AMD library', '')
+opts.Add('ufconfig_path', 'Add path of AMD library', '')
+
+
+Help(opts.GenerateHelpText(env))
+
+#env.Append(CCFLAGS = ARGUMENTS.get('debug', ''))
+# Hack, for test only
+# env.Replace(CXX = '/san/intel/cc/9.0/bin/icc')
+
+# To test with conceptg++
+# requires
+if int(ARGUMENTS.get('with_concepts', 0)):
+    env['CXX'] = '~/bin/conceptg++'
+    env.Append(CCFLAGS = macro_flag(env, 'CONCEPTS_WITHOUT_OVERLOADED_REQUIREMENTS'))
+
+
+# Full warnings and treat as errors
+if int(ARGUMENTS.get('full_warnings', 0)):
+    if env['CXX'] == 'g++':
+        env.Append(CXXFLAGS = '-Werror -Wall -Wreorder -pedantic -Wpointer-arith -Wcast-qual -Wcast-align -Wwrite-strings')
+full_warnings = ARGUMENTS.get('full_warnings', 0)
+
+
+# add user-defined CC flags
+add_ccflag = ARGUMENTS.get('add_ccflag', '')
+if add_ccflag:
+     env.Append(CCFLAGS = add_ccflag)
+
+# add user-defined CXX flags
+add_cxxflag = ARGUMENTS.get('add_cxxflag', '')
+if add_cxxflag:
+     env.Append(CXXFLAGS = add_cxxflag)
+
+add_platform_flags(env)
+
+#create alternative environments
+debug_env = env.Copy()
+opt_env = env.Copy()
+high_opt_env = env.Copy()
+
+# add debugging flags to appropriate environment
+debug = ARGUMENTS.get('debug', 0)
+if int(debug):
+    add_debug(env)
+add_debug(debug_env)
+
+
+# add optimization flags to appropriate environment
+opt = ARGUMENTS.get('opt', 0)
+
+# Add optimization flags and then copy C flags into C++ flags
+add_opt(debug_env, 0)
+add_opt(env, int(opt))
+add_opt(opt_env, 1)
+add_opt(high_opt_env, 2)
+
+# add user-defined optimization flags
+add_optflag = ARGUMENTS.get('add_optflag', '')
+if add_optflag:
+     add_user_opt(env, add_optflag)
+     add_user_opt(opt_env, add_optflag)
+     add_user_opt(high_opt_env, add_optflag)
+
+# whether test programs should be ran
+check = ARGUMENTS.get('check', 0)
+
+#####################
+# Sub-directories
+#####################
+
+Export('env debug_env opt_env high_opt_env check full_warnings')
+
+SConscript(['libs/numeric/mtl/build/SConscript',
+            'libs/numeric/mtl/test/SConscript',
+            'libs/numeric/mtl/test_with_optimization/SConscript',
+            'libs/numeric/mtl/examples/SConscript',
+            'libs/numeric/mtl/experimental/SConscript',
+            'libs/numeric/mtl/timing/SConscript',
+            'libs/numeric/linear_algebra/test/SConscript',
+            'libs/numeric/itl/test/SConscript'])
+

AMDiS/lib/mtl4/boost/numeric/itl/iteration/basic_iteration.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef ITL_BASIC_ITERATION_INCLUDE
+#define ITL_BASIC_ITERATION_INCLUDE
+
+#include <iostream>
+#include <complex>
+#include <string>
+
+namespace itl {
+
+
+template <class Real>
+class basic_iteration
+{
+  public:
+    typedef basic_iteration self;
+    typedef Real            real;
+
+    template <class Vector>
+    basic_iteration(const Vector& r0, int max_iter_, Real t, Real a = Real(0))
+      : error(0), i(0), norm_r0(std::abs(two_norm(r0))),
+	max_iter(max_iter_), rtol_(t), atol_(a), is_finished(false), my_quite(false), my_suppress(false) { }
+
+    basic_iteration(Real nb, int max_iter_, Real t, Real a = Real(0))
+      : error(0), i(0), norm_r0(nb), max_iter(max_iter_), rtol_(t), atol_(a), is_finished(false),
+	my_quite(false), my_suppress(false) {}
+
+    virtual ~basic_iteration() {}
+
+    bool check_max()
+    {
+	if (i >= max_iter)
+	    error= 1, is_finished= true, err_msg= "Too many iterations.";
+	return is_finished;
+    }
+
+    template <class Vector>
+    bool finished(const Vector& r)
+    {
+	if (converged(two_norm(r)))
+	    return is_finished= true;
+	return check_max();
+    }
+
+    bool finished(const Real& r)
+    {
+	if (converged(r))
+	    return is_finished= true;
+	return check_max();
+    }
+
+    template <typename T>
+    bool finished(const std::complex<T>& r)
+    {
+	if (converged(std::abs(r)))
+	    return is_finished= true;
+	return check_max();
+    }
+
+    bool finished() const { return is_finished; }
+
+    template <class T>
+    int terminate(const T& r) { finished(r); return error; }
+
+    bool converged(const Real& r) { resid_= r; return converged(); }
+
+    bool converged() const
+    {
+	if (norm_r0 == 0)
+	    return resid_ <= atol_;  // ignore relative tolerance if |r0| is zero
+	return resid_ <= rtol_ * norm_r0 || resid_ <= atol_;
+    }
+
+    self& operator++() { ++i; return *this; }
+
+    self& operator+=(int n) { i+= n; return *this; }
+
+    bool first() const { return i <= 1; }
+
+    virtual operator int() const { return error; }
+
+    virtual int error_code() const { return error; }
+
+    bool is_converged() const { return is_finished && error == 0; }
+
+    int iterations() const { return i; }
+
+    int max_iterations() const { return max_iter; }
+
+    void set_max_iterations(int m) { max_iter= m; }
+
+    Real resid() const { return resid_; }
+
+    Real relresid() const { return resid_ / norm_r0; }
+
+    Real normb() const { return norm_r0; }
+
+    Real tol() const { return rtol_; }
+    Real atol() const { return atol_; }
+
+    int fail(int err_code) { error = err_code; return error_code(); }
+
+    int fail(int err_code, const std::string& msg)
+    { error = err_code; err_msg = msg; return error_code(); }
+
+    void set(Real v) { norm_r0 = v; }
+
+    void set_quite(bool q) { my_quite= q; }
+
+    bool is_quite() const { return my_quite; }
+
+    void suppress_resume(bool s) { my_suppress= s; }
+
+    bool resume_suppressed() const { return my_suppress; }
+
+    void update_progress(const basic_iteration& that)
+    {
+	i= that.i;
+	resid_= that.resid_;
+	if (that.error > 1) { // copy error except too many iterations
+	    error= that.error;
+	    err_msg= that.err_msg;
+	    is_finished= true;
+	} else
+	    finished(resid_);
+    }
+
+  protected:
+    int          error, i;
+    Real         norm_r0;
+    int          max_iter;
+    Real         rtol_, atol_, resid_;
+    std::string  err_msg;
+    bool         is_finished, my_quite, my_suppress;
+};
+
+
+} // namespace itl
+
+#endif // ITL_BASIC_ITERATION_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/iteration/cyclic_iteration.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library2
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef ITL_CYCLIC_ITERATION_INCLUDE
+#define ITL_CYCLIC_ITERATION_INCLUDE
+
+#include <iostream>
+#include <boost/numeric/itl/iteration/basic_iteration.hpp>
+
+namespace itl {
+
+  /// Class for iteration control that cyclically prints residual
+  template <class Real, class OStream = std::ostream>
+  class cyclic_iteration : public basic_iteration<Real>
+  {
+      typedef basic_iteration<Real> super;
+      typedef cyclic_iteration self;
+
+      void print_resid()
+      {
+	  if (!this->my_quite && this->i % cycle == 0)
+	      if (multi_print || this->i != last_print) { // Avoid multiple print-outs in same iteration
+		  out << "iteration " << this->i << ": resid " << this->resid()
+		      // << " / " << this->norm_r0 << " = " << this->resid() / this->norm_r0 << " (rel. error)"
+		      << std::endl;
+		  last_print= this->i;
+	      }
+      }
+
+    public:
+
+      template <class Vector>
+      cyclic_iteration(const Vector& r0, int max_iter_, Real tol_, Real atol_ = Real(0), int cycle_ = 100,
+		       OStream& out = std::cout)
+	: super(r0, max_iter_, tol_, atol_), cycle(cycle_), last_print(-1), multi_print(false), out(out)
+      {}
+
+      cyclic_iteration(Real r0, int max_iter_, Real tol_, Real atol_ = Real(0), int cycle_ = 100,
+		       OStream& out = std::cout)
+	: super(r0, max_iter_, tol_, atol_), cycle(cycle_), last_print(-1), multi_print(false), out(out)
+      {}
+
+
+      bool finished() { return super::finished(); }
+
+      template <typename T>
+      bool finished(const T& r)
+      {
+	  bool ret= super::finished(r);
+	  print_resid();
+	  return ret;
+      }
+
+      inline self& operator++() { ++this->i; return *this; }
+
+      inline self& operator+=(int n) { this->i+= n; return *this; }
+
+      operator int() const { return error_code(); }
+
+      /// Whether the residual is printed multiple times in iteration
+      bool is_multi_print() const { return multi_print; }
+
+      /// Set whether the residual is printed multiple times in iteration
+      void set_multi_print(bool m) { multi_print= m; }
+
+      int error_code() const
+      {
+	  if (!this->my_suppress)
+	      out << "finished! error code = " << this->error << '\n'
+		  << this->iterations() << " iterations\n"
+		  << this->resid() << " is actual final residual. \n"
+		  << this->relresid() << " is actual relative tolerance achieved. \n"
+		  << "Relative tol: " << this->rtol_ << "  Absolute tol: " << this->atol_ << '\n'
+		  << "Convergence:  " << pow(this->relresid(), 1.0 / double(this->iterations())) << std::endl;
+	  return this->error;
+      }
+    protected:
+      int        cycle, last_print;
+      bool       multi_print;
+      OStream&   out;
+  };
+
+
+
+} // namespace itl
+
+#endif // ITL_CYCLIC_ITERATION_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/iteration/noisy_iteration.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef ITL_NOISY_ITERATION_INCLUDE
+#define ITL_NOISY_ITERATION_INCLUDE
+
+#include <iostream>
+#include <boost/numeric/itl/iteration/cyclic_iteration.hpp>
+
+namespace itl {
+
+  template <class Real, class OStream = std::ostream>
+  class noisy_iteration : public cyclic_iteration<Real, OStream>
+  {
+    public:
+      template <class Vector>
+      noisy_iteration(const Vector& r0, int max_iter_, Real tol_, Real atol_ = Real(0),
+		      OStream& out = std::cout)
+	: cyclic_iteration<Real, OStream>(r0, max_iter_, tol_, atol_, 1, out)
+      {}
+  };
+
+} // namespace itl
+
+#endif // ITL_NOISY_ITERATION_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/itl.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef ITL_ITL_INCLUDE
+#define ITL_ITL_INCLUDE
+
+
+#include <boost/numeric/itl/iteration/basic_iteration.hpp>
+#include <boost/numeric/itl/iteration/cyclic_iteration.hpp>
+#include <boost/numeric/itl/iteration/noisy_iteration.hpp>
+
+#include <boost/numeric/itl/krylov/cg.hpp>
+#include <boost/numeric/itl/krylov/cgs.hpp>
+#include <boost/numeric/itl/krylov/bicg.hpp>
+#include <boost/numeric/itl/krylov/bicgstab.hpp>
+#include <boost/numeric/itl/krylov/bicgstab_2.hpp>
+#include <boost/numeric/itl/krylov/bicgstab_ell.hpp>
+#include <boost/numeric/itl/krylov/fsm.hpp>
+#include <boost/numeric/itl/krylov/idr_s.hpp>
+#include <boost/numeric/itl/krylov/gmres.hpp>
+#include <boost/numeric/itl/krylov/tfqmr.hpp>
+#include <boost/numeric/itl/krylov/qmr.hpp>
+
+#include <boost/numeric/itl/krylov/repeating_solver.hpp>
+
+#include <boost/numeric/itl/minimization/quasi_newton.hpp>
+
+#include <boost/numeric/itl/pc/identity.hpp>
+#include <boost/numeric/itl/pc/is_identity.hpp>
+#include <boost/numeric/itl/pc/diagonal.hpp>
+#include <boost/numeric/itl/pc/ilu.hpp>
+#include <boost/numeric/itl/pc/ilu_0.hpp>
+#include <boost/numeric/itl/pc/ilut.hpp>
+#include <boost/numeric/itl/pc/ic_0.hpp>
+
+#include <boost/numeric/itl/pc/imf_preconditioner.hpp>
+#include <boost/numeric/itl/pc/imf_algorithms.hpp>
+
+#include <boost/numeric/itl/pc/sub_matrix_pc.hpp>
+#include <boost/numeric/itl/pc/concat.hpp>
+
+#include <boost/numeric/itl/smoother/gauss_seidel.hpp>
+
+#include <boost/numeric/itl/stepper/armijo.hpp>
+#include <boost/numeric/itl/stepper/wolf.hpp>
+
+#include <boost/numeric/itl/updater/bfgs.hpp>
+#include <boost/numeric/itl/updater/broyden.hpp>
+#include <boost/numeric/itl/updater/dfp.hpp>
+#include <boost/numeric/itl/updater/psb.hpp>
+#include <boost/numeric/itl/updater/sr1.hpp>
+
+#endif // ITL_ITL_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/itl_fwd.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef ITL_ITL_FWD_INCLUDE
+#define ITL_ITL_FWD_INCLUDE
+
+#include <boost/numeric/mtl/concept/collection.hpp>
+
+namespace itl {
+
+    template <class Real>                class basic_iteration;
+    template <class Real, class OStream> class cyclic_iteration;
+    template <class Real, class OStream> class noisy_iteration;
+
+    template <typename Solver, typename VectorIn, bool trans> class solver_proxy;
+
+    namespace pc {
+
+	template <typename PC, typename Vector, bool> struct solver;
+	template <typename Matrix, typename Value> class identity;
+	// template <typename Matrix, typename Value, typename Vector> Vector solve(const identity<Matrix>&, const Vector& x);
+	// template <typename Matrix, typename Vector> Vector adjoint_solve(const identity<Matrix>&, const Vector& x);
+
+	template <typename Matrix, typename Value> class diagonal;
+	// template <typename Matrix, typename Vector> Vector solve(const diagonal<Matrix>& P, const Vector& x);
+	// template <typename Matrix, typename Vector> Vector adjoint_solve(const diagonal<Matrix>& P, const Vector& x);
+
+	template <typename Matrix, typename Factorizer, typename Value> class ilu;
+	template <typename Matrix, typename Value> class ilu_0; // Maybe we should declare the default here???
+	template <typename Matrix, typename Value> class ilut; // Maybe we should declare the default here???
+	// template <typename Matrix, typename Vector> Vector solve(const ilu_0<Matrix>& P, const Vector& x);
+	// template <typename Matrix, typename Vector> Vector adjoint_solve(const ilu_0<Matrix>& P, const Vector& x);
+
+	template <typename Matrix, typename Value> class ic_0; // Maybe we should declare the default here???
+	// template <typename Matrix, typename Vector> Vector solve(const ic_0<Matrix>& P, const Vector& x);
+	// template <typename Matrix, typename Value, typename Vector> Vector adjoint_solve(const ic_0<Matrix, Value>& P, const Vector& x);
+
+    } //  namespace pc
+
+    template < typename LinearOperator, typename HilbertSpaceX, typename HilbertSpaceB,
+	       typename Preconditioner, typename Iteration >
+    int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
+	   const Preconditioner& M, Iteration& iter);
+
+    template < typename LinearOperator, typename Preconditioner= pc::identity<LinearOperator, double>,
+	       typename RightPreconditioner= pc::identity<LinearOperator, double> >
+    class cg_solver;
+
+
+    template < typename LinearOperator, typename Vector,
+	       typename Preconditioner, typename Iteration >
+    int bicg(const LinearOperator &A, Vector &x, const Vector &b,
+	     const Preconditioner &M, Iteration& iter);
+
+    template < class LinearOperator, class HilbertSpaceX, class HilbertSpaceB,
+	       class Preconditioner, class Iteration >
+    int bicgstab(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
+		 const Preconditioner& M, Iteration& iter);
+
+    template < class LinearOperator, class HilbertSpaceX, class HilbertSpaceB,
+	       class Preconditioner, class Iteration >
+    int bicgstab_2(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
+		   const Preconditioner& M, Iteration& iter);
+
+    template < typename LinearOperator, typename Vector,
+	       typename LeftPreconditioner, typename RightPreconditioner,
+	       typename Iteration >
+    int bicgstab_ell(const LinearOperator &A, Vector &x, const Vector &b,
+		     const LeftPreconditioner &L, const RightPreconditioner &R,
+		     Iteration& iter, size_t l);
+
+    template < typename LinearOperator, typename Preconditioner= pc::identity<LinearOperator, double>,
+	       typename RightPreconditioner= pc::identity<LinearOperator, double> >
+    class bicgstab_ell_solver;
+
+    template <typename Solver, unsigned N, bool Stored= false>
+    class repeating_solver;
+
+} // namespace itl
+
+#endif // ITL_ITL_FWD_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/krylov/bicg.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef ITL_BICG_INCLUDE
+#define ITL_BICG_INCLUDE
+
+#include <complex>
+#include <boost/numeric/itl/itl_fwd.hpp>
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/operation/conj.hpp>
+#include <boost/numeric/mtl/operation/resource.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+namespace itl {
+
+/// Bi-Conjugate Gradient
+template < typename LinearOperator, typename Vector,
+	   typename Preconditioner, typename Iteration >
+int bicg(const LinearOperator &A, Vector &x, const Vector &b,
+	 const Preconditioner &M, Iteration& iter)
+{
+    mtl::vampir_trace<7003> tracer;
+    using mtl::conj;
+    typedef typename mtl::Collection<Vector>::value_type Scalar;
+    Scalar     rho_1(0), rho_2(0), alpha(0), beta(0);
+    Vector     r(b - A * x), z(resource(x)), p(resource(x)), q(resource(x)),
+ 	       r_tilde(r), z_tilde(resource(x)), p_tilde(resource(x)), q_tilde(resource(x));
+
+    while ( ! iter.finished(r)) {
+	++iter;
+	z= solve(M, r);
+	z_tilde= adjoint_solve(M, r_tilde);
+	rho_1= dot(z_tilde, z);
+
+	if (rho_1 == 0.) return iter.fail(2, "bicg breakdown");
+	if (iter.first()) {
+	    p= z;
+	    p_tilde= z_tilde;
+	} else {
+	    beta= rho_1 / rho_2;
+	    p= z + beta * p;
+	    p_tilde= z_tilde + conj(beta) * p_tilde;
+	}
+
+	q= A * p;
+	q_tilde= adjoint(A) * p_tilde;
+	alpha= rho_1 / dot(p_tilde, q);
+
+	x+= alpha * p;
+	r-= alpha * q;
+	r_tilde-= conj(alpha) * q_tilde;
+
+	rho_2= rho_1;
+    }
+    return iter;
+}
+
+/// Solver class for BiCG method; right preconditioner ignored (prints warning if not identity)
+template < typename LinearOperator, typename Preconditioner= pc::identity<LinearOperator>,
+	   typename RightPreconditioner= pc::identity<LinearOperator> >
+class bicg_solver
+{
+  public:
+    /// Construct solver from a linear operator; generate (left) preconditioner from it
+    explicit bicg_solver(const LinearOperator& A) : A(A), L(A)
+    {
+	if (!pc::static_is_identity<RightPreconditioner>::value)
+	    std::cerr << "Right Preconditioner ignored!" << std::endl;
+    }
+
+    /// Construct solver from a linear operator and (left) preconditioner
+    bicg_solver(const LinearOperator& A, const Preconditioner& L) : A(A), L(L)
+    {
+	if (!pc::static_is_identity<RightPreconditioner>::value)
+	    std::cerr << "Right Preconditioner ignored!" << std::endl;
+    }
+
+    /// Solve linear system approximately as specified by \p iter
+    template < typename HilbertSpaceB, typename HilbertSpaceX, typename Iteration >
+    int solve(const HilbertSpaceB& b, HilbertSpaceX& x, Iteration& iter) const
+    {
+	return bicg(A, x, b, L, iter);
+    }
+
+    /// Perform one bicg iteration on linear system
+    template < typename HilbertSpaceB, typename HilbertSpaceX >
+    int solve(const HilbertSpaceB& b, HilbertSpaceX& x) const
+    {
+	itl::basic_iteration<double> iter(x, 1, 0, 0);
+	return bicg(A, x, b, L, iter);
+    }
+
+  private:
+    const LinearOperator& A;
+    Preconditioner        L;
+};
+
+} // namespace itl
+
+#endif // ITL_BICG_INCLUDE
+
+
+
+
+
+

AMDiS/lib/mtl4/boost/numeric/itl/krylov/bicgstab.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef ITL_BICGSTAB_INCLUDE
+#define ITL_BICGSTAB_INCLUDE
+
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/itl/utility/exception.hpp>
+#include <boost/numeric/mtl/operation/resource.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+namespace itl {
+
+///  Bi-Conjugate Gradient Stabilized
+template < class LinearOperator, class HilbertSpaceX, class HilbertSpaceB,
+	   class Preconditioner, class Iteration >
+int bicgstab(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
+	     const Preconditioner& M, Iteration& iter)
+{
+  typedef typename mtl::Collection<HilbertSpaceX>::value_type Scalar;
+  typedef HilbertSpaceX                                       Vector;
+  mtl::vampir_trace<7004> tracer;
+
+  Scalar     rho_1(0), rho_2(0), alpha(0), beta(0), gamma, omega(0);
+  Vector     p(resource(x)), phat(resource(x)), s(resource(x)), shat(resource(x)),
+             t(resource(x)), v(resource(x)), r(resource(x)), rtilde(resource(x));
+
+  r = b - A * x;
+  rtilde = r;
+
+  while (! iter.finished(r)) {
+    ++iter;
+    rho_1 = dot(rtilde, r);
+    MTL_THROW_IF(rho_1 == 0.0, unexpected_orthogonality());
+
+    if (iter.first())
+      p = r;
+    else {
+      MTL_THROW_IF(omega == 0.0, unexpected_orthogonality());
+      beta = (rho_1 / rho_2) * (alpha / omega);
+      p = r + beta * (p - omega * v);
+    }
+    phat = solve(M, p);
+    v = A * phat;
+
+    gamma = dot(rtilde, v);
+    MTL_THROW_IF(gamma == 0.0, unexpected_orthogonality());
+
+    alpha = rho_1 / gamma;
+    s = r - alpha * v;
+
+    if (iter.finished(s)) {
+      x += alpha * phat;
+      break;
+    }
+    shat = solve(M, s);
+    t = A * shat;
+    omega = dot(t, s) / dot(t, t);
+
+    x += omega * shat + alpha * phat;
+    r = s - omega * t;
+
+    rho_2 = rho_1;
+  }
+  return iter;
+}
+
+/// Solver class for BiCGStab method; right preconditioner ignored (prints warning if not identity)
+template < typename LinearOperator, typename Preconditioner= pc::identity<LinearOperator>,
+	   typename RightPreconditioner= pc::identity<LinearOperator> >
+class bicgstab_solver
+{
+  public:
+    /// Construct solver from a linear operator; generate (left) preconditioner from it
+    explicit bicgstab_solver(const LinearOperator& A) : A(A), L(A)
+    {
+	if (!pc::static_is_identity<RightPreconditioner>::value)
+	    std::cerr << "Right Preconditioner ignored!" << std::endl;
+    }
+
+    /// Construct solver from a linear operator and (left) preconditioner
+    bicgstab_solver(const LinearOperator& A, const Preconditioner& L) : A(A), L(L)
+    {
+	if (!pc::static_is_identity<RightPreconditioner>::value)
+	    std::cerr << "Right Preconditioner ignored!" << std::endl;
+    }
+
+    /// Solve linear system approximately as specified by \p iter
+    template < typename HilbertSpaceB, typename HilbertSpaceX, typename Iteration >
+    int solve(const HilbertSpaceB& b, HilbertSpaceX& x, Iteration& iter) const
+    {
+	return bicgstab(A, x, b, L, iter);
+    }
+
+    /// Perform one BiCGStab iteration on linear system
+    template < typename HilbertSpaceB, typename HilbertSpaceX >
+    int solve(const HilbertSpaceB& b, HilbertSpaceX& x) const
+    {
+	itl::basic_iteration<double> iter(x, 1, 0, 0);
+	return solve(b, x, iter);
+    }
+
+  private:
+    const LinearOperator& A;
+    Preconditioner        L;
+};
+
+} // namespace itl
+
+#endif // ITL_BICGSTAB_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/krylov/bicgstab_2.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef ITL_BICGSTAB_2_INCLUDE
+#define ITL_BICGSTAB_2_INCLUDE
+
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/linear_algebra/identity.hpp>
+#include <boost/numeric/mtl/operation/resource.hpp>
+#include <boost/numeric/mtl/operation/size.hpp>
+#include <boost/numeric/itl/pc/identity.hpp>
+#include <boost/numeric/mtl/vector/dense_vector.hpp>
+#include <boost/numeric/mtl/matrix/strict_upper.hpp>
+#include <boost/numeric/mtl/matrix/dense2D.hpp>
+#include <boost/numeric/mtl/utility/irange.hpp>
+#include <boost/numeric/mtl/operation/orth.hpp>
+#include <boost/numeric/mtl/operation/lazy.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+namespace itl {
+
+/// Bi-Conjugate Gradient Stabilized(2)
+template < typename LinearOperator, typename Vector,
+	   typename Preconditioner, typename Iteration >
+int bicgstab_2(const LinearOperator &A, Vector &x, const Vector &b,
+	       const Preconditioner &L, Iteration& iter)
+{
+    mtl::vampir_trace<7005> tracer;
+    using mtl::size; using mtl::irange; using mtl::imax; using mtl::matrix::strict_upper; using mtl::lazy;
+    typedef typename mtl::Collection<Vector>::value_type Scalar;
+    typedef typename mtl::Collection<Vector>::size_type  Size;
+
+    if (size(b) == 0) throw mtl::logic_error("empty rhs vector");
+
+    const size_t                l= 2;
+    const Scalar                zero= math::zero(Scalar()), one= math::one(Scalar());
+    Vector                      x0(resource(x)), y(resource(x));
+    mtl::vector::dense_vector<Vector>   r_hat(l+1,Vector(resource(x))), u_hat(l+1,Vector(resource(x)));
+
+    // shift problem
+    x0= zero;
+    r_hat[0]= b;
+    if (two_norm(x) != zero) {
+	r_hat[0]-= A * x;
+	x0= x;
+	x= zero;
+    }
+
+    Vector  r0_tilde(r_hat[0]/two_norm(r_hat[0]));
+    y= solve(L, r_hat[0]);
+    r_hat[0]= y;
+    u_hat[0]= zero;
+
+    Scalar                      rho_0(one), rho_1(zero), alpha(zero), Gamma(zero), beta(zero), omega(one);
+    mtl::matrix::dense2D<Scalar>        tau(l+1, l+1);
+    mtl::vector::dense_vector<Scalar>   sigma(l+1), gamma(l+1), gamma_a(l+1), gamma_aa(l+1);
+
+    while (! iter.finished(r_hat[0])) {
+	++iter;
+	rho_0= -omega * rho_0;
+
+	for (Size j= 0; j < 2; ++j) {
+	    rho_1= dot(r0_tilde, r_hat[j]);
+	    beta= alpha * rho_1/rho_0; rho_0= rho_1;
+
+	    for (Size i= 0; i <= j; ++i)
+		u_hat[i]= r_hat[i] - beta * u_hat[i];
+
+	    y= A * u_hat[j];
+	    u_hat[j+1]= solve(L, y);
+	    Gamma= dot(r0_tilde, u_hat[j+1]);
+	    alpha= rho_0 / Gamma;
+
+	    for (Size i= 0; i <= j; ++i)
+		r_hat[i]-= alpha * u_hat[i+1];
+
+	    if (iter.finished(r_hat[j])) {
+		x+= x0;
+		return iter;
+	    }
+
+	    y= A * r_hat[j];
+	    r_hat[j+1]= solve(L, y);
+	    x+= alpha * u_hat[0];
+	}
+
+	// mod GS (MR part)
+	irange  i1m(1, imax);
+	mtl::vector::dense_vector<Vector>   r_hat_tail(r_hat[i1m]);
+	tau[i1m][i1m]= orthogonalize_factors(r_hat_tail);
+	for (Size j= 1; j <= l; ++j)
+	    gamma_a[j]= dot(r_hat[j], r_hat[0]) / tau[j][j];
+
+	gamma[l]= gamma_a[l]; omega= gamma[l];
+	if (omega == zero) return iter.fail(3, "bicg breakdown #2");
+
+	// is this something like a tri-solve?
+	for (Size j= l-1; j > 0; --j) {
+	    Scalar sum= zero;
+	    for (Size i=j+1;i<=l;++i)
+		sum += tau[j][i] * gamma[i];
+	    gamma[j] = gamma_a[j] - sum;
+	}
+
+	gamma_aa[irange(1, l)]= strict_upper(tau[irange(1, l)][irange(1, l)]) * gamma[irange(2, l+1)] + gamma[irange(2, l+1)];
+
+	x+= gamma[1] * r_hat[0];
+	r_hat[0]-= gamma_a[l] * r_hat[l];
+	u_hat[0]-= gamma[l] * u_hat[l];
+	for (Size j=1; j < l; ++j) {
+	    u_hat[0] -= gamma[j] * u_hat[j];
+	    x+= gamma_aa[j] * r_hat[j];
+	    r_hat[0] -= gamma_a[j] * r_hat[j];
+	}
+    }
+    x+= x0; // convert to real solution and undo shift
+    return iter;
+}
+
+
+
+/// Solver class for BiCGStab(2) method; right preconditioner ignored (prints warning if not identity)
+template < typename LinearOperator, typename Preconditioner= pc::identity<LinearOperator>,
+	   typename RightPreconditioner= pc::identity<LinearOperator> >
+class bicgstab_2_solver
+{
+  public:
+    /// Construct solver from a linear operator; generate (left) preconditioner from it
+    explicit bicgstab_2_solver(const LinearOperator& A) : A(A), L(A)
+    {
+	if (!pc::static_is_identity<RightPreconditioner>::value)
+	    std::cerr << "Right Preconditioner ignored!" << std::endl;
+    }
+
+    /// Construct solver from a linear operator and (left) preconditioner
+    bicgstab_2_solver(const LinearOperator& A, const Preconditioner& L) : A(A), L(L)
+    {
+	if (!pc::static_is_identity<RightPreconditioner>::value)
+	    std::cerr << "Right Preconditioner ignored!" << std::endl;
+    }
+
+    /// Solve linear system approximately as specified by \p iter
+    template < typename HilbertSpaceB, typename HilbertSpaceX, typename Iteration >
+    int solve(const HilbertSpaceB& b, HilbertSpaceX& x, Iteration& iter) const
+    {
+	return bicgstab_2(A, x, b, L, iter);
+    }
+
+    /// Perform one BiCGStab(2) iteration on linear system
+    template < typename HilbertSpaceB, typename HilbertSpaceX >
+    int solve(const HilbertSpaceB& b, HilbertSpaceX& x) const
+    {
+	itl::basic_iteration<double> iter(x, 1, 0, 0);
+	return bicgstab_2(A, x, b, L, iter);
+    }
+
+  private:
+    const LinearOperator& A;
+    Preconditioner        L;
+};
+
+
+
+} // namespace itl
+
+#endif // ITL_BICGSTAB_2_INCLUDE
+
+
+
+
+
+

AMDiS/lib/mtl4/boost/numeric/itl/krylov/bicgstab_ell.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+// Written by Jan Bos
+// Edited  by Peter Gottschling
+
+#ifndef ITL_BICGSTAB_ELL_INCLUDE
+#define ITL_BICGSTAB_ELL_INCLUDE
+
+#include <boost/numeric/itl/itl_fwd.hpp>
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/vector/dense_vector.hpp>
+#include <boost/numeric/mtl/matrix/strict_upper.hpp>
+#include <boost/numeric/mtl/matrix/dense2D.hpp>
+#include <boost/numeric/mtl/utility/irange.hpp>
+#include <boost/numeric/mtl/operation/orth.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/linear_algebra/identity.hpp>
+#include <boost/numeric/mtl/operation/resource.hpp>
+#include <boost/numeric/mtl/operation/size.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+namespace itl {
+
+/// Bi-Conjugate Gradient Stabilized(ell)
+template < typename LinearOperator, typename Vector,
+	   typename LeftPreconditioner, typename RightPreconditioner,
+	   typename Iteration >
+int bicgstab_ell(const LinearOperator &A, Vector &x, const Vector &b,
+		 const LeftPreconditioner &L, const RightPreconditioner &R,
+		 Iteration& iter, size_t l)
+{
+    mtl::vampir_trace<7006> tracer;
+    using mtl::size; using mtl::irange; using mtl::imax; using mtl::matrix::strict_upper;
+    typedef typename mtl::Collection<Vector>::value_type Scalar;
+    typedef typename mtl::Collection<Vector>::size_type  Size;
+
+    if (size(b) == 0) throw mtl::logic_error("empty rhs vector");
+
+    const Scalar                zero= math::zero(Scalar()), one= math::one(Scalar());
+    Vector                      x0(resource(x)), y(resource(x));
+    mtl::vector::dense_vector<Vector>   r_hat(l+1,Vector(resource(x))), u_hat(l+1,Vector(resource(x)));
+
+    // shift problem
+    x0= zero;
+    r_hat[0]= b;
+    if (two_norm(x) != zero) {
+	r_hat[0]-= A * x;
+	x0= x;
+	x= zero;
+    }
+
+    Vector  r0_tilde(r_hat[0]/two_norm(r_hat[0]));
+    y= solve(L, r_hat[0]);
+    r_hat[0]= y;
+    u_hat[0]= zero;
+
+    Scalar                      rho_0(one), rho_1(zero), alpha(zero), Gamma(zero), beta(zero), omega(one);
+    mtl::matrix::dense2D<Scalar>        tau(l+1, l+1);
+    mtl::vector::dense_vector<Scalar>   sigma(l+1), gamma(l+1), gamma_a(l+1), gamma_aa(l+1);
+
+    while (! iter.finished(r_hat[0])) {
+	++iter;
+	rho_0= -omega * rho_0;
+
+	for (Size j= 0; j < l; ++j) {
+	    rho_1= dot(r0_tilde, r_hat[j]);
+	    beta= alpha * rho_1/rho_0; rho_0= rho_1;
+
+	    for (Size i= 0; i <= j; ++i)
+		u_hat[i]= r_hat[i] - beta * u_hat[i];
+
+	    y= A * Vector(solve(R, u_hat[j]));
+	    u_hat[j+1]= solve(L, y);
+	    Gamma= dot(r0_tilde, u_hat[j+1]);
+	    alpha= rho_0 / Gamma;
+
+	    for (Size i= 0; i <= j; ++i)
+		r_hat[i]-= alpha * u_hat[i+1];
+
+	    if (iter.finished(r_hat[j])) {
+		x= solve(R, x);
+		x+= x0;
+		return iter;
+	    }
+
+	    r_hat[j+1]= solve(R, r_hat[j]);
+	    y= A * r_hat[j+1];
+	    r_hat[j+1]= solve(L, y);
+	    x+= alpha * u_hat[0];
+	}
+
+	// mod GS (MR part)
+	irange  i1m(1, imax);
+	mtl::vector::dense_vector<Vector>   r_hat_tail(r_hat[i1m]);
+	tau[i1m][i1m]= orthogonalize_factors(r_hat_tail);
+	for (Size j= 1; j <= l; ++j)
+	    gamma_a[j]= dot(r_hat[j], r_hat[0]) / tau[j][j];
+
+	gamma[l]= gamma_a[l]; omega= gamma[l];
+	if (omega == zero) return iter.fail(3, "bicg breakdown #2");
+
+	// is this something like a tri-solve?
+	for (Size j= l-1; j > 0; --j) {
+	    Scalar sum= zero;
+	    for (Size i=j+1;i<=l;++i)
+		sum += tau[j][i] * gamma[i];
+	    gamma[j] = gamma_a[j] - sum;
+	}
+
+	gamma_aa[irange(1, l)]= strict_upper(tau[irange(1, l)][irange(1, l)]) * gamma[irange(2, l+1)] + gamma[irange(2, l+1)];
+
+	x+= gamma[1] * r_hat[0];
+	r_hat[0]-= gamma_a[l] * r_hat[l];
+	u_hat[0]-= gamma[l] * u_hat[l];
+	for (Size j=1; j < l; ++j) {
+	    u_hat[0] -= gamma[j] * u_hat[j];
+	    x+= gamma_aa[j] * r_hat[j];
+	    r_hat[0] -= gamma_a[j] * r_hat[j];
+	}
+    }
+    x= solve(R, x); x+= x0; // convert to real solution and undo shift
+    return iter;
+}
+
+
+/// Solver class for BiCGStab(ell) method; right preconditioner ignored (prints warning if not identity)
+template < typename LinearOperator, typename Preconditioner,
+	   typename RightPreconditioner>
+class bicgstab_ell_solver
+{
+  public:
+    /// Construct solver from a linear operator; generate (left) preconditioner from it
+    explicit bicgstab_ell_solver(const LinearOperator& A, size_t l= 8) : A(A), l(l), L(A), R(A) {}
+
+    /// Construct solver from a linear operator and left preconditioner
+    bicgstab_ell_solver(const LinearOperator& A, size_t l, const Preconditioner& L) : A(A), l(l), L(L), R(A) {}
+
+    /// Construct solver from a linear operator and left preconditioner
+    bicgstab_ell_solver(const LinearOperator& A, size_t l, const Preconditioner& L, const RightPreconditioner& R)
+      : A(A), l(l), L(L), R(R) {}
+
+    /// Solve linear system approximately as specified by \p iter
+    template < typename HilbertSpaceB, typename HilbertSpaceX, typename Iteration >
+    int solve(const HilbertSpaceB& b, HilbertSpaceX& x, Iteration& iter) const
+    {
+	return bicgstab_ell(A, x, b, L, R, iter, l);
+    }
+
+    /// Perform one BiCGStab(ell) iteration on linear system
+    template < typename HilbertSpaceB, typename HilbertSpaceX >
+    int solve(const HilbertSpaceB& b, HilbertSpaceX& x) const
+    {
+	itl::basic_iteration<double> iter(x, 1, 0, 0);
+	return solve(b, x, iter);
+    }
+
+  private:
+    const LinearOperator& A;
+    size_t                l;
+    Preconditioner        L;
+    RightPreconditioner   R;
+};
+
+
+} // namespace itl
+
+#endif // ITL_BICGSTAB_ELL_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/krylov/cg.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef ITL_CG_INCLUDE
+#define ITL_CG_INCLUDE
+
+#include <cmath>
+#include <cassert>
+#include <iostream>
+#include <boost/mpl/bool.hpp>
+
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/itl/itl_fwd.hpp>
+#include <boost/numeric/itl/iteration/basic_iteration.hpp>
+#include <boost/numeric/itl/pc/identity.hpp>
+#include <boost/numeric/itl/pc/is_identity.hpp>
+#include <boost/numeric/mtl/operation/dot.hpp>
+#include <boost/numeric/mtl/operation/unary_dot.hpp>
+#include <boost/numeric/mtl/operation/conj.hpp>
+#include <boost/numeric/mtl/operation/resource.hpp>
+#include <boost/numeric/mtl/operation/lazy.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+namespace itl {
+
+/// Conjugate Gradients without preconditioning
+template < typename LinearOperator, typename HilbertSpaceX, typename HilbertSpaceB,
+	   typename Iteration >
+int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
+       Iteration& iter)
+{
+    mtl::vampir_trace<7001> tracer;
+    using std::abs; using mtl::conj; using mtl::lazy;
+    typedef HilbertSpaceX Vector;
+    typedef typename mtl::Collection<HilbertSpaceX>::value_type Scalar;
+    typedef typename Iteration::real                            Real;
+
+    Scalar rho(0), rho_1(0), alpha(0), alpha_1(0);
+    Vector p(resource(x)), q(resource(x)), r(resource(x)), z(resource(x));
+
+    r = b - A*x;
+    rho = dot(r, r);
+    while (! iter.finished(Real(sqrt(abs(rho))))) {
+	++iter;
+	if (iter.first())
+	    p = r;
+	else
+	    p = r + (rho / rho_1) * p;
+
+	// q = A * p; alpha = rho / dot(p, q);
+	(lazy(q)= A * p) || (lazy(alpha_1)= lazy_dot(p, q));
+	alpha= rho / alpha_1;
+
+	x += alpha * p;
+	rho_1 = rho;
+	(lazy(r) -= alpha * q) || (lazy(rho) = lazy_unary_dot(r));
+    }
+
+    return iter;
+}
+
+/// Conjugate Gradients
+template < typename LinearOperator, typename HilbertSpaceX, typename HilbertSpaceB,
+	   typename Preconditioner, typename Iteration >
+int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
+       const Preconditioner& L, Iteration& iter)
+{
+    using pc::is_identity;
+    if (is_identity(L))
+	return cg(A, x, b, iter);
+
+    mtl::vampir_trace<7002> tracer;
+    using std::abs; using mtl::conj; using mtl::lazy;
+    typedef HilbertSpaceX Vector;
+    typedef typename mtl::Collection<HilbertSpaceX>::value_type Scalar;
+    typedef typename Iteration::real                            Real;
+
+    Scalar rho(0), rho_1(0), rr, alpha(0), alpha_1;
+    Vector p(resource(x)), q(resource(x)), r(resource(x)), z(resource(x));
+
+    r = b - A*x;
+    rr = dot(r, r);
+    while (! iter.finished(Real(sqrt(abs(rr))))) {
+	++iter;
+	(lazy(z)= solve(L, r)) || (lazy(rho)= lazy_dot(r, z));
+
+	if (iter.first())
+	    p = z;
+	else
+	    p = z + (rho / rho_1) * p;
+
+	(lazy(q)= A * p) || (lazy(alpha_1)= lazy_dot(p, q));
+	alpha= rho / alpha_1;
+
+	x += alpha * p;
+	rho_1 = rho;
+	(lazy(r) -= alpha * q) || (lazy(rr) = lazy_unary_dot(r));
+    }
+    return iter;
+}
+
+/// Conjugate Gradients with ignored right preconditioner to unify interface
+template < typename LinearOperator, typename HilbertSpaceX, typename HilbertSpaceB,
+	   typename Preconditioner, typename RightPreconditioner, typename Iteration >
+int cg(const LinearOperator& A, HilbertSpaceX& x, const HilbertSpaceB& b,
+       const Preconditioner& L, const RightPreconditioner&, Iteration& iter)
+{
+    return cg(A, x, b, L, iter);
+}
+
+/// Solver class for CG method; right preconditioner ignored (prints warning if not identity)
+template < typename LinearOperator, typename Preconditioner,
+	   typename RightPreconditioner>
+class cg_solver
+{
+  public:
+    /// Construct solver from a linear operator; generate (left) preconditioner from it
+    explicit cg_solver(const LinearOperator& A) : A(A), L(A)
+    {
+	if (!pc::static_is_identity<RightPreconditioner>::value)
+	    std::cerr << "Right Preconditioner ignored!" << std::endl;
+    }
+
+    /// Construct solver from a linear operator and (left) preconditioner
+    cg_solver(const LinearOperator& A, const Preconditioner& L) : A(A), L(L)
+    {
+	if (!pc::static_is_identity<RightPreconditioner>::value)
+	    std::cerr << "Right Preconditioner ignored!" << std::endl;
+    }
+
+    /// Solve linear system approximately as specified by \p iter
+    template < typename HilbertSpaceB, typename HilbertSpaceX, typename Iteration >
+    int solve(const HilbertSpaceB& b, HilbertSpaceX& x, Iteration& iter) const
+    {
+	return cg(A, x, b, L, iter);
+    }
+
+    /// Perform one CG iteration on linear system
+    template < typename HilbertSpaceB, typename HilbertSpaceX >
+    int solve(const HilbertSpaceB& b, HilbertSpaceX& x) const
+    {
+	itl::basic_iteration<double> iter(x, 1, 0, 0);
+	return solve(b, x, iter);
+    }
+
+  private:
+    const LinearOperator& A;
+    Preconditioner        L;
+};
+
+
+} // namespace itl
+
+#endif // ITL_CG_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/krylov/cgs.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef ITL_CGS_INCLUDE
+#define ITL_CGS_INCLUDE
+
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/operation/resource.hpp>
+#include <boost/numeric/mtl/operation/dot.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+namespace itl {
+
+/// Conjugate Gradient Squared
+template < typename LinearOperator, typename Vector,
+	   typename Preconditioner, typename Iteration >
+int cgs(const LinearOperator &A, Vector &x, const Vector &b,
+	const Preconditioner &M, Iteration& iter)
+{
+    mtl::vampir_trace<7007> tracer;
+    typedef typename mtl::Collection<Vector>::value_type Scalar;
+    Scalar     rho_1(0), rho_2(0), alpha(0), beta(0);
+    Vector     p(resource(x)), phat(resource(x)), q(resource(x)), qhat(resource(x)), vhat(resource(x)),
+	       u(resource(x)), uhat(resource(x)), r(b - A * x), rtilde= r;
+
+    while (! iter.finished(r)) {
+	++iter;
+	rho_1= dot(rtilde, r);
+
+	if (rho_1 == 0.) iter.fail(2, "cgs breakdown");
+
+	if (iter.first())
+	    p= u= r;
+	else {
+	    beta = rho_1 / rho_2;
+	    u= r + beta * q;
+	    p= u + beta * (q + beta * p);
+	}
+
+        vhat= A * Vector(solve(M, p));
+	alpha = rho_1 / dot(rtilde, vhat);
+	q= u - alpha * vhat;
+
+	u+= q;
+	uhat= solve(M, u);
+
+	x+= alpha * uhat;
+	qhat= A * uhat;
+	r-= alpha * qhat;
+
+	rho_2= rho_1;
+    }
+    return iter;
+}
+
+/// Solver class for CGS method; right preconditioner ignored (prints warning if not identity)
+template < typename LinearOperator, typename Preconditioner= pc::identity<LinearOperator>,
+	   typename RightPreconditioner= pc::identity<LinearOperator> >
+class cgs_solver
+{
+  public:
+    /// Construct solver from a linear operator; generate (left) preconditioner from it
+    explicit cgs_solver(const LinearOperator& A) : A(A), L(A)
+    {
+	if (!pc::static_is_identity<RightPreconditioner>::value)
+	    std::cerr << "Right Preconditioner ignored!" << std::endl;
+    }
+
+    /// Construct solver from a linear operator and (left) preconditioner
+    cgs_solver(const LinearOperator& A, const Preconditioner& L) : A(A), L(L)
+    {
+	if (!pc::static_is_identity<RightPreconditioner>::value)
+	    std::cerr << "Right Preconditioner ignored!" << std::endl;
+    }
+
+    /// Solve linear system approximately as specified by \p iter
+    template < typename HilbertSpaceB, typename HilbertSpaceX, typename Iteration >
+    int solve(const HilbertSpaceB& b, HilbertSpaceX& x, Iteration& iter) const
+    {
+	return cgs(A, x, b, L, iter);
+    }
+
+    /// Perform one CGS iteration on linear system
+    template < typename HilbertSpaceB, typename HilbertSpaceX >
+    int solve(const HilbertSpaceB& b, HilbertSpaceX& x) const
+    {
+	itl::basic_iteration<double> iter(x, 1, 0, 0);
+	return solve(b, x, iter);
+    }
+
+  private:
+    const LinearOperator& A;
+    Preconditioner        L;
+};
+
+} // namespace itl
+
+#endif // ITL_CGS_INCLUDE
+
+
+

AMDiS/lib/mtl4/boost/numeric/itl/krylov/fsm.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+#ifndef ITL_FSM_INCLUDE
+#define ITL_FSM_INCLUDE
+
+#include <boost/numeric/mtl/operation/two_norm.hpp>
+
+namespace itl {
+
+/// Folded spectrum method
+/** Computed and named as in http://en.wikipedia.org/wiki/Folded_spectrum_method **/
+template < typename LinearOperator, typename VectorSpace, typename EigenValue,
+	   typename Damping, typename Iteration >
+int fsm(const LinearOperator& H, VectorSpace& phi, EigenValue eps, Damping alpha, Iteration& iter)
+{
+    VectorSpace v1(H * phi - eps * phi);
+    for (; !iter.finished(v1); ++iter) {
+	VectorSpace v2(H * v1 - eps * v1);
+	phi-= alpha * v2;
+	phi/= two_norm(phi);
+	v1= H * phi - eps * phi;
+    }
+    return iter;
+}
+
+
+} // namespace itl
+
+#endif // ITL_FSM_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/krylov/gmres.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+// Written by Cornelius Steinhardt
+
+
+#ifndef ITL_GMRES_INCLUDE
+#define ITL_GMRES_INCLUDE
+
+#include <algorithm>
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/vector/dense_vector.hpp>
+#include <boost/numeric/mtl/matrix/dense2D.hpp>
+#include <boost/numeric/mtl/matrix/multi_vector.hpp>
+#include <boost/numeric/mtl/operation/givens.hpp>
+#include <boost/numeric/mtl/operation/two_norm.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/utility/irange.hpp>
+
+namespace itl {
+
+/// Generalized Minimal Residual method (without restart)
+/** It computes at most kmax_in iterations (or size(x) depending on what is smaller)
+    regardless on whether the termination criterion is reached or not.   **/
+template < typename Matrix, typename Vector, typename LeftPreconditioner, typename RightPreconditioner, typename Iteration >
+int gmres_full(const Matrix &A, Vector &x, const Vector &b,
+               LeftPreconditioner &L, RightPreconditioner &R, Iteration& iter)
+{
+    using mtl::size; using mtl::irange; using mtl::iall; using std::abs; using std::sqrt;
+    typedef typename mtl::Collection<Vector>::value_type Scalar;
+    typedef typename mtl::Collection<Vector>::size_type  Size;
+
+    if (size(b) == 0) throw mtl::logic_error("empty rhs vector");
+
+    const Scalar                zero= math::zero(Scalar());
+    Scalar                      rho, nu, hr;
+    Size                        k, kmax(std::min(size(x), Size(iter.max_iterations() - iter.iterations())));
+    Vector                      r0(b - A *x), r(solve(L,r0)), va(resource(x)), va0(resource(x)), va00(resource(x));
+    mtl::matrix::multi_vector<Vector>   V(Vector(resource(x), zero), kmax+1);
+    mtl::vector::dense_vector<Scalar>   s(kmax+1, zero), c(kmax+1, zero), g(kmax+1, zero), y(kmax, zero);  // replicated in distributed solvers
+    mtl::matrix::dense2D<Scalar>        H(kmax+1, kmax);                                             // dito
+    H= 0;
+
+    rho= g[0]= two_norm(r);
+    if (iter.finished(rho))
+	return iter;
+    V.vector(0)= r / rho;
+    H= zero;
+
+    // GMRES iteration
+    for (k= 0; k < kmax ; ++k, ++iter) {
+        va0= A * Vector(solve(R, V.vector(k)));
+        V.vector(k+1)= va= solve(L,va0);
+	// orth(V, V[k+1], false);
+        // modified Gram Schmidt method
+        for (Size j= 0; j < k+1; j++) {
+	    H[j][k]= dot(V.vector(j), V.vector(k+1));
+	    V.vector(k+1)-= H[j][k] * V.vector(j);
+        }
+        H[k+1][k]= two_norm(V.vector(k+1));
+        //reorthogonalize
+        for(Size j= 0; j < k+1; j++) {
+	    hr= dot(V.vector(k+1), V.vector(j));
+            H[j][k]+= hr;
+            V.vector(k+1)-= hr * V.vector(j);
+        }
+        H[k+1][k]= two_norm(V.vector(k+1));
+	if (H[k+1][k] != zero)                // watch for breakdown
+            V.vector(k+1)*= 1. / H[k+1][k];
+
+        // k Given's rotations
+	for(Size i= 0; i < k; i++)
+	    mtl::matrix::givens<mtl::matrix::dense2D<Scalar> >(H, H[i][k-1], H[i+1][k-1]).trafo(i);
+
+       nu= sqrt(H[k][k]*H[k][k]+H[k+1][k]*H[k+1][k]);
+       if(nu != zero){
+            c[k]=  H[k][k]/nu;
+            s[k]= -H[k+1][k]/nu;
+            H[k][k]=c[k]*H[k][k]-s[k]*H[k+1][k];
+            H[k+1][k]=0;
+ 	    mtl::vector::givens<mtl::vector::dense_vector<Scalar> >(g, c[k], s[k]).trafo(k);
+        }
+	rho= abs(g[k+1]);
+    }
+
+    //reduce k, to get regular matrix
+    while (k > 0 && abs(g[k-1]<= iter.atol())) k--;
+
+    // iteration is finished -> compute x: solve H*y=g as far as rank of H allows
+    irange                  range(k);
+    for (; !range.empty(); --range) {
+	try {
+	    y[range]= lu_solve(H[range][range], g[range]);
+	} catch (mtl::matrix_singular) { continue; } // if singular then try with sub-matrix
+	break;
+    }
+
+    if (range.finish() < k)
+  	std::cerr << "GMRES orhogonalized with " << k << " vectors but matrix singular, can only use "
+		  << range.finish() << " vectors!\n";
+    if (range.empty())
+        return iter.fail(2, "GMRES did not find any direction to correct x");
+    x+= Vector(solve(R, Vector(V.vector(range)*y[range])));
+
+    r= b - A*x;
+    return iter.terminate(r);
+}
+
+/// Generalized Minimal Residual method with restart
+template < typename Matrix, typename Vector, typename LeftPreconditioner,
+           typename RightPreconditioner, typename Iteration >
+int gmres(const Matrix &A, Vector &x, const Vector &b,
+          LeftPreconditioner &L, RightPreconditioner &R,
+	  Iteration& iter, typename mtl::Collection<Vector>::size_type restart)
+{
+     do {
+	 Iteration inner(iter);
+	 inner.set_max_iterations(std::min(int(iter.iterations()+restart), iter.max_iterations()));
+	 inner.suppress_resume(true);
+	 gmres_full(A, x, b, L, R, inner);
+	 iter.update_progress(inner);
+     } while (!iter.finished());
+
+     return iter;
+}
+
+/// Solver class for GMRES; right preconditioner ignored (prints warning if not identity)
+template < typename LinearOperator, typename Preconditioner= pc::identity<LinearOperator>,
+	   typename RightPreconditioner= pc::identity<LinearOperator> >
+class gmres_solver
+{
+  public:
+    /// Construct solver from a linear operator; generate (left) preconditioner from it
+    explicit gmres_solver(const LinearOperator& A, size_t restart= 8)
+      : A(A), restart(restart), L(A), R(A) {}
+
+    /// Construct solver from a linear operator and left preconditioner
+    gmres_solver(const LinearOperator& A, size_t restart, const Preconditioner& L)
+      : A(A), restart(restart), L(L), R(A) {}
+
+    /// Construct solver from a linear operator and left preconditioner
+    gmres_solver(const LinearOperator& A, size_t restart, const Preconditioner& L, const RightPreconditioner& R)
+      : A(A), restart(restart), L(L), R(R) {}
+
+    /// Solve linear system approximately as specified by \p iter
+    template < typename HilbertSpaceB, typename HilbertSpaceX, typename Iteration >
+    int solve(const HilbertSpaceB& b, HilbertSpaceX& x, Iteration& iter) const
+    {
+	return gmres(A, x, b, L, R, iter, restart);
+    }
+
+    /// Perform one GMRES iteration on linear system
+    template < typename HilbertSpaceB, typename HilbertSpaceX >
+    int solve(const HilbertSpaceB& b, HilbertSpaceX& x) const
+    {
+	itl::basic_iteration<double> iter(x, 1, 0, 0);
+	return solve(b, x, iter);
+    }
+
+  private:
+    const LinearOperator& A;
+    size_t                restart;
+    Preconditioner        L;
+    RightPreconditioner   R;
+};
+
+
+} // namespace itl
+
+#endif // ITL_GMRES_INCLUDE
+
+

AMDiS/lib/mtl4/boost/numeric/itl/krylov/idr_s.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+// Peter Sonneveld and Martin B. van Gijzen, IDR(s): a family of simple and fast algorithms for solving large nonsymmetric linear systems.
+// SIAM J. Sci. Comput. Vol. 31, No. 2, pp. 1035-1062 (2008). (copyright SIAM)
+
+#ifndef ITL_IDR_S_INCLUDE
+#define ITL_IDR_S_INCLUDE
+
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/vector/dense_vector.hpp>
+#include <boost/numeric/mtl/operation/random.hpp>
+#include <boost/numeric/mtl/operation/orth.hpp>
+#include <boost/numeric/mtl/operation/resource.hpp>
+#include <boost/numeric/mtl/matrix/strict_upper.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/mtl/utility/irange.hpp>
+#include <boost/numeric/linear_algebra/identity.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+namespace itl {
+
+/// Induced Dimension Reduction on s dimensions (IDR(s))
+template < typename LinearOperator, typename Vector,
+	   typename LeftPreconditioner, typename RightPreconditioner,
+	   typename Iteration >
+int idr_s(const LinearOperator &A, Vector &x, const Vector &b,
+	  const LeftPreconditioner &, const RightPreconditioner &,
+	  Iteration& iter, size_t s)
+{
+    mtl::vampir_trace<7010> tracer;
+    using mtl::size; using mtl::iall; using mtl::matrix::strict_upper;
+    typedef typename mtl::Collection<Vector>::value_type Scalar;
+    typedef typename mtl::Collection<Vector>::size_type  Size;
+
+    if (size(b) == 0) throw mtl::logic_error("empty rhs vector");
+    if (s < 1) s= 1;
+
+    const Scalar                zero= math::zero(Scalar());
+    Scalar                      omega(zero);
+    Vector                      x0(x), y(resource(x)), v(resource(x)), t(resource(x)), q(resource(x)), r(b - A * x);
+    mtl::matrix::multi_vector<Vector>   dR(Vector(resource(x), zero), s), dX(Vector(resource(x), zero), s), P(Vector(resource(x), zero), s);
+    mtl::vector::dense_vector<Scalar>   m(s), c(s), dm(s);   // replicated in distributed solvers
+    mtl::matrix::dense2D<Scalar>        M(s, s);             // dito
+
+    random(P);
+    P.vector(0)= r;
+    orth(P);
+
+    for (size_t k= 0; k < s; k++) {
+	v= A * r;
+	omega= dot(v, r) / dot(v, v);
+	dX.vector(k)= omega * r;
+	dR.vector(k)= -omega * v;
+	x+= dX.vector(k);
+	r+= dR.vector(k);
+	if ((++iter).finished(r)) return iter;
+	M[iall][k]= trans(P) * dR.vector(k);
+    }
+
+    Size oldest= 0;
+    m= trans(P) * r;
+
+    while (! iter.finished(r)) {
+
+	for (size_t k= 0; k < s; k++) {
+	    c= lu_solve(M, m);
+	    q= dR * -c;
+	    v= r + q;
+	    if (k == 0) {
+		t= A * v;
+		omega= dot(t, v) / dot(t, t);
+		dR.vector(oldest)= q - omega * t;
+		dX.vector(oldest)= omega * v - dX * c;
+	    } else {
+		dX.vector(oldest)= omega * v - dX * c;
+		dR.vector(oldest)= A * -dX.vector(oldest);
+	    }
+	    r+= dR.vector(oldest);
+	    x+= dX.vector(oldest);
+
+	    if ((++iter).finished(r))
+		return iter;
+
+	    dm= trans(P) * dR.vector(oldest);
+	    M[iall][oldest]= dm;
+	    m+= dm;
+	    oldest= (oldest + 1) % s;
+	}
+    }
+    return iter;
+}
+
+/// Solver class for IDR(s) method; right preconditioner ignored (prints warning if not identity)
+template < typename LinearOperator, typename Preconditioner= pc::identity<LinearOperator>,
+	   typename RightPreconditioner= pc::identity<LinearOperator> >
+class idr_s_solver
+{
+  public:
+    /// Construct solver from a linear operator; generate (left) preconditioner from it
+    explicit idr_s_solver(const LinearOperator& A, size_t s= 8) : A(A), s(s), L(A), R(A) {}
+
+    /// Construct solver from a linear operator and left preconditioner
+    idr_s_solver(const LinearOperator& A, size_t s, const Preconditioner& L) : A(A), s(s), L(L), R(A) {}
+
+    /// Construct solver from a linear operator and left preconditioner
+    idr_s_solver(const LinearOperator& A, size_t s, const Preconditioner& L, const RightPreconditioner& R)
+      : A(A), s(s), L(L), R(R) {}
+
+    /// Solve linear system approximately as specified by \p iter
+    template < typename HilbertSpaceB, typename HilbertSpaceX, typename Iteration >
+    int solve(const HilbertSpaceB& b, HilbertSpaceX& x, Iteration& iter) const
+    {
+	return idr_s(A, x, b, L, R, iter, s);
+    }
+
+    /// Perform one IDR(s) iteration on linear system
+    template < typename HilbertSpaceB, typename HilbertSpaceX >
+    int solve(const HilbertSpaceB& b, HilbertSpaceX& x) const
+    {
+	itl::basic_iteration<double> iter(x, 1, 0, 0);
+	return solve(b, x, iter);
+    }
+
+  private:
+    const LinearOperator& A;
+    size_t                s;
+    Preconditioner        L;
+    RightPreconditioner   R;
+};
+
+
+} // namespace itl
+
+#endif // ITL_IDR_S_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/krylov/qmr.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+// Written by Cornelius Steinhardt
+
+
+#ifndef ITL_QMR_INCLUDE
+#define ITL_QMR_INCLUDE
+
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/operation/trans.hpp>
+#include <boost/numeric/mtl/operation/resource.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+namespace itl {
+
+/// Quasi-Minimal Residual method
+template < typename Matrix, typename Vector,typename LeftPreconditioner,
+	   typename RightPreconditioner, typename Iteration >
+int qmr(const Matrix& A, Vector& x, const Vector& b, LeftPreconditioner& L,
+	const RightPreconditioner& R, Iteration& iter)
+{
+    mtl::vampir_trace<7008> tracer;
+    using mtl::size;
+    typedef typename mtl::Collection<Vector>::value_type Scalar;
+    if (size(b) == 0) throw mtl::logic_error("empty rhs vector");
+
+    const Scalar         zero= math::zero(Scalar()), one= math::one(Scalar());
+    Scalar               beta, gamma(one), gamma_1, delta, eta(-one), ep(one), rho_1, theta(zero), theta_1;
+    Vector               r(b - A * x), v_tld(r), y(solve(L, v_tld)), w_tld(r), z(adjoint_solve(R,w_tld)), v(resource(x)), w(resource(x)),
+                         y_tld(resource(x)), z_tld(resource(x)), p(resource(x)), q(resource(x)), p_tld(resource(x)), d(resource(x)), s(resource(x));
+
+    if (iter.finished(r))
+	return iter;
+
+    Scalar rho = two_norm(y), xi = two_norm(z);
+    while(! iter.finished(rho)) {
+	++iter;
+        if (rho == zero)
+	    return iter.fail(1, "qmr breakdown #1, rho=0");
+        if (xi == zero)
+            return iter.fail(2, "qmr breakdown #2, xi=0");
+
+        v= v_tld / rho;
+        y/= rho;
+        w= w_tld / xi;
+        z/= xi;
+
+        delta = dot(z,y);
+        if (delta == zero)
+            return iter.fail(3, "qmr breakdown, delta=0 #3");
+
+        y_tld = solve(R,y);
+        z_tld = adjoint_solve(L,z);
+
+	if (iter.first()) {
+            p = y_tld;
+            q = z_tld;
+	} else {
+            p = y_tld - ((xi * delta) / ep) * p;
+            q = z_tld - ((rho* delta) / ep) * q;
+        }
+
+        p_tld = A * p;
+        ep = dot(q, p_tld);
+        if (ep == zero)
+            return iter.fail(4, "qmr breakdown ep=0 #4");
+        beta= ep / delta;
+        if (beta == zero)
+            return iter.fail(5, "qmr breakdown beta=0 #5");
+        v_tld = p_tld - beta * v;
+        y = solve(L,v_tld);
+        rho_1 = rho;
+	rho = two_norm(y);
+        w_tld= trans(A)*q  - beta*w;
+        z = adjoint_solve(R, w_tld);
+        xi = two_norm(z);
+        gamma_1 = gamma;
+        theta_1 = theta;
+        theta = rho / (gamma_1 * beta);
+        gamma = one / (sqrt(one + theta * theta));
+
+        if (gamma == zero)
+            return iter.fail(6, "qmr breakdown gamma=0 #6");
+
+        eta= -eta * rho_1 * gamma * gamma / (beta * gamma_1 * gamma_1);
+	if (iter.first()) {
+           d= eta * p;
+	   s= eta * p_tld;
+	} else {
+            d= eta * p + (theta_1 * theta_1 * gamma * gamma) * d;
+            s= eta * p_tld + (theta_1 * theta_1 * gamma * gamma) * s;
+        }
+        x += d;
+        r -= s;
+    }
+    return iter;
+}
+
+/// Solver class for Quasi-minimal residual method; right preconditioner ignored (prints warning if not identity)
+template < typename LinearOperator, typename Preconditioner= pc::identity<LinearOperator>,
+	   typename RightPreconditioner= pc::identity<LinearOperator> >
+class qmr_solver
+{
+  public:
+    /// Construct solver from a linear operator; generate (left) preconditioner from it
+    explicit qmr_solver(const LinearOperator& A) : A(A), L(A), R(A) {}
+
+    /// Construct solver from a linear operator and left preconditioner
+    qmr_solver(const LinearOperator& A, const Preconditioner& L) : A(A), L(L), R(A) {}
+
+    /// Construct solver from a linear operator and left preconditioner
+    qmr_solver(const LinearOperator& A, const Preconditioner& L, const RightPreconditioner& R)
+      : A(A), L(L), R(R) {}
+
+    /// Solve linear system approximately as specified by \p iter
+    template < typename HilbertSpaceB, typename HilbertSpaceX, typename Iteration >
+    int solve(const HilbertSpaceB& b, HilbertSpaceX& x, Iteration& iter) const
+    {
+	return qmr(A, x, b, L, R, iter);
+    }
+
+    /// Perform one QMR iteration on linear system
+    template < typename HilbertSpaceB, typename HilbertSpaceX >
+    int solve(const HilbertSpaceB& b, HilbertSpaceX& x) const
+    {
+	itl::basic_iteration<double> iter(x, 1, 0, 0);
+	return solve(b, x, iter);
+    }
+
+  private:
+    const LinearOperator& A;
+    Preconditioner        L;
+    RightPreconditioner   R;
+};
+
+
+} // namespace itl
+
+#endif // ITL_QMR_INCLUDE
+

AMDiS/lib/mtl4/boost/numeric/itl/krylov/repeating_solver.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG, www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also tools/license/license.mtl.txt in the distribution.
+
+#ifndef ITL_REPEATING_SOLVER_INCLUDE
+#define ITL_REPEATING_SOLVER_INCLUDE
+
+#include <boost/mpl/if.hpp>
+#include <boost/static_assert.hpp>
+#include <boost/numeric/itl/itl_fwd.hpp>
+
+namespace itl {
+
+/// Class for calling \tparam N iterations of the given \tparam Solver
+/** If \tparam Stored is true then the \tparam Solver object is stored (i.e. possibly copied) here.
+    Otherwise it is only referred and passing temporary objects to the constructor will cause errors. **/
+template <typename Solver, unsigned N, bool Stored>
+class repeating_solver
+{
+    typedef typename boost::mpl::if_c<Stored, Solver, const Solver&>::type solver_type;
+  public:
+    explicit repeating_solver(const Solver& s) : s(s) {}
+
+    template <typename Matrix>
+    explicit repeating_solver(const Matrix& A) : s(A)
+    {
+	BOOST_STATIC_ASSERT((Stored)); // if matrix is passed class must own solver
+    }
+
+    /// Perform one GMRES iteration on linear system
+    template < typename VectorIn, typename VectorOut >
+    int solve(const VectorIn& b, VectorOut& x) const
+    {
+	itl::basic_iteration<double> iter(x, N, 0, 0);
+	return s.solve(b, x, iter);
+    }
+
+  private:
+    solver_type   s;
+};
+
+} // namespace itl
+
+#endif // ITL_REPEATING_SOLVER_INCLUDE

AMDiS/lib/mtl4/boost/numeric/itl/krylov/tfqmr.hpp

+// Software License for MTL
+//
+// Copyright (c) 2007 The Trustees of Indiana University.
+//               2008 Dresden University of Technology and the Trustees of Indiana University.
+//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
+// All rights reserved.
+// Authors: Peter Gottschling and Andrew Lumsdaine
+//
+// This file is part of the Matrix Template Library
+//
+// See also license.mtl.txt in the distribution.
+
+// Written by Cornelius Steinhardt
+
+
+#ifndef ITL_TFQMR_INCLUDE
+#define ITL_TFQMR_INCLUDE
+
+#include <boost/numeric/mtl/concept/collection.hpp>
+#include <boost/numeric/mtl/utility/exception.hpp>
+#include <boost/numeric/linear_algebra/identity.hpp>
+#include <boost/numeric/linear_algebra/inverse.hpp>
+#include <boost/numeric/mtl/utility/irange.hpp>
+#include <boost/numeric/mtl/operation/resource.hpp>
+#include <boost/numeric/mtl/interface/vpt.hpp>
+
+namespace itl {
+
+/// Transposed-free Quasi-minimal residual
+template < typename Matrix, typename Vector,
+	   typename LeftPreconditioner, typename RightPreconditioner, typename Iteration >
+int tfqmr(const Matrix &A, Vector &x, const Vector &b, const LeftPreconditioner &L,
+	  const RightPreconditioner &R, Iteration& iter)
+{
+    mtl::vampir_trace<7009> tracer;
+    using math::reciprocal; using mtl::size;
+    typedef typename mtl::Collection<Vector>::value_type Scalar;
+
+    if (size(b) == 0) throw mtl::logic_error("empty rhs vector");
+
+    const Scalar                zero= math::zero(Scalar()), one= math::one(Scalar());
+    Scalar                      theta(zero), eta(zero), tau, rho, rhon, sigma, alpha, beta, c;
+    Vector                      rt(b - A*Vector(solve(R, x))) /* shift x= R*x */, r(solve(L, rt)), u1(resource(x)), u2(resource(x)),
+                                y1(resource(x)), y2(resource(x)), w(resource(x)), d(resource(x), zero), v(resource(x));
+
+    if (iter.finished(rt))
+	return iter;
+    y1= w= r;
+    rt= A * Vector(solve(R, y1));
+    u1= v= solve(L,rt);
+    tau= two_norm(r);
+    rho= tau * tau;
+
+    // TFQMR iteration
+    while (! iter.finished(tau)) {
+	++iter;
+	sigma= dot(r,v);
+        if (sigma == zero)
+	    return iter.fail(1, "tfgmr breakdown, sigma=0 #1");
+        alpha= rho / sigma;
+
+        // inner loop
+        for(int j=1; j < 3; j++) {
+            if (j == 1) {
+                w-= alpha * u1;
+                d= y1+ (theta * theta * eta / alpha) * d;
+	    } else {
+                y2= y1 - alpha * v;
+                rt= A * Vector(solve(R, y2));
+                u2= solve(L, rt);
+                w-= alpha * u2;
+                d= y2 + (theta * theta * eta / alpha) * d;
+            }
+            theta= two_norm(w) / tau;
+            c= reciprocal(sqrt(one + theta*theta));
+            tau*= theta * c;
+            eta= c * c * alpha;
+            x+= eta * d;
+        } // end inner loop
+        if (rho == zero)
+            return iter.fail(1, "tfgmr breakdown, rho=0 #2");
+        rhon= dot(r,w);
+        beta= rhon/rho;
+        rho= rhon;
+        y1= w + beta*y2;
+        rt= A * Vector(solve(R, y1));
+        u1= solve(L, rt);
+        v= u1 + beta*(u2 + beta*v);
+        rt= A * x - b;
+    }
+    //shift back
+    x= solve(R, x);
+    return iter;
+}
+
+/// Solver class for Transposed-free quasi-minimal residual method; right preconditioner ignored (prints warning if not identity)
+template < typename LinearOperator, typename Preconditioner= pc::identity<LinearOperator>,
+	   typename RightPreconditioner= pc::identity<LinearOperator> >
+class tfqmr_solver
+{
+  public:
+    /// Construct solver from a linear operator; generate (left) preconditioner from it
+    explicit tfqmr_solver(const LinearOperator& A) : A(A), L(A), R(A) {}
+
+    /// Construct solver from a linear operator and left preconditioner
+    tfqmr_solver(const LinearOperator& A, const Preconditioner& L) : A(A), L(L), R(A) {}
+
+    /// Construct solver from a linear operator and left preconditioner
+    tfqmr_solver(const LinearOperator& A, const Preconditioner& L, const RightPreconditioner& R)
+      : A(A), L(L), R(R) {}
+
+    /// Solve linear system approximately as specified by \p iter
+    template < typename HilbertSpaceB, typename HilbertSpaceX, typename Iteration >
+    int solve(const HilbertSpaceB& b, HilbertSpaceX& x, Iteration& iter) const
+    {
+	return tfqmr(A, x, b, L, R, iter);
+    }
+
+    /// Perform one TFQMR iteration on linear system
+    template < typename HilbertSpaceB, typename HilbertSpaceX >
+    int solve(const HilbertSpaceB& b, HilbertSpaceX& x) const
+    {
+	itl::basic_iteration<double> iter(x, 1, 0, 0);
+	return solve(b, x, iter);
+    }
+
+  private:
+    const LinearOperator& A;
+    Preconditioner        L;
+    RightPreconditioner   R;
+};
+
+} // namespace itl
+
+#endif // ITL_TFQMR_INCLUDE