%
% @file tpardiso.m
%
% @brief
%
% @copyright
% Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
% SPC (Swiss Plasma Center)
%
% SPClibs is free software: you can redistribute it and/or modify it under
% the terms of the GNU Lesser General Public License as published by the Free
% Software Foundation, either version 3 of the License, or (at your option)
% any later version.
%
% SPClibs 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 General Public License for more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with this program. If not, see <https://www.gnu.org/licenses/>.
%
% @author
% (in alphabetical order)
% @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
%
mat='/MAT1';
gbmat;
clear S gb_mat;

% Diagonal balancing of matrix
dbal = 1./sqrt(diag(a));
a = diag(dbal)*a*diag(dbal);

file='tpardiso.h5';
mat='/MAT';
n=hdf5read(file,strcat(mat,'/N'));
nz=hdf5read(file,strcat(mat,'/NZ'));

irow=hdf5read(file,strcat(mat,'/irow'));
cols=hdf5read(file,strcat(mat,'/cols'));
val=hdf5read(file,strcat(mat,'/val'));
perm=hdf5read(file,strcat(mat,'/perm'));

amat=zeros(n,n);

% Check PARDISO mat
for i=1:n
    for k=irow(i):irow(i+1)-1
        j=cols(k);
        amat(i,j) = val(k);
        amat(j,i) = val(k);
    end
end
err = a-amat;
errmx = max(max(abs(err)));
fprintf(1,'Max. error = %e\n', errmx);

figure
spy(sparse(amat(perm,perm)),'r.');
LABEL=sprintf('n = %d, nz =%d', n, nz);
title(LABEL)

% $$$ pmat=zeros(n);
% $$$ for i=1:n
% $$$     pmat(i,perm(i))=1;
% $$$ end
% $$$ amod=pmat*amat*pmat';
% $$$ S=sparse(amod);
% $$$ figure
% $$$ spy(S,'r.');
% $$$ LABEL=sprintf('n = %d', n);
% $$$ title(LABEL)