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 by Power System Engineering Research Center (PSERC)
0032 %   by Ray Zimmerman, PSERC Cornell
0033 %
0034 %   $Id: mplinsolve.m 2660 2015-03-20 02:10:18Z ray $
0035 %
0036 %   This file is part of MIPS.
0037 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0038 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0039 
0040 if nargin < 4
0041     opt = [];
0042     if nargin < 3
0043         solver = ''
0044     end
0045 end
0046 
0047 info = [];
0048 
0049 switch upper(solver)
0050     case {'', '\'}
0051         x = A \ b;
0052     case {'PARDISO'}
0053         %% get number of threads from OpenMP env variable
0054         persistent num_threads;
0055         if isempty(num_threads)
0056             num_threads = str2num(getenv('OMP_NUM_THREADS'));
0057             if ~num_threads
0058                 num_threads = 1;
0059             end
0060         end
0061 
0062         %% set default options
0063         verbose = false;
0064         mtype = 11;
0065         solver = 0;
0066         
0067         %% override if provided via opt
0068         if ~isempty(opt)
0069             if isfield(opt, 'verbose')
0070                 verbose = opt.verbose;
0071             end
0072             if isfield(opt, 'mtype')
0073                 mtype = opt.mtype;
0074             end
0075             if isfield(opt, 'solver')
0076                 solver = opt.solver;
0077             end
0078         end
0079 
0080         %% begin setup and solve
0081         info = pardisoinit(mtype, solver);
0082         info.iparm(3) = num_threads;
0083         if ~isempty(opt)
0084             if isfield(opt, 'iparm') && ~isempty(opt.iparm)
0085                 info.iparm(opt.iparm(:, 1)) = opt.iparm(:, 2);
0086             end
0087             if isfield(opt, 'dparm') && ~isempty(opt.dparm)
0088                 info.iparm(opt.dparm(:, 1)) = opt.dparm(:, 2);
0089             end
0090         end
0091         info = pardisoreorder(A, info, verbose);
0092         info = pardisofactor(A, info, verbose);
0093         [x, info] = pardisosolve(A, b, info, verbose);
0094         pardisofree(info);
0095     otherwise
0096         warning('mplinsolve: ''%s'' is not a valid value for SOLVER, using default.', solver);
0097         x = A \ b;
0098 end