function varargout = rbfemsolve(Amat, rhs, method, varargin)
%
% res=rbfemsolve(Amat, rhs, method, options ...)
%
% Solving a linear system defined by Amat*res=rhs using various methods
%
% author: Qianqian Fang (q.fang <at> neu.edu)
%
% input:
%     Amat: the left-hand-size matrix, can be sparse
%     rhs:  the right-hand-size vector or matrix (multiple rhs)
%     method: (optional) a string specifying the solving method
%            'mldivide': use the left-divide method: res=Amat\rhs
%            'blqmr': use the block-QMR iterative method for entire RHS
%
%          other supported methods include
%            qmr, tfqmr, cgs, gmres, pcg, cgs, minres, symmlq, bicgstab;
%
%          if method is a positive integer, it calls blqmr.m (part of the
%          blit toolbox) with block size defined by "method".
%
%          for all non-block solvers, adding "par" prefix calls parfor to
%          solve each RHS in parallel, one must call matlabpool or parpool
%          (matlab 2016 or newer) to create the pool first
%     options: (optional) additional parameters are passed to the solver
%          function to specify solver options, please help qmr, cgs, etc to
%          see the additionally accepted parameters
% output:
%
%
% license:
%     GPL version 3, see LICENSE_GPLv3.txt files for details
%
% -- this function is part of Redbird-m toolbox
%

if (nargin < 3)
    method = 'mldivide';
end
if (nargout < 1)
    error('output can not be empty');
end
outputnum = nargout;

if (~ischar(method) && nargout > 0)
    res = zeros(size(rhs));
    if (method < 0)
        method = -method;
        len = 1:method:size(rhs, 2);
        res = cell(1, length(len) - 1);
        parfor i = 1:length(res)
            vo = cell(1, outputnum);
            maxidx = min([size(rhs, 2), len(i) + method - 1]);
            [vo{:}] = blqmr(Amat, full(rhs(:, len(i):maxidx)), varargin{:});
            res{i} = vo{1};
        end
        res = cell2mat(res);
    else
        vo = cell(1, nargout);
        for i = 1:method:size(rhs, 2)
            maxidx = min([size(rhs, 2), i + method - 1]);
            [vo{:}] = blqmr(Amat, full(rhs(:, i:maxidx)), varargin{:});
            res(:, i:maxidx) = vo{1};
        end
    end
    vo{1} = res;
    varargout = vo;
    return
end

% block solvers can handle multiple RHSs
switch method
    case 'blqmr'
        [varargout{1:nargout}] = blqmr(Amat, full(rhs), varargin{:});
        return
    case 'mldivide'
        varargout{1} = full(Amat \ rhs);
        for i = 2:nargout
            varargout{i} = [];
        end
        return
end

% non-block solvers have to solve one RHS at a time

sol = cell(1, size(rhs, 2));
for i = 1:size(rhs, 2)
    sol{i} = cell(1, nargout);
end

if (regexp(method, '^par'))
    method = regexprep(method, '^par', '');
    parfor i = 1:size(rhs, 2)
        switch method
            case 'qmr'
                [sol{i}{:}] = qmr(Amat, rhs(:, i), varargin{:});
            case 'tfqmr'
                [sol{i}{:}] = tfqmr(Amat, rhs(:, i), varargin{:});
            case 'cgs'
                [sol{i}{:}] = cgs(Amat, rhs(:, i), varargin{:});
            case 'lsqr'
                [sol{i}{:}] = lsqr(Amat, rhs(:, i), varargin{:});
            case 'gmres'
                [sol{i}{:}] = gmres(Amat, rhs(:, i), varargin{:});
            case 'gmres'
                [sol{i}{:}] = gmres(Amat, rhs(:, i), varargin{:});
            case 'pcg'
                [sol{i}{:}] = pcg(Amat, rhs(:, i), varargin{:});
            case 'minres'
                [sol{i}{:}] = minres(Amat, rhs(:, i), varargin{:});
            case 'symmlq'
                [sol{i}{:}] = symmlq(Amat, rhs(:, i), varargin{:});
            case 'bicgstab'
                [sol{i}{:}] = bicgstab(Amat, rhs(:, i), varargin{:});
            otherwise
                error(['method "', method, '" is not supported']);
        end
    end
else
    % do not use parfor
    for i = 1:size(rhs, 2)
        switch method
            case 'qmr'
                [sol{i}{:}] = qmr(Amat, rhs(:, i), varargin{:});
            case 'tfqmr'
                [sol{i}{:}] = tfqmr(Amat, rhs(:, i), varargin{:});
            case 'cgs'
                [sol{i}{:}] = cgs(Amat, rhs(:, i), varargin{:});
            case 'lsqr'
                [sol{i}{:}] = lsqr(Amat, rhs(:, i), varargin{:});
            case 'gmres'
                [sol{i}{:}] = gmres(Amat, rhs(:, i), varargin{:});
            case 'gmres'
                [sol{i}{:}] = gmres(Amat, rhs(:, i), varargin{:});
            case 'pcg'
                [sol{i}{:}] = pcg(Amat, rhs(:, i), varargin{:});
            case 'minres'
                [sol{i}{:}] = minres(Amat, rhs(:, i), varargin{:});
            case 'symmlq'
                [sol{i}{:}] = symmlq(Amat, rhs(:, i), varargin{:});
            case 'bicgstab'
                [sol{i}{:}] = bicgstab(Amat, rhs(:, i), varargin{:});
            otherwise
                error(['method "', method, '" is not supported']);
        end
    end
end

res = cell(1, nargout);
for i = 1:nargout
    try
        res{i} = cell2mat(cellfun(@(x) x{i}, sol, 'uniformoutput', true));
    catch
        res{i} = cell2mat(cellfun(@(x) x{i}, sol, 'uniformoutput', false));
    end
end

varargout = res;
