test_packages_dke_yp

LUKE - Test function of matrix inversion techniques for LUKE

 by Y. Peysson (CEA-IRFM) <yves.peysson@cea.fr> and J. Decker (CEA-IRFM) <joan.decker@cea.fr>

0001 %LUKE - Test function of matrix inversion techniques for LUKE
0002 %
0003 % by Y. Peysson (CEA-IRFM) <yves.peysson@cea.fr> and J. Decker (CEA-IRFM) <joan.decker@cea.fr>
0004 %
0005 clc
0006 clear all
0007 close all
0008 %
0009 % Flags of the tested methods (0: not testes, 1: tested). WARNING: it tests
0010 % all compiled version of a matrix inversion method. Check before the
0011 % compatibility with the working release of Matlab.
0012 %
0013 MUMPS = 1;
0014 PETSC = 0;
0015 SUPERLU = 1;
0016 PARDISO = 0;
0017 %
0018 [lukeversion,matver] = LUKEversion_jd;
0019 %
0020 try
0021     GPUstatus = gpuDevice;%data from the GPU device
0022 catch
0023     GPUstatus.DeviceSupported = false;
0024 end    
0025 %
0026 if GPUstatus.DeviceSupported,
0027     reset(parallel.gpu.GPUDevice.current());%reset memory on GPU device (should be done to have the largest memory available)
0028     GPUstatus = gpuDevice;%data from the GPU device
0029 end
0030 %
0031 id_dkeparam = 'NONUNIFORM10010020';%For DKE code parameters
0032 path_dkeparam = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0033 %
0034 [dkepath,dkeparam] = load_structures_yp('dkepath','','','dkeparam',id_dkeparam,path_dkeparam);
0035 dirinit = pwd;%initial directory
0036 cd([dkepath.luke_root,'Project_DKE/Packages']);
0037 %
0038 disp(' ');
0039 disp('Test performances of matrix inversion routines on your computer')
0040 disp('---------------------------------------------------------------')
0041 if GPUstatus.DeviceSupported,
0042     mode = input_dke_yp('2-D (0), 2-D for GPU (1) or 3-D (2) input matrix',0);
0043 else
0044     mode = input_dke_yp('2-D (0) or 3-D (1) input matrix',0);
0045     if mode == 1,
0046         mode = 2;
0047     end
0048 end
0049 %
0050 npass = input_dke_yp('Number of passes',1);
0051 %
0052 if GPUstatus.DeviceSupported && mode==1,
0053     gpu_mode = input_dke_yp('Calculation with GPU',1);disp(' ');
0054 elseif GPUstatus.DeviceSupported && mode~=1,
0055     gpu_mode = input_dke_yp('Calculation with GPU',0);disp(' ');    
0056 end
0057 %
0058 disp(' ');
0059 disp('---------------------------------------------------------------')
0060 %
0061 if mode ==0,
0062    load('testLU2D.mat','MMXPC_f_t','MMXR_t');%,'dkeparam'
0063    if exist('testmatinv-2D.txt','file'),
0064        delete testmatinv-2D.txt;
0065    end;
0066    diary testmatinv-2D.txt;
0067    disp(['2-D LUKE matrix inversion performances on ',computer]);disp('----------------------------------------------')
0068 elseif mode == 1
0069    load('testLU2D_GPU.mat','MMXPC_f_t','MMXR_t');%,'dkeparam'
0070    if exist('testmatinv-2D_GPU.txt','file'),
0071        delete testmatinv-2D_GPU.txt;
0072    end;
0073    diary testmatinv-2D_GPU.txt;
0074    disp(['2-D (GPU) LUKE matrix inversion performances on ',computer]);disp('----------------------------------------------')
0075 elseif mode == 2
0076    load('testLU3D.mat','MMXPC_f_t','MMXR_t');%,'dkeparam'
0077    if exist('testmatinv-3D.txt','file'),
0078        delete testmatinv-3D.txt;
0079    end;
0080    diary testmatinv-3D.txt;
0081    disp(['3-D LUKE matrix inversion performances on ',computer]);disp('----------------------------------------------')
0082 end
0083 %
0084 if GPUstatus.DeviceSupported && gpu_mode,
0085     addpath([dkepath.luke_root,'/Project_DKE/Packages/MDCE/GPU']);
0086 end
0087 %
0088 dkepath.temppath = [dkepath.temppathroot,filesep,timeid_jd(clock)];
0089 if isfield(dkepath,'local_root'),
0090     local_root = dkepath.local_root;
0091 else
0092     local_root = dkepath.luke_root;
0093 end
0094 %
0095 % Matrix inversion tests
0096 %
0097 disp(' ');disp(['Matrix size: (',int2str(length(MMXR_t)),'x',int2str(length(MMXR_t)),')']);disp(' ');disp('---------------------------------------------------------------')
0098 %
0099 % MatLab built-in solvers
0100 %
0101 ipass = 0;
0102 %
0103 ilu_param.type = 'crout';%equivalent to old matlab function luinc
0104 ilu_param.droptol = 1e-6;
0105 ilu_param.milu = 'row';
0106 %
0107 if mode < 2,%not too large matrices for GPU tests and direct matrix inversions
0108     if GPUstatus.DeviceSupported && gpu_mode,
0109         MMXPC_f_t_full = single(full(MMXPC_f_t));%full matrix for GPU tests (single for memory limitations)
0110         MMXR_t_full = single(full(MMXR_t));
0111     else
0112         MMXPC_f_t_full = full(MMXPC_f_t);%full matrix for GPU tests
0113         MMXR_t_full = full(MMXR_t);
0114     end
0115 end
0116 %
0117 while ipass < npass,
0118     %
0119     isim = 0;  
0120     %
0121     % Direct methods
0122     %
0123     isim = isim + 1;disp('- MATLAB-Direct inversion (\) by CPU with sparse matrix'),tic,MMXS_t = MMXPC_f_t\MMXR_t;t(isim) = toc;disp(['Elapsed time is ',num2str(t(isim)),' seconds per inversion']);
0124     Zname{isim} = ['MATLAB_DIRECT_INVERSION'];
0125     %
0126     if mode < 2,%not too large matrices for GPU tests and direct matrix inversions
0127         isim = isim + 1;
0128         if GPUstatus.DeviceSupported && gpu_mode,
0129             disp('- MATLAB-Direct inversion (\) by CPU with full matrix (single precision)'),
0130             tic,MMXS_t_full = invmatGPU(MMXPC_f_t_full,MMXR_t_full);t(isim) = toc;disp(['Elapsed time is ',num2str(t(isim)),' seconds per inversion']);
0131             Zname{isim} = ['MATLAB_DIRECT_INVERSION_FULL_MATRIX_SINGLE'];
0132         else
0133             disp('- MATLAB-Direct inversion (\) by CPU with full matrix'),
0134             tic,MMXS_t_full = MMXPC_f_t_full\MMXR_t_full;t(isim) = toc;disp(['Elapsed time is ',num2str(t(isim)),' seconds per inversion']);
0135             Zname{isim} = ['MATLAB_DIRECT_INVERSION_FULL_MATRIX'];
0136         end
0137         %
0138         if GPUstatus.DeviceSupported && gpu_mode,
0139             isim = isim + 1;tic,dummy = gather(invmatGPU(gpuArray(rand(2,2,'single')),gpuArray(rand(2,1,'single'))));t(isim) = toc;%GPU initialisation
0140             Zname{isim} = ['MATLAB_GPU_INITIALISATION'];
0141             disp('- MATLAB-Direct inversion (\) by GPU with full matrix'),
0142             isim = isim + 1;tic,MMXS_t_full_GPU = gather(invmatGPU(gpuArray(MMXPC_f_t_full),gpuArray(MMXR_t_full)));t(isim) = toc;disp(['Elapsed time is ',num2str(t(isim)),' seconds per inversion + GPU initialisation (once): ',num2str(t(isim-1)),' seconds.']);
0143             Zname{isim} = ['MATLAB_DIRECT_GPU_INVERSION_FULL_MATRIX_SINGLE'];
0144             %
0145             rmpath([dkepath.luke_root,'/Project_DKE/Packages/MDCE/GPU']);
0146         end
0147     end
0148     %
0149     disp('---------------------------------------------------------------')
0150     %
0151     % Iterative methods
0152     %
0153     isim = isim + 1;disp('- MATLAB-ILU by CPU with sparse matrix (drop_tol = 1e-6)'),tic,[MMXL_f_t,MMXU_f_t] = ilu(MMXPC_f_t,ilu_param);t(isim) = toc;
0154     Zname{isim} = ['MATLAB_LU_FACTORISATION'];
0155     isim = isim + 1;disp('- MATLAB-CGS by CPU with sparse matrix'),tic,[MMXS_t_cgs,flaginv_f] = cgs(MMXPC_f_t,MMXR_t,dkeparam.mlprecinv,dkeparam.mlmaxitinv,MMXL_f_t,MMXU_f_t);t(isim) = toc;disp(['Elapsed time is ',num2str(t(isim)),' seconds per inversion + ILU time (once): ',num2str(t(isim-1)),' seconds.']);
0156     ZMMXS_t{isim} = MMXS_t_cgs;
0157     Zname{isim} = ['MATLAB_GCS'];
0158     isim = isim + 1;disp('- MATLAB-BICG by CPU with sparse matrix'),tic,[MMXS_t_bicg,flaginv_f] = bicg(MMXPC_f_t,MMXR_t,dkeparam.mlprecinv,dkeparam.mlmaxitinv,MMXL_f_t,MMXU_f_t);t(isim) = toc;disp(['Elapsed time is ',num2str(t(isim)),' seconds per inversion + ILU time (once): ',num2str(t(isim-2)),' seconds.']);
0159     ZMMXS_t{isim} = MMXS_t_bicg;
0160     Zname{isim} = ['MATLAB_BICG'];
0161     isim = isim + 1;disp('- MATLAB-QMR by CPU with sparse matrix'),tic,[MMXS_t_qmr,flaginv_f] = qmr(MMXPC_f_t,MMXR_t,dkeparam.mlprecinv,dkeparam.mlmaxitinv,MMXL_f_t,MMXU_f_t);t(isim) = toc;disp(['Elapsed time is ',num2str(t(isim)),' seconds per inversion + ILU time (once): ',num2str(t(isim-3)),' seconds.']);
0162     ZMMXS_t{isim} = MMXS_t_qmr;
0163     Zname{isim} = ['MATLAB_QMR'];    
0164     disp('---------------------------------------------------------------')
0165 %
0166 % MUMPS solver
0167 %
0168 % MUMPS packages (MEX application)
0169 %
0170     MUMPSDIRS = {[local_root,'Packages/Bin/MUMPS/'],local_root};
0171     MUMPSDIR = '';
0172     for idir = 1:length(MUMPSDIRS),
0173         if exist(MUMPSDIRS{idir},'dir'),
0174             MUMPSDIR = MUMPSDIRS{idir};
0175             break
0176         end
0177     end
0178     %
0179     if ~isempty(MUMPSDIR) & MUMPS,
0180         list_of_packages = dir(MUMPSDIR);
0181         for ifile = 1:length(list_of_packages),
0182             if length(list_of_packages(ifile).name) >= 5 && strcmp(list_of_packages(ifile).name(1:5),'MUMPS') && ~isempty(strfind(list_of_packages(ifile).name,'_S_')),
0183                 MUMPSMATLAB =[MUMPSDIR,list_of_packages(ifile).name,filesep,'MATLAB',matver,filesep];
0184                 if exist(MUMPSMATLAB,'dir') && exist([MUMPSMATLAB,'dmumpsmex.',mexext],'file'),
0185                     addpath(MUMPSMATLAB)
0186                     isim = isim + 1;disp(['- ',list_of_packages(ifile).name,'-ILU (MEX)']),tic,job_MUMPSMEX = 1;MUMPSMEX_dke_yp;job_MUMPSMEX = 2;MUMPSMEX_dke_yp;t(isim) = toc;
0187                     Zname{isim} = ['MEX_',list_of_packages(ifile).name,'_LU_FACTORISATION'];
0188                     isim = isim + 1;disp(['- ',list_of_packages(ifile).name,'-SOLVE (MEX)']),tic,job_MUMPSMEX = 3;MUMPSMEX_dke_yp;t(isim) = toc;MMXS_t_MUMPS = id_MUMPSMEX.SOL;disp(['Elapsed time is ',num2str(t(isim)),' seconds per inversion + ILU time (once): ',num2str(t(isim-1)),' seconds.']);job_MUMPSMEX = -2;MUMPSMEX_dke_yp;
0189                     ZMMXS_t{isim} = MMXS_t_MUMPS;
0190                     Zname{isim} = ['MEX_',list_of_packages(ifile).name];
0191                     rmpath(MUMPSMATLAB)  
0192                 end    
0193             end
0194         end
0195     end
0196 %
0197     disp('---------------------------------------------------------------')
0198 %
0199 % MUMPS packages (stand-alone application)
0200 %
0201     if exist([local_root,'Packages/Bin/MUMPS']) & MUMPS,
0202         list_of_packages = dir([local_root,'Packages/Bin/MUMPS']);
0203         for ifile = 1:length(list_of_packages),
0204             if length(list_of_packages(ifile).name) >= 5 & list_of_packages(ifile).name(1:5) == 'MUMPS',
0205 %
0206 % Sequential processing (stand-alone application)
0207 %
0208                 if ~isempty(strfind(list_of_packages(ifile).name,'_S_'))
0209                     dkepath.MUMPS.seq = [local_root,'Packages/Bin/MUMPS/',list_of_packages(ifile).name,'/dke'];
0210                     if isfield(dkepath.MUMPS,'par'),rmfield(dkepath.MUMPS,'par');end
0211                     if exist(dkepath.MUMPS.seq),  
0212                         try
0213                             isim = isim + 1;disp(['- ',list_of_packages(ifile).name,' by CPU with sparse matrix (sequential processing, stand-alone application)']),tic,[MMXS_t_MUMPS,flaginv_f] = MUMPS_dke_yp(dkeparam.MUMPSparam,MMXPC_f_t,full(MMXR_t),0,dkepath);t(isim) = toc;disp(['Elapsed time is ',num2str(t(isim)),' seconds.']);  
0214                         catch
0215                             MMXS_t_MUMPS = [];
0216                         end
0217                         ZMMXS_t{isim} = MMXS_t_MUMPS;
0218                         Zname{isim} = list_of_packages(ifile).name;
0219                         clear MMXS_t_MUMPS
0220                     end
0221                 end
0222 %
0223 % Parallel processing (stand-alone application)
0224 %
0225                 if ~isempty(strfind(list_of_packages(ifile).name,'_P_'))
0226                     dkepath.MUMPS.par = [local_root,'Packages/Bin/MUMPS/',list_of_packages(ifile).name,'/dke'];
0227                     if isfield(dkepath.MUMPS,'seq'),rmfield(dkepath.MUMPS,'seq');end
0228                     if exist(dkepath.MUMPS.par),                                   
0229                         dkeparam.MUMPSparam.nproc = 4;%case with 4 processors
0230                         try
0231                             isim = isim + 1;disp(['- ',list_of_packages(ifile).name,', np = 4 (parallel processing, stand-alone application)']),tic,[MMXS_t_MUMPS,flaginv_f] = MUMPS_dke_yp(dkeparam.MUMPSparam,MMXPC_f_t,full(MMXR_t),0,dkepath);t(isim) = toc;disp(['Elapsed time is ',num2str(t(isim)),' seconds.']);
0232                         catch 
0233                             MMXS_t_MUMPS = [];
0234                         end
0235                         ZMMXS_t{isim} = MMXS_t_MUMPS;
0236                         Zname{isim} = [list_of_packages(ifile).name,'_np_4'];
0237                     %
0238                         dkeparam.MUMPSparam.nproc = 1;%case with 1 processor
0239                         isim = isim + 1;disp(['- ',list_of_packages(ifile).name,', np = 1 (parallel processing, stand-alone application)']),tic,[MMXS_t_MUMPS,flaginv_f] = MUMPS_dke_yp(dkeparam.MUMPSparam,MMXPC_f_t,full(MMXR_t),0,dkepath);t(isim) = toc;disp(['Elapsed time is ',num2str(t(isim)),' seconds.']);
0240                         ZMMXS_t{isim} = MMXS_t_MUMPS;
0241                         Zname{isim} = [list_of_packages(ifile).name,'_np_1'];
0242                     end
0243                 end
0244             end
0245         end
0246     end           
0247     disp('---------------------------------------------------------------')
0248 %
0249 % PETSc packages (stand-alone application)
0250 %
0251     if exist([local_root,'Packages/Bin/PETSC']) & PETSC,
0252         list_of_packages = dir([local_root,'Packages/Bin/PETSC']);
0253         for ifile = 1:length(list_of_packages),   
0254             if length(list_of_packages(ifile).name) >= 5 & list_of_packages(ifile).name(1:5) == 'petsc',
0255 
0256 %
0257 % Sequential processing (stand-alone application)
0258 %
0259                 if ~isempty(strfind(list_of_packages(ifile).name,'_S_'))
0260                     dkepath.PETSc.seq = [local_root,'Packages/Bin/PETSC/',list_of_packages(ifile).name,'/DKE'];
0261                     if isfield(dkepath.PETSc,'par'),rmfield(dkepath.PETSc,'par');end
0262                     if exist(dkepath.PETSc.seq),
0263                         dkeparam.PETScparam.nproc = 1;%case with 1 processor
0264                         dkeparam.PETScparam.matrixtype = '';%Matrix type used by PETSc. For SUPERLU use '-matload_type superlu_dist' otherwise ' '. See PETSc documentation
0265                         dkeparam.PETScparam.kspmethod = '-ksp_type bcgs';%KSP method used by PETSc. '-ksp_type preonly': no iteration, direct matrix solve. Otherwise '-ksp_type bicg' or '-ksp_type bcgs'  for example. See PETSc documentation
0266                         dkeparam.PETScparam.pcmethod = '';%Preconditioning method used by PETSc. '-pc_type lu' for SUPERLU or '-pc-type hypre' for HYPRE.  See PETSc documentation
0267                         dkeparam.PETScparam.cloop = 0;%First order collision selfconsistency is done internaly in C (1) or in MatLab (0)
0268                         isim = isim + 1;disp(['- ',list_of_packages(ifile).name,' bcgs (sequential processing, stand-alone application)']),tic,[MMXS_t_PETSc,flaginv_f] = PETSc_dke_ft(dkeparam.PETScparam,MMXPC_f_t,full(MMXR_t),0,dkepath);t(isim) = toc;disp(['Elapsed time is ',num2str(t(isim)),' seconds.']);
0269                         ZMMXS_t{isim} = MMXS_t_PETSc;
0270                         Zname{isim} = [list_of_packages(ifile).name,'_np_1_bgcs_Mloop'];
0271                     end
0272                 end
0273 %
0274 % Parallel processing (stand-alone application)
0275 %
0276                 if ~isempty(strfind(list_of_packages(ifile).name,'_P_'))
0277                     dkepath.PETSc.par = [local_root,'Packages/Bin/PETSC/',list_of_packages(ifile).name,'/DKE'];
0278                     if isfield(dkepath.PETSc,'seq'),rmfield(dkepath.PETSc,'seq');end
0279                     if exist(dkepath.PETSc.par),
0280                         dkeparam.PETScparam.nproc = 4;%case with 4 processors
0281                         dkeparam.PETScparam.matrixtype = '';%Matrix type used by PETSc. For SUPERLU use '-matload_type superlu_dist' otherwise ' '. See PETSc documentation
0282                         dkeparam.PETScparam.kspmethod = '-ksp_type bcgs';%KSP method used by PETSc. '-ksp_type preonly': no iteration, direct matrix solve. Otherwise '-ksp_type bicg' or '-ksp_type bcgs'  for example. See PETSc documentation
0283                         dkeparam.PETScparam.pcmethod = '';%Preconditioning method used by PETSc. '-pc_type lu' for SUPERLU or '-pc-type hypre' for HYPRE.  See PETSc documentation
0284                         dkeparam.PETScparam.cloop = 0;%First order collision selfconsistency is done internaly in C (1) or in MatLab (0)
0285                         isim = isim + 1;disp(['- ',list_of_packages(ifile).name,', np = 4 bcgs (parallel processing, stand-alone application)']),tic,[MMXS_t_PETSc,flaginv_f] = PETSc_dke_ft(dkeparam.PETScparam,MMXPC_f_t,full(MMXR_t),0,dkepath);t(isim) = toc;disp(['Elapsed time is ',num2str(t(isim)),' seconds.']);
0286                         ZMMXS_t{isim} = MMXS_t_PETSc;
0287                         Zname{isim} = [list_of_packages(ifile).name,'_np_4_bgcs_Mloop'];
0288                     %
0289                         dkeparam.PETScparam.nproc = 1;%case with 1 processor
0290                         dkeparam.PETScparam.matrixtype = '';%Matrix type used by PETSc. For SUPERLU use '-matload_type superlu_dist' otherwise ' '. See PETSc documentation
0291                         dkeparam.PETScparam.kspmethod = '-ksp_type bcgs';%KSP method used by PETSc. '-ksp_type preonly': no iteration, direct matrix solve. Otherwise '-ksp_type bicg' or '-ksp_type bcgs'  for example. See PETSc documentation
0292                         dkeparam.PETScparam.pcmethod = '';%Preconditioning method used by PETSc. '-pc_type lu' for SUPERLU or '-pc-type hypre' for HYPRE.  See PETSc documentation
0293                         dkeparam.PETScparam.cloop = 0;%First order collision selfconsistency is done internaly in C (1) or in MatLab (0)
0294                         isim = isim + 1;disp(['- ',list_of_packages(ifile).name,', np = 1 bcgs (parallel processing, stand-alone application)']),tic,[MMXS_t_PETSc,flaginv_f] = PETSc_dke_ft(dkeparam.PETScparam,MMXPC_f_t,full(MMXR_t),0,dkepath);t(isim) = toc;disp(['Elapsed time is ',num2str(t(isim)),' seconds.']);
0295                         ZMMXS_t{isim} = MMXS_t_PETSc;
0296                         Zname{isim} = [list_of_packages(ifile).name,'_np_1_bgcs_Mloop'];
0297 %
0298                         [flag,message] = unix('uname -n');
0299                         if strfind(message,'saturne'),%this option runs only on the SATURNE CEA-DRFC cluster
0300                             dkeparam.PETScparam.nproc = 1;%case with 1 processor
0301                             dkeparam.PETScparam.matrixtype = '-matload_type superlu_dist';%Matrix type used by PETSc. For SUPERLU use '-matload_type superlu_dist' otherwise ' '. See PETSc documentation
0302                             dkeparam.PETScparam.kspmethod = '-ksp_type preonly';%KSP method used by PETSc. '-ksp_type preonly': no iteration, direct matrix solve. Otherwise '-ksp_type bicg' or '-ksp_type bcgs'  for example. See PETSc documentation
0303                             dkeparam.PETScparam.pcmethod = '-pc_type lu';%Preconditioning method used by PETSc. '-pc_type lu' for SUPERLU or '-pc-type hypre' for HYPRE.  See PETSc documentation
0304                             dkeparam.PETScparam.cloop = 0;%First order collision selfconsistency is done internaly in C (1) or in MatLab (0)
0305                             isim = isim + 1;disp(['- ',list_of_packages(ifile).name,', np = 1 superlu_dist (parallel processing, stand-alone application)']),tic,[MMXS_t_PETSc,flaginv_f] = PETSc_dke_ft(dkeparam.PETScparam,MMXPC_f_t,full(MMXR_t),0,dkepath);t(isim) = toc;disp(['Elapsed time is ',num2str(t(isim)),' seconds.']);
0306                             ZMMXS_t{isim} = MMXS_t_PETSc;
0307                             Zname{isim} = [list_of_packages(ifile).name,'_np_1_superludist_Mloop'];
0308 %
0309                             dkeparam.PETScparam.nproc = 4;%case with 1 processor
0310                             dkeparam.PETScparam.matrixtype = '-matload_type superlu_dist';%Matrix type used by PETSc. For SUPERLU use '-matload_type superlu_dist' otherwise ' '. See PETSc documentation
0311                             dkeparam.PETScparam.kspmethod = '-ksp_type preonly';%KSP method used by PETSc. '-ksp_type preonly': no iteration, direct matrix solve. Otherwise '-ksp_type bicg' or '-ksp_type bcgs'  for example. See PETSc documentation
0312                             dkeparam.PETScparam.pcmethod = '-pc_type lu';%Preconditioning method used by PETSc. '-pc_type lu' for SUPERLU or '-pc-type hypre' for HYPRE.  See PETSc documentation
0313                             dkeparam.PETScparam.cloop = 0;%First order collision selfconsistency is done internaly in C (1) or in MatLab (0)
0314                             isim = isim + 1;disp(['- ',list_of_packages(ifile).name,', np = 4 superlu_dist (parallel processing, stand-alone application)']),tic,[MMXS_t_PETSc,flaginv_f] = PETSc_dke_ft(dkeparam.PETScparam,MMXPC_f_t,full(MMXR_t),0,dkepath);t(isim) = toc;disp(['Elapsed time is ',num2str(t(isim)),' seconds.']);
0315                             ZMMXS_t{isim} = MMXS_t_PETSc;
0316                             Zname{isim} = [list_of_packages(ifile).name,'_np_4_superludist_Mloop'];
0317                         end
0318                     end
0319                 end
0320             end
0321         end
0322     end           
0323     disp('---------------------------------------------------------------')  
0324 %
0325 % SUPERLU packages (MEX application)
0326 %
0327     if exist([local_root,'Packages/Bin/SUPERLU']) & SUPERLU,
0328         list_of_packages = dir([local_root,'Packages/Bin/SUPERLU']);
0329         list_of_packages = dir([local_root,'Packages/Bin/SUPERLU']);
0330         for ifile = 1:length(list_of_packages),   
0331             if length(list_of_packages(ifile).name) >= 7 
0332                 if list_of_packages(ifile).name(1:7) == 'SuperLU',
0333                     dkepath.SUPERLU.mex = [local_root,'Packages/Bin/SUPERLU/',list_of_packages(ifile).name,'/MATLAB',matver];
0334                     if exist(dkepath.SUPERLU.mex),
0335                         addpath(dkepath.SUPERLU.mex)
0336                         if exist('mexlusolve');isim = isim + 1;disp(['- ',list_of_packages(ifile).name,'-SOLVE (MEX)']),tic,[MMXS_t_SUPERLU] = lusolve(MMXPC_f_t,MMXR_t);t(isim) = toc;disp(['Elapsed time is ',num2str(t(isim)),' seconds.']);end
0337                         ZMMXS_t{isim} = MMXS_t_SUPERLU;
0338                         Zname{isim} = [list_of_packages(ifile).name,'_np_1_SUPERLUMEX'];  
0339                         rmpath(dkepath.SUPERLU.mex)   
0340                     end
0341                 end
0342             end
0343         end
0344     end
0345 %
0346     disp('---------------------------------------------------------------')  
0347 %
0348 % PARDISO packages (MEX application)
0349 %
0350     if exist([local_root,'Packages/Bin/PARDISO']) & PARDISO,
0351         list_of_packages = dir([local_root,'Packages/Bin/PARDISO']);
0352         for ifile = 1:length(list_of_packages),
0353             if strfind(list_of_packages(ifile).name,'MATLAB'),
0354                 dkepath.PARDISO.mex = [local_root,'Packages/Bin/PARDISO/MATLAB',matver];
0355                 if exist([dkepath.PARDISO.mex,filesep,'pardisosolve.',mexext]),
0356                     addpath(dkepath.PARDISO.mex);
0357                     addpath([local_root,'Packages/Bin/PARDISO']);
0358                     cd([local_root,'Packages/Bin/PARDISO']);
0359                     isim = isim + 1;disp(['- PARDISO-ILU (MEX) by CPU with sparse matrix']),tic,job_PARDISOMEX = 1;PARDISOMEX_dke_yp;job_PARDISOMEX = 2;PARDISOMEX_dke_yp;t(isim) = toc;
0360                     Zname{isim} = ['MEX_PARDISO_',list_of_packages(ifile).name,'_LU_FACTORISATION by CPU with sparse matrix'];
0361                     isim = isim + 1;disp(['- PARDISO-SOLVE (MEX) by CPU with sparse matrix']),tic,job_PARDISOMEX = 3;PARDISOMEX_dke_yp;t(isim) = toc;MMXS_t_PARDISO = X_PARDISOMEX;disp(['Elapsed time is ',num2str(t(isim)),' seconds per inversion + ILU time (once): ',num2str(t(isim-1)),' seconds.']);job_PARDISOMEX = -2;PARDISOMEX_dke_yp;
0362                     ZMMXS_t{isim} = MMXS_t_PARDISO;
0363                     Zname{isim} = ['MEX_PARDISO_',list_of_packages(ifile).name];  
0364                 end 
0365                 rmpath(dkepath.PARDISO.mex);
0366                 rmpath([local_root,'Packages/Bin/PARDISO']);
0367             end
0368         end
0369     end
0370     %
0371     ipass = ipass + 1;
0372     for ii = 1:isim
0373         tt(ipass,ii) = t(ii);
0374     end
0375 end
0376 %
0377 if mode==0, dirfig = 'figures/2D/';end
0378 if mode>=1, dirfig = 'figures/3D/';end
0379 %
0380 if npass > 1,
0381     meantt = mean(tt);%Mean value on several runs
0382     stdtt = std(tt);%Standard deviation on several runs
0383 %
0384     for iisim = 1:isim,
0385         figure('Name',Zname{iisim}),graph1D_jd([1:npass],tt(:,iisim),0,0,'Number of passes','Computing time (s)',[strrep(Zname{iisim},'_','-'),' on ',computer],[num2str(meantt(iisim)),'(s) +/- ',num2str(stdtt(iisim))],[1,npass],[0,2*max(tt(:,iisim))],'-','non','r',2,20,gca,0.85,0.7,0.5);hold on    
0386         savefig_dke_yp(Zname{iisim},dirfig);
0387     end
0388 end
0389 %
0390 for iisim = 1:isim,
0391     if ~isempty(ZMMXS_t{iisim})
0392         figure('Name',Zname{iisim}),graph1D_jd([1:length(ZMMXS_t{iisim})],ZMMXS_t{iisim},0,0,'Vector index','Vector coefficient',[strrep(Zname{iisim},'_','-'),' on ',computer],NaN,'','','-','non','r',1,20,gca,0.85,0.7,0.5);hold on    
0393         savefig_dke_yp(strrep(Zname{iisim},'.','_'),dirfig);
0394     end
0395 end
0396 %
0397 diary off
0398 %
0399 cd(dirinit);%back to the initial directory
0400 %