MUMPSMEX_dke_yp

PURPOSE ^

Link with the parallel matrix solver package MUMPS written in FORTAN 95

SYNOPSIS ^

This is a script file.

DESCRIPTION ^

Link with the parallel matrix solver package MUMPS written in FORTAN 95
for MEX file (single processor).  

Solve the linear system of equation AX=B

 INPUTS:

    - A: Matrix [n,n]
   - B: Vector for solving the linear system of equations AX=B [n,1]  
   - job: MUMPS job type (1,2,3,5,6,-2) see MUMPS documentation
  
 OUTPUTS:

   - X: Solution of the linear system of equation [n,1]
   - flag: execution flag (see the MUMPS documentation) [1,1]

By Joan Decker (CEA-DRFC, joan.decker@cea.fr) and Yves Peysson (CEA-DRFC, yves.peysson@cea.fr)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 %
0002 %Link with the parallel matrix solver package MUMPS written in FORTAN 95
0003 %for MEX file (single processor).
0004 %
0005 %Solve the linear system of equation AX=B
0006 %
0007 % INPUTS:
0008 %
0009 %    - A: Matrix [n,n]
0010 %   - B: Vector for solving the linear system of equations AX=B [n,1]
0011 %   - job: MUMPS job type (1,2,3,5,6,-2) see MUMPS documentation
0012 %
0013 % OUTPUTS:
0014 %
0015 %   - X: Solution of the linear system of equation [n,1]
0016 %   - flag: execution flag (see the MUMPS documentation) [1,1]
0017 %
0018 %By Joan Decker (CEA-DRFC, joan.decker@cea.fr) and Yves Peysson (CEA-DRFC, yves.peysson@cea.fr)
0019 %
0020 if ~isfield(dkeparam,'opsplit'),dkeparam.opsplit = 0;end
0021 %
0022 if dkeparam.opsplit == 0,
0023     A_MUMPSMEX = MMXPC_f_t;
0024     B_MUMPSMEX = MMXR_t;
0025     %
0026     if job_MUMPSMEX == 1,
0027         id_MUMPSMEX = initmumps;
0028         id_MUMPSMEX = dmumps(id_MUMPSMEX);
0029         id_MUMPSMEX.ICNTL(1:4) = [6,0,0,1];%no display during runtime
0030         id_MUMPSMEX.JOB = job_MUMPSMEX;
0031         id_MUMPSMEX = dmumps(id_MUMPSMEX,A_MUMPSMEX);% Call to mumps
0032         if isfield(dkeparam.MUMPSparam,'thresholdpivoting'),
0033             if ~isempty(dkeparam.MUMPSparam.thresholdpivoting) | ~isnan(dkeparam.MUMPSparam.thresholdpivoting),
0034                 id_MUMPSMEX.CNTL(1) = dkeparam.MUMPSparam.thresholdpivoting;%Relative threshold for numerical pivoting in the MUMPS partial LU decomposition (default MUMPS = 1e-2 for non-symmetric matrices)
0035             end
0036         end
0037         if isfield(dkeparam.MUMPSparam,'precinv'),
0038             if ~isempty(dkeparam.MUMPSparam.precinv) | ~isnan(dkeparam.MUMPSparam.precinv),
0039                 id_MUMPSMEX.CNTL(2) = dkeparam.MUMPSparam.precinv;%Drop tolerance for the MUMPS partial LU decomposition
0040             end
0041         end    
0042     elseif job_MUMPSMEX == 2,
0043        id_MUMPSMEX.JOB = job_MUMPSMEX;
0044        id_MUMPSMEX = dmumps(id_MUMPSMEX,A_MUMPSMEX);% Call to mumps
0045        if id_MUMPSMEX.INFOG(1) == -13,
0046             error('Error in a Fortran ALLOCATE statement in MUMPS factorization. Try reducing the problem size (grid size)');
0047        end
0048     elseif job_MUMPSMEX == 3 | job_MUMPSMEX == 6,
0049        id_MUMPSMEX.JOB = job_MUMPSMEX;
0050        id_MUMPSMEX.RHS = full(B_MUMPSMEX);
0051        id_MUMPSMEX = dmumps(id_MUMPSMEX,A_MUMPSMEX);% Call to mumps
0052        if(norm(A_MUMPSMEX*id_MUMPSMEX.SOL - id_MUMPSMEX.RHS,'inf') > sqrt(eps))
0053           warning('Precision of the matrix inversion may not be OK (MUMPSMEX). Try to increase the dkeparam.MUMPSparam.thresholdpivoting parameters (see MUMPS documentation).');
0054        end
0055        X_MUMPSMEX = id_MUMPSMEX.SOL;
0056        flag_MUMPSMEX = 0;
0057     elseif job_MUMPSMEX == -2,
0058        id_MUMPSMEX.JOB = job_MUMPSMEX;
0059        id_MUMPSMEX = dmumps(id_MUMPSMEX);% destroy mumps instance
0060     end
0061 elseif dkeparam.opsplit >= 1,
0062     if mode_MUMPSMEX,
0063         A_MUMPSMEX = MMXPC_P_f_t;
0064         B_MUMPSMEX = MMXR_t;
0065         %
0066         if job_MUMPSMEX == 1,
0067             id_MUMPSMEX_P = initmumps;
0068             id_MUMPSMEX_P = dmumps(id_MUMPSMEX_P);
0069             id_MUMPSMEX_P.ICNTL(1:4) = [6,0,0,1];%no display during runtime
0070             id_MUMPSMEX_P.JOB = job_MUMPSMEX;
0071             id_MUMPSMEX_P = dmumps(id_MUMPSMEX_P,A_MUMPSMEX);% Call to mumps
0072             if isfield(dkeparam.MUMPSparam,'thresholdpivoting'),
0073                 if ~isempty(dkeparam.MUMPSparam.thresholdpivoting) | ~isnan(dkeparam.MUMPSparam.thresholdpivoting),
0074                     id_MUMPSMEX_P.CNTL(1) = dkeparam.MUMPSparam.thresholdpivoting;%Relative threshold for numerical pivoting in the MUMPS partial LU decomposition (default MUMPS = 1e-2 for non-symmetric matrices)
0075                 end
0076             end
0077             if isfield(dkeparam.MUMPSparam,'precinv'),
0078                 if ~isempty(dkeparam.MUMPSparam.precinv) | ~isnan(dkeparam.MUMPSparam.precinv),
0079                     id_MUMPSMEX_P.CNTL(2) = dkeparam.MUMPSparam.precinv;%Drop tolerance for the MUMPS partial LU decomposition
0080                 end
0081             end    
0082         elseif job_MUMPSMEX == 2,
0083            id_MUMPSMEX_P.JOB = job_MUMPSMEX;
0084            id_MUMPSMEX_P = dmumps(id_MUMPSMEX_P,A_MUMPSMEX);% Call to mumps
0085            if id_MUMPSMEX_P.INFOG(1) == -13,
0086                 error('Error in a Fortran ALLOCATE statement in MUMPS factorization. Try reducing the problem size (grid size)');
0087            end
0088         elseif job_MUMPSMEX == 3 | job_MUMPSMEX == 6,
0089            id_MUMPSMEX_P.JOB = job_MUMPSMEX;
0090            id_MUMPSMEX_P.RHS = full(B_MUMPSMEX);
0091            id_MUMPSMEX_P = dmumps(id_MUMPSMEX_P,A_MUMPSMEX);% Call to mumps
0092            if(norm(A_MUMPSMEX*id_MUMPSMEX_P.SOL - id_MUMPSMEX_P.RHS,'inf') > sqrt(eps))
0093               warning('Precision of the matrix inversion may not be OK (MUMPSMEX). Try to increase the dkeparam.MUMPSparam.thresholdpivoting parameters (see MUMPS documentation).');
0094            end
0095            X_MUMPSMEX = id_MUMPSMEX_P.SOL;
0096            flag_MUMPSMEX = 0;
0097         elseif job_MUMPSMEX == -2,
0098            id_MUMPSMEX_P.JOB = job_MUMPSMEX;
0099            id_MUMPSMEX_P = dmumps(id_MUMPSMEX_P);% destroy mumps instance
0100         end 
0101     else
0102         A_MUMPSMEX = MMXPC_R_f_t;
0103         B_MUMPSMEX = MMXR_t;
0104         %
0105         if job_MUMPSMEX == 1,
0106             id_MUMPSMEX_R = initmumps;
0107             id_MUMPSMEX_R = dmumps(id_MUMPSMEX_R);
0108             id_MUMPSMEX_R.ICNTL(1:4) = [6,0,0,1];%no display during runtime
0109             id_MUMPSMEX_R.JOB = job_MUMPSMEX;
0110             id_MUMPSMEX_R = dmumps(id_MUMPSMEX_R,A_MUMPSMEX);% Call to mumps
0111             if isfield(dkeparam.MUMPSparam,'thresholdpivoting'),
0112                 if ~isempty(dkeparam.MUMPSparam.thresholdpivoting) | ~isnan(dkeparam.MUMPSparam.thresholdpivoting),
0113                     id_MUMPSMEX_R.CNTL(1) = dkeparam.MUMPSparam.thresholdpivoting;%Relative threshold for numerical pivoting in the MUMPS partial LU decomposition (default MUMPS = 1e-2 for non-symmetric matrices)
0114                 end
0115             end
0116             if isfield(dkeparam.MUMPSparam,'precinv'),
0117                 if ~isempty(dkeparam.MUMPSparam.precinv) | ~isnan(dkeparam.MUMPSparam.precinv),
0118                     id_MUMPSMEX_R.CNTL(2) = dkeparam.MUMPSparam.precinv;%Drop tolerance for the MUMPS partial LU decomposition
0119                 end
0120             end    
0121         elseif job_MUMPSMEX == 2,
0122            id_MUMPSMEX_R.JOB = job_MUMPSMEX;
0123            id_MUMPSMEX_R = dmumps(id_MUMPSMEX_R,A_MUMPSMEX);% Call to mumps
0124            if id_MUMPSMEX_R.INFOG(1) == -13,
0125                 error('Error in a Fortran ALLOCATE statement in MUMPS factorization. Try reducing the problem size (grid size)');
0126            end
0127         elseif job_MUMPSMEX == 3 | job_MUMPSMEX == 6,
0128            id_MUMPSMEX_R.JOB = job_MUMPSMEX;
0129            id_MUMPSMEX_R.RHS = full(B_MUMPSMEX);
0130            id_MUMPSMEX_R = dmumps(id_MUMPSMEX_R,A_MUMPSMEX);% Call to mumps
0131            if(norm(A_MUMPSMEX*id_MUMPSMEX_R.SOL - id_MUMPSMEX_R.RHS,'inf') > sqrt(eps))
0132               warning('Precision of the matrix inversion may not be OK (MUMPSMEX). Try to increase the dkeparam.MUMPSparam.thresholdpivoting parameters (see MUMPS documentation).');
0133            end
0134            X_MUMPSMEX = id_MUMPSMEX_R.SOL;
0135            flag_MUMPSMEX = 0;
0136         elseif job_MUMPSMEX == -2,
0137            id_MUMPSMEX_R.JOB = job_MUMPSMEX;
0138            id_MUMPSMEX_R = dmumps(id_MUMPSMEX_R);% destroy mumps instance
0139         end 
0140     end
0141 end
0142   
0143   
0144
Community support and wiki are available on Redmine. Last update: 18-Apr-2019.