mplinsolve

MPLINSOLVE  Solves A * x = b using specified solver
   X = MPLINSOLVE(A, B)
   X = MPLINSOLVE(A, B, SOLVER)
   X = MPLINSOLVE(A, B, SOLVER, OPT)
   [X, INFO] = MPLINSOLVE(...)

   Solves the linear system of equations A * x = b, using the selected
   solver.

   Inputs:
       A      : sparse matrix
       B      : full vector or matrix
       SOLVER : ('') selected linear system solver
           ''  - use default solver, currently this is
                   always the built-in backslash operator
           '\' - built-in backslash operator
           'PARDISO' - PARDISO
       OPT    : struct of options, with the following fields
                (currently used only by PARDISO, default shown in parens,
                 see PARDISO documentation for details)
           verbose (0) - true or false
           mtype (11)  - matrix type (default is real and nonsymmetric)
           solver (0)  - solver method (default is sparse direct)
           iparm ([])  - n x 2 matrix of integer parameters
               1st, 2nd columns are index, value of parameter respectively
           dparm ([])  - n x 2 matrix of double parameters
               1st, 2nd columns are index, value of parameter respectively

0001 function [x, info] = mplinsolve(A, b, solver, opt)
0002 %MPLINSOLVE  Solves A * x = b using specified solver
0003 %   X = MPLINSOLVE(A, B)
0004 %   X = MPLINSOLVE(A, B, SOLVER)
0005 %   X = MPLINSOLVE(A, B, SOLVER, OPT)
0006 %   [X, INFO] = MPLINSOLVE(...)
0007 %
0008 %   Solves the linear system of equations A * x = b, using the selected
0009 %   solver.
0010 %
0011 %   Inputs:
0012 %       A      : sparse matrix
0013 %       B      : full vector or matrix
0014 %       SOLVER : ('') selected linear system solver
0015 %           ''  - use default solver, currently this is
0016 %                   always the built-in backslash operator
0017 %           '\' - built-in backslash operator
0018 %           'PARDISO' - PARDISO
0019 %       OPT    : struct of options, with the following fields
0020 %                (currently used only by PARDISO, default shown in parens,
0021 %                 see PARDISO documentation for details)
0022 %           verbose (0) - true or false
0023 %           mtype (11)  - matrix type (default is real and nonsymmetric)
0024 %           solver (0)  - solver method (default is sparse direct)
0025 %           iparm ([])  - n x 2 matrix of integer parameters
0026 %               1st, 2nd columns are index, value of parameter respectively
0027 %           dparm ([])  - n x 2 matrix of double parameters
0028 %               1st, 2nd columns are index, value of parameter respectively
0029 
0030 %   MIPS
0031 %   Copyright (c) 2015-2016, Power Systems Engineering Research Center (PSERC)
0032 %   by Ray Zimmerman, PSERC Cornell
0033 %
0034 %   This file is part of MIPS.
0035 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0036 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0037 
0038 if nargin < 4
0039     opt = [];
0040     if nargin < 3
0041         solver = ''
0042     end
0043 end
0044 
0045 info = [];
0046 
0047 switch upper(solver)
0048     case {'', '\'}
0049         x = A \ b;
0050     case {'PARDISO'}
0051         %% get number of threads from OpenMP env variable
0052         persistent num_threads;
0053         if isempty(num_threads)
0054             num_threads = str2num(getenv('OMP_NUM_THREADS'));
0055             if ~num_threads
0056                 num_threads = 1;
0057             end
0058         end
0059 
0060         %% set default options
0061         verbose = false;
0062         mtype = 11;
0063         solver = 0;
0064         
0065         %% override if provided via opt
0066         if ~isempty(opt)
0067             if isfield(opt, 'verbose')
0068                 verbose = opt.verbose;
0069             end
0070             if isfield(opt, 'mtype')
0071                 mtype = opt.mtype;
0072             end
0073             if isfield(opt, 'solver')
0074                 solver = opt.solver;
0075             end
0076         end
0077 
0078         %% begin setup and solve
0079         info = pardisoinit(mtype, solver);
0080         info.iparm(3) = num_threads;
0081         if ~isempty(opt)
0082             if isfield(opt, 'iparm') && ~isempty(opt.iparm)
0083                 info.iparm(opt.iparm(:, 1)) = opt.iparm(:, 2);
0084             end
0085             if isfield(opt, 'dparm') && ~isempty(opt.dparm)
0086                 info.iparm(opt.dparm(:, 1)) = opt.dparm(:, 2);
0087             end
0088         end
0089         info = pardisoreorder(A, info, verbose);
0090         info = pardisofactor(A, info, verbose);
0091         [x, info] = pardisosolve(A, b, info, verbose);
0092         pardisofree(info);
0093     otherwise
0094         warning('mplinsolve: ''%s'' is not a valid value for SOLVER, using default.', solver);
0095         x = A \ b;
0096 end