als.h

Go to the documentation of this file.

// Xerus - A General Purpose Tensor Library
// Copyright (C) 2014-2017 Benjamin Huber and Sebastian Wolf. 
// 
// Xerus is free software: you can redistribute it and/or modify
// it under the terms of the GNU Affero General Public License as published
// by the Free Software Foundation, either version 3 of the License,
// or (at your option) any later version.
// 
// Xerus is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU Affero General Public License for more details.
// 
// You should have received a copy of the GNU Affero General Public License
// along with Xerus. If not, see <http://www.gnu.org/licenses/>.
//
// For further information on Xerus visit https://libXerus.org 
// or contact us at contact@libXerus.org.

#pragma once

#include "../ttNetwork.h"
#include "../performanceData.h"

namespace xerus {

    class ALSVariant {
    public:
        enum Direction { Increasing, Decreasing };
        
    protected:
        double solve(const TTOperator *_Ap, TTTensor &_x, const TTTensor &_b, size_t _numHalfSweeps, value_t _convergenceEpsilon, PerformanceData &_perfData = NoPerfData) const;
    
        struct ALSAlgorithmicData {
            struct ContractedTNCache {
                std::vector<Tensor> left, right;
            };
            const ALSVariant &ALS; 
            const TTOperator *A; 
            TTTensor &x; 
            const TTTensor &b; 
            std::vector<size_t> targetRank; 
            ContractedTNCache localOperatorCache; 
            ContractedTNCache rhsCache; 
            value_t normB; 
            std::pair<size_t, size_t> optimizedRange; 
            bool canonicalizeAtTheEnd; 
            size_t corePosAtTheEnd; 
            std::function<value_t()> energy_f; 
            std::function<value_t()> residual_f; 
            
            size_t currIndex; 
            value_t lastEnergy2; 
            value_t lastEnergy; 
            value_t energy; 
            size_t halfSweepCount; 
            Direction direction; 
            
            void prepare_x_for_als();
            
            TensorNetwork localOperatorSlice(size_t _pos);
            TensorNetwork localRhsSlice(size_t _pos);
            
            void prepare_stacks();
            
            void choose_energy_functional();
            
            ALSAlgorithmicData(const ALSVariant &, const TTOperator *, TTTensor &, const TTTensor &);
            
            void move_to_next_index();
        };
        
        TensorNetwork construct_local_operator(ALSAlgorithmicData &_data) const;
        TensorNetwork construct_local_RHS(ALSAlgorithmicData &_data) const;
        bool check_for_end_of_sweep(ALSAlgorithmicData& _data, size_t _numHalfSweeps, value_t _convergenceEpsilon, PerformanceData &_perfData) const;
    public:
        const size_t FLAG_FINISHED_HALFSWEEP = 1;
        const size_t FLAG_FINISHED_FULLSWEEP = 3; // contains the flag for a new halfsweep!
        
        uint sites; 
        size_t numHalfSweeps; 
        value_t convergenceEpsilon; 
//      value_t minimumLocalResidual; ///< below this bound for the local residual, no local solver will be called.
        bool useResidualForEndCriterion; 
        bool preserveCorePosition; 
        bool assumeSPD; 
        
        // TODO std::function endCriterion
        
        using LocalSolver = std::function<void(const TensorNetwork &, std::vector<Tensor> &, const TensorNetwork &, const ALSAlgorithmicData &)>;
        LocalSolver localSolver;
        
        static void lapack_solver(const TensorNetwork &_A, std::vector<Tensor> &_x, const TensorNetwork &_b, const ALSAlgorithmicData &_data);
        static void ASD_solver(const TensorNetwork &_A, std::vector<Tensor> &_x, const TensorNetwork &_b, const ALSAlgorithmicData &_data);
        
        //TODO add local CG solver
        //TODO add AMEn solver
        
        ALSVariant(uint _sites, size_t _numHalfSweeps, LocalSolver _localSolver, bool _assumeSPD,
            bool _useResidual=false
        ) 
                : sites(_sites), numHalfSweeps(_numHalfSweeps), convergenceEpsilon(1e-6), 
                useResidualForEndCriterion(_useResidual), preserveCorePosition(true), assumeSPD(_assumeSPD), localSolver(_localSolver)
        {
            XERUS_REQUIRE(_sites>0, "");
        }
        
        double operator()(const TTOperator &_A, TTTensor &_x, const TTTensor &_b, value_t _convergenceEpsilon, PerformanceData &_perfData = NoPerfData) const {
            return solve(&_A, _x, _b, numHalfSweeps, _convergenceEpsilon, _perfData);
        }
        
        double operator()(const TTOperator &_A, TTTensor &_x, const TTTensor &_b, size_t _numHalfSweeps, PerformanceData &_perfData = NoPerfData) const {
            return solve(&_A, _x, _b, _numHalfSweeps, convergenceEpsilon, _perfData);
        }
        
        double operator()(const TTOperator &_A, TTTensor &_x, const TTTensor &_b, PerformanceData &_perfData = NoPerfData) const {
            return solve(&_A, _x, _b, numHalfSweeps, convergenceEpsilon, _perfData);
        }
        
        double operator()(TTTensor &_x, const TTTensor &_b, value_t _convergenceEpsilon, PerformanceData &_perfData = NoPerfData) const {
            return solve(nullptr, _x, _b, numHalfSweeps, _convergenceEpsilon, _perfData);
        }
        
        double operator()(TTTensor &_x, const TTTensor &_b, size_t _numHalfSweeps, PerformanceData &_perfData = NoPerfData) const {
            return solve(nullptr, _x, _b, _numHalfSweeps, convergenceEpsilon, _perfData);
        }
        
        double operator()(TTTensor &_x, const TTTensor &_b, PerformanceData &_perfData = NoPerfData) const {
            return solve(nullptr, _x, _b, numHalfSweeps, convergenceEpsilon, _perfData);
        }
    };
    
    extern const ALSVariant ALS;
    extern const ALSVariant ALS_SPD;
    
    extern const ALSVariant DMRG;
    extern const ALSVariant DMRG_SPD;
    
    extern const ALSVariant ASD;
    extern const ALSVariant ASD_SPD;
}