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
'LU' - use explicit LU decomposition and back substitution
The following are also provided as short-cuts, with less
overhead and thus better performance on small systems, as
alternatives for 'LU' with the following options:
OPT.lu.nout OPT.lu.vec OPT.lu.thresh
'LU3' 3 1 1.0
'LU3a' 3 1
'LU4' 4 1
'LU5' 5 1
'LU3m' 3 0 1.0
'LU3am' 3 0
'LU4m' 4 0
'LU5m' 5 0
'PARDISO' - PARDISO
OPT : struct of options for certain solvers (e.g. LU and PARDISO)
lu : struct of options to determine form of call to LU solver,
with the following possible fields (default value in parens):
nout (4) - number of output args for call to LU, UMFPACK is
used for 4 or 5 output args, Gilbert-Peierls algorithm
with AMD ordering for 3 output args.
vec (1) - use permutation vectors instead of matrices
(permutation matrices used by default for MATLAB < 7.3)
thresh - pivot threshold, see 'help lu' for details
pardiso : struct of PARDISO options (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 respectively0001 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 % 'LU' - use explicit LU decomposition and back substitution 0019 % The following are also provided as short-cuts, with less 0020 % overhead and thus better performance on small systems, as 0021 % alternatives for 'LU' with the following options: 0022 % OPT.lu.nout OPT.lu.vec OPT.lu.thresh 0023 % 'LU3' 3 1 1.0 0024 % 'LU3a' 3 1 0025 % 'LU4' 4 1 0026 % 'LU5' 5 1 0027 % 'LU3m' 3 0 1.0 0028 % 'LU3am' 3 0 0029 % 'LU4m' 4 0 0030 % 'LU5m' 5 0 0031 % 'PARDISO' - PARDISO 0032 % OPT : struct of options for certain solvers (e.g. LU and PARDISO) 0033 % lu : struct of options to determine form of call to LU solver, 0034 % with the following possible fields (default value in parens): 0035 % nout (4) - number of output args for call to LU, UMFPACK is 0036 % used for 4 or 5 output args, Gilbert-Peierls algorithm 0037 % with AMD ordering for 3 output args. 0038 % vec (1) - use permutation vectors instead of matrices 0039 % (permutation matrices used by default for MATLAB < 7.3) 0040 % thresh - pivot threshold, see 'help lu' for details 0041 % pardiso : struct of PARDISO options (default shown in parens), 0042 % see PARDISO documentation for details 0043 % verbose (0) - true or false 0044 % mtype (11) - matrix type (default is real and nonsymmetric) 0045 % solver (0) - solver method (default is sparse direct) 0046 % iparm ([]) - n x 2 matrix of integer parameters 0047 % 1st, 2nd columns are index, value of parameter respectively 0048 % dparm ([]) - n x 2 matrix of double parameters 0049 % 1st, 2nd columns are index, value of parameter respectively 0050 0051 % MIPS 0052 % Copyright (c) 2015-2018, Power Systems Engineering Research Center (PSERC) 0053 % by Ray Zimmerman, PSERC Cornell 0054 % 0055 % This file is part of MIPS. 0056 % Covered by the 3-clause BSD License (see LICENSE file for details). 0057 % See https://github.com/MATPOWER/mips for more info. 0058 0059 if nargin < 4 0060 opt = []; 0061 if nargin < 3 0062 solver = '' 0063 end 0064 end 0065 0066 info = []; 0067 0068 switch solver 0069 case {'', '\'} %% use built-in \ operator 0070 x = A \ b; 0071 case 'LU3' %% 3 output LU: Gilbert-Peierls alg, perm vec, 1.0 piv thresh 0072 q = amd(A); %% permutation vector for AMD reordering 0073 [L, U, p] = lu(A(q,q), 1.0, 'vector'); 0074 x = zeros(size(A, 1), 1); 0075 x(q) = U \ ( L \ b(q(p)) ); 0076 case 'LU3a' %% 3 output LU: Gilbert-Peierls alg, permutation vectors 0077 q = amd(A); %% permutation vector for AMD reordering 0078 [L, U, p] = lu(A(q,q), 'vector'); 0079 x = zeros(size(A, 1), 1); 0080 x(q) = U \ ( L \ b(q(p)) ); 0081 case 'LU4' %% 4 output LU: UMFPACK, permutation vectors 0082 [L, U, p, q] = lu(A, 'vector'); 0083 x = zeros(size(A, 1), 1); 0084 x(q) = U \ ( L \ b(p)); 0085 case 'LU5' %% 5 output LU: UMFPACK w/row scaling, permutation vectors 0086 [L, U, p, q, R] = lu(A, 'vector'); 0087 x = zeros(size(A, 1), 1); 0088 x(q) = U \ ( L \ (R(:, p) \ b)); 0089 case 'LU3m' %% 3 output LU: Gilbert-Peierls alg, perm mat, 1.0 piv thresh 0090 Q = sparse(amd(A), 1:size(A, 1), 1); %% permutation matrix for AMD reordering 0091 [L, U, P] = lu(Q'*A*Q, 1.0); 0092 x = Q * ( U \ ( L \ (P * Q' * b)) ); 0093 case 'LU3am' %% 3 output LU: Gilbert-Peierls alg, permutation matrices 0094 Q = sparse(amd(A), 1:size(A, 1), 1); %% permutation matrix for AMD reordering 0095 [L, U, P] = lu(Q'*A*Q); 0096 x = Q * ( U \ ( L \ (P * Q' * b)) ); 0097 case 'LU4m' %% 4 output LU: UMFPACK, permutation matrices 0098 [L, U, P, Q] = lu(A); 0099 x = Q * (U \ (L \ (P * b))); 0100 case 'LU5m' %% 5 output LU: UMFPACK w/row scaling, permutation matrices 0101 [L, U, P, Q, R] = lu(A); 0102 x = Q * (U \ (L \ (P * (R \ b)))); 0103 case 'LU' %% explicit LU, with options struct 0104 %% default options 0105 nout = 4; %% 4 output args, UMFPACK 0106 vec = have_lu_vec(); %% use permulation vectors, if available 0107 thresh = []; %% use default pivot threshold 0108 if isfield(opt, 'lu') 0109 opt_lu = opt.lu; 0110 if isfield(opt_lu, 'nout') 0111 nout = opt_lu.nout; 0112 end 0113 if isfield(opt_lu, 'vec') 0114 vec = opt_lu.vec; 0115 end 0116 if isfield(opt_lu, 'thresh') 0117 thresh = opt_lu.thresh; 0118 end 0119 end 0120 %% call the appropriate form 0121 switch nout 0122 case 3 %% 3 output args: Gilbert-Peierls algorithm, with AMD reordering 0123 q = amd(A); %% permutation vector for AMD reordering 0124 n = size(A, 1); 0125 if vec 0126 if isempty(thresh) 0127 [L, U, p] = lu(A(q,q), 'vector'); 0128 else 0129 [L, U, p] = lu(A(q,q), thresh, 'vector'); 0130 end 0131 x = zeros(n, 1); 0132 x(q) = U \ ( L \ b(q(p)) ); 0133 else 0134 Q = sparse(q, 1:n, 1); %% permutation matrix for AMD reordering 0135 if isempty(thresh) 0136 [L, U, P] = lu(Q'*A*Q); 0137 else 0138 [L, U, P] = lu(Q'*A*Q, thresh); 0139 end 0140 x = Q * ( U \ ( L \ (P * Q' * b)) ); 0141 end 0142 case 4 %% 4 output args: UMFPACK 0143 if vec 0144 [L, U, p, q] = lu(A, 'vector'); 0145 x = zeros(size(A, 1), 1); 0146 x(q) = U \ ( L \ b(p)); 0147 else 0148 [L, U, P, Q] = lu(A); 0149 x = Q * (U \ (L \ (P * b))); 0150 end 0151 case 5 %% 5 output args: UMFPACK w/row scaling 0152 if vec 0153 [L, U, p, q, R] = lu(A, 'vector'); 0154 x = zeros(size(A, 1), 1); 0155 x(q) = U \ ( L \ (R(:, p) \ b)); 0156 else 0157 [L, U, P, Q, R] = lu(A); 0158 x = Q * (U \ (L \ (P * (R \ b)))); 0159 end 0160 end 0161 case {'PARDISO'} 0162 %% set default options 0163 verbose = false; 0164 mtype = 11; 0165 pardiso_solver = 0; 0166 0167 %% override if provided via opt 0168 if ~isempty(opt) && isfield(opt, 'pardiso') 0169 if isfield(opt.pardiso, 'verbose') && opt.pardiso.verbose 0170 verbose = true; 0171 end 0172 if isfield(opt.pardiso, 'mtype') 0173 mtype = opt.pardiso.mtype; 0174 end 0175 if isfield(opt.pardiso, 'solver') 0176 pardiso_solver = opt.pardiso.solver; 0177 end 0178 end 0179 0180 %% begin setup and solve 0181 v6 = have_pardiso_object(); 0182 if v6 %% PARDISO v6+ 0183 id = 1; 0184 p = pardiso(id, mtype, pardiso_solver); 0185 if verbose 0186 p.verbose(); 0187 end 0188 else %% PARDISO v5 0189 p = pardisoinit(mtype, pardiso_solver); 0190 end 0191 if ~isempty(opt) && isfield(opt, 'pardiso') 0192 if isfield(opt.pardiso, 'iparm') && ~isempty(opt.pardiso.iparm) 0193 p.iparm(opt.pardiso.iparm(:, 1)) = opt.pardiso.iparm(:, 2); 0194 end 0195 if isfield(opt.pardiso, 'dparm') && ~isempty(opt.pardiso.dparm) 0196 p.dparm(opt.pardiso.dparm(:, 1)) = opt.pardiso.dparm(:, 2); 0197 end 0198 end 0199 if v6 || abs(mtype) == 2 || mtype == 6 %% need non-zero diagonal 0200 nx = size(A, 1); 0201 if abs(mtype) == 2 || mtype == 6 %% symmetric 0202 myeps = 1e-14; 0203 A = tril(A); 0204 else %% non-symmetric 0205 myeps = 1e-8; 0206 end 0207 A = A + myeps * speye(nx, nx); 0208 end 0209 if v6 0210 p.factorize(id, A); 0211 x = p.solve(id, A, b); 0212 p.free(id); 0213 p.clear(); 0214 else 0215 p = pardisoreorder(A, p, verbose); 0216 p = pardisofactor(A, p, verbose); 0217 [x, p] = pardisosolve(A, b, p, verbose); 0218 pardisofree(p); 0219 end 0220 otherwise 0221 warning('mplinsolve: ''%s'' is not a valid value for SOLVER, using default.', solver); 0222 x = A \ b; 0223 end 0224 0225 0226 function TorF = have_lu_vec() 0227 % Checks whether or not LU supports lu(..., 'vector') syntax 0228 persistent lu_vec; %% cache the result for performance reasons 0229 if isempty(lu_vec) 0230 lu_vec = 1; %% assume it does, unless this is MATLAB ver < 7.3 0231 v = ver('matlab'); 0232 if length(v) > 1 0233 warning('The built-in VER command is behaving strangely, probably as a result of installing a 3rd party toolbox in a directory named ''matlab'' on your path. Check each element of the output of ver(''matlab'') to find the offending toolbox, then move the toolbox to a more appropriately named directory.'); 0234 v = v(1); 0235 end 0236 if ~isempty(v) && isfield(v, 'Version') && ~isempty(v.Version) 0237 vstr = v.Version; 0238 if ~isempty(vstr) && vstr2num_(vstr) < 7.003 0239 lu_vec = 0; 0240 end 0241 end 0242 end 0243 TorF = lu_vec; 0244 0245 0246 function TorF = have_pardiso_object() 0247 % Checks for availability of PARDISO 6 0248 persistent pardiso_object; %% cache the result for performance reasons 0249 if isempty(pardiso_object) 0250 pardiso_object = exist('pardiso', 'file') == 2; 0251 end 0252 TorF = pardiso_object; 0253 0254 0255 function num = vstr2num_(vstr) 0256 % Converts version string to numerical value suitable for < or > comparisons 0257 % E.g. '3.11.4' --> 3.011004 0258 pat = '\.?(\d+)'; 0259 [s,e,tE,m,t] = regexp(vstr, pat); 0260 b = 1; 0261 num = 0; 0262 for k = 1:length(t) 0263 num = num + b * str2num(t{k}{1}); 0264 b = b / 1000; 0265 end