tsgLinearSolvers.hpp

/*
 * Copyright (c) 2017, Miroslav Stoyanov
 *
 * This file is part of
 * Toolkit for Adaptive Stochastic Modeling And Non-Intrusive ApproximatioN: TASMANIAN
 *
 * Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions
 *    and the following disclaimer in the documentation and/or other materials provided with the distribution.
 *
 * 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse
 *    or promote products derived from this software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
 * IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
 * OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
 * OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
 * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *
 * UT-BATTELLE, LLC AND THE UNITED STATES GOVERNMENT MAKE NO REPRESENTATIONS AND DISCLAIM ALL WARRANTIES, BOTH EXPRESSED AND IMPLIED.
 * THERE ARE NO EXPRESS OR IMPLIED WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, OR THAT THE USE OF THE SOFTWARE WILL NOT INFRINGE ANY PATENT,
 * COPYRIGHT, TRADEMARK, OR OTHER PROPRIETARY RIGHTS, OR THAT THE SOFTWARE WILL ACCOMPLISH THE INTENDED RESULTS OR THAT THE SOFTWARE OR ITS USE WILL NOT RESULT IN INJURY OR DAMAGE.
 * THE USER ASSUMES RESPONSIBILITY FOR ALL LIABILITIES, PENALTIES, FINES, CLAIMS, CAUSES OF ACTION, AND COSTS AND EXPENSES, CAUSED BY, RESULTING FROM OR ARISING OUT OF,
 * IN WHOLE OR IN PART THE USE, STORAGE OR DISPOSAL OF THE SOFTWARE.
 */
 
#ifndef __TASMANIAN_LINEAR_SOLVERS_HPP
#define __TASMANIAN_LINEAR_SOLVERS_HPP
 
#include "tsgAcceleratedDataStructures.hpp"


 
namespace TasGrid{

namespace TasmanianDenseSolver{
    void solveLeastSquares(int n, int m, const double A[], double b[], double *x);
    void solveLeastSquares(AccelerationContext const *acceleration, int n, int m, double A[], double b[], double *x);
    template<typename scalar_type>
    void solvesLeastSquares(AccelerationContext const *acceleration, int n, int m, scalar_type A[], int nrhs, scalar_type B[]);
    template<typename scalar_type>
    void solvesLeastSquaresGPU(AccelerationContext const *acceleration, int n, int m, scalar_type A[], int nrhs, scalar_type B[]);
}

namespace TasmanianTridiagonalSolver{
    std::vector<double> getSymmetricEigenvalues(int n, std::vector<double> const &diag, std::vector<double> const &offdiag);

    static constexpr int decompose_version = 1;

    void decompose(std::vector<double> &diag, std::vector<double> &off_diag, const double mu0, std::vector<double> &nodes,
                   std::vector<double> &weights);

    void decompose1(int n, std::vector<double> &d, std::vector<double> &e, std::vector<double> &z);

    void decompose2(std::vector<double> &diag, std::vector<double> &off_diag, const double mu0, std::vector<double> &nodes,
                    std::vector<double> &weights);
}

namespace TasmanianFourierTransform{

    void fast_fourier_transform(std::vector<std::vector<std::complex<double>>> &data, std::vector<int> &num_points);


    void fast_fourier_transform1D(std::vector<std::vector<std::complex<double>>> &data, std::vector<int> &indexes);
}

namespace TasSparse{

class WaveletBasisMatrix{
public:
    WaveletBasisMatrix() : tol(Maths::num_tol), num_rows(0){}
    WaveletBasisMatrix(AccelerationContext const *acceleration, const std::vector<int> &lpntr, const std::vector<std::vector<int>> &lindx, const std::vector<std::vector<double>> &lvals);
    WaveletBasisMatrix(AccelerationContext const *acceleration, int cnum_rows, GpuVector<double> &&matrix)
        : tol(Maths::num_tol), num_rows(cnum_rows), gpu_dense(std::move(matrix)) { factorize(acceleration); }
    ~WaveletBasisMatrix() = default;

    WaveletBasisMatrix(WaveletBasisMatrix const&) = delete;
    WaveletBasisMatrix(WaveletBasisMatrix &&) = default;
    WaveletBasisMatrix& operator =(WaveletBasisMatrix const&) = delete;
    WaveletBasisMatrix& operator =(WaveletBasisMatrix &&) = default;

    static bool useDense(AccelerationContext const *acceleration, int nrows){
        return ((acceleration->algorithm_select != AccelerationContext::algorithm_sparse)
                and not (acceleration->algorithm_select == AccelerationContext::algorithm_autoselect and nrows > 10000));
    }

    bool isSparse() const{ return (num_rows > 0) and dense.empty(); }
    bool isDense() const{ return (num_rows > 0) and not dense.empty(); }

    int getNumRows() const{ return num_rows; }

    void invertTransposed(AccelerationContext const *acceleration, double b[]) const;

    void invert(AccelerationContext const *acceleration, int num_colums, double B[]);

    template<bool transpose, bool blas> void solve(const double b[], double x[]) const;
 
protected:
    void factorize(AccelerationContext const *acceleration);
    void computeILU();
    template<bool transpose> void applyILU(double x[]) const;
    template<bool transpose> void apply(double const x[], double r[]) const;
    void residual(double const x[], double const b[], double r[]) const;
    static constexpr bool use_transpose = true;
    static constexpr bool no_transpose = false;
    static constexpr bool use_blas = true;
    static constexpr bool no_blas = false;
 
private:
    double tol;
    int num_rows;
    std::vector<int> pntr, indx, indxD;
    std::vector<double> vals, ilu;
    std::vector<double> dense; // if using the dense format
    std::vector<int> ipiv; // pivots for the dense factorize
 
    GpuVector<double> gpu_dense; // GPU factors
    GpuVector<int_gpu_lapack> gpu_ipiv;
};
 
}
 
}
 
#endif