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 if issparse(A) 0074 [L, U, p] = lu(A(q,q), 1.0, 'vector'); 0075 else 0076 [L, U, p] = lu(A(q,q), 'vector'); 0077 end 0078 x = zeros(size(A, 1), 1); 0079 x(q) = U \ ( L \ b(q(p)) ); 0080 case 'LU3a' %% 3 output LU: Gilbert-Peierls alg, permutation vectors 0081 q = amd(A); %% permutation vector for AMD reordering 0082 [L, U, p] = lu(A(q,q), 'vector'); 0083 x = zeros(size(A, 1), 1); 0084 x(q) = U \ ( L \ b(q(p)) ); 0085 case 'LU4' %% 4 output LU: UMFPACK, permutation vectors 0086 [L, U, p, q] = lu(A, 'vector'); 0087 x = zeros(size(A, 1), 1); 0088 x(q) = U \ ( L \ b(p)); 0089 case 'LU5' %% 5 output LU: UMFPACK w/row scaling, permutation vectors 0090 [L, U, p, q, R] = lu(A, 'vector'); 0091 x = zeros(size(A, 1), 1); 0092 x(q) = U \ ( L \ (R(:, p) \ b)); 0093 case 'LU3m' %% 3 output LU: Gilbert-Peierls alg, perm mat, 1.0 piv thresh 0094 Q = sparse(amd(A), 1:size(A, 1), 1); %% permutation matrix for AMD reordering 0095 if issparse(A) 0096 [L, U, P] = lu(Q'*A*Q, 1.0); 0097 else 0098 [L, U, P] = lu(Q'*A*Q); 0099 end 0100 x = Q * ( U \ ( L \ (P * Q' * b)) ); 0101 case 'LU3am' %% 3 output LU: Gilbert-Peierls alg, permutation matrices 0102 Q = sparse(amd(A), 1:size(A, 1), 1); %% permutation matrix for AMD reordering 0103 [L, U, P] = lu(Q'*A*Q); 0104 x = Q * ( U \ ( L \ (P * Q' * b)) ); 0105 case 'LU4m' %% 4 output LU: UMFPACK, permutation matrices 0106 [L, U, P, Q] = lu(A); 0107 x = Q * (U \ (L \ (P * b))); 0108 case 'LU5m' %% 5 output LU: UMFPACK w/row scaling, permutation matrices 0109 [L, U, P, Q, R] = lu(A); 0110 x = Q * (U \ (L \ (P * (R \ b)))); 0111 case 'LU' %% explicit LU, with options struct 0112 %% default options 0113 nout = 4; %% 4 output args, UMFPACK 0114 if ~issparse(A) 0115 nout = 3; 0116 end 0117 vec = have_feature('lu_vec'); %% use permulation vectors, if available 0118 thresh = []; %% use default pivot threshold 0119 if isfield(opt, 'lu') 0120 opt_lu = opt.lu; 0121 if isfield(opt_lu, 'nout') 0122 nout = opt_lu.nout; 0123 end 0124 if isfield(opt_lu, 'vec') 0125 vec = opt_lu.vec; 0126 end 0127 if isfield(opt_lu, 'thresh') 0128 thresh = opt_lu.thresh; 0129 end 0130 end 0131 %% call the appropriate form 0132 switch nout 0133 case 3 %% 3 output args: Gilbert-Peierls algorithm, with AMD reordering 0134 q = amd(A); %% permutation vector for AMD reordering 0135 n = size(A, 1); 0136 if vec 0137 if isempty(thresh) 0138 [L, U, p] = lu(A(q,q), 'vector'); 0139 else 0140 [L, U, p] = lu(A(q,q), thresh, 'vector'); 0141 end 0142 x = zeros(n, 1); 0143 x(q) = U \ ( L \ b(q(p)) ); 0144 else 0145 Q = sparse(q, 1:n, 1); %% permutation matrix for AMD reordering 0146 if isempty(thresh) 0147 [L, U, P] = lu(Q'*A*Q); 0148 else 0149 [L, U, P] = lu(Q'*A*Q, thresh); 0150 end 0151 x = Q * ( U \ ( L \ (P * Q' * b)) ); 0152 end 0153 case 4 %% 4 output args: UMFPACK 0154 if vec 0155 [L, U, p, q] = lu(A, 'vector'); 0156 x = zeros(size(A, 1), 1); 0157 x(q) = U \ ( L \ b(p)); 0158 else 0159 [L, U, P, Q] = lu(A); 0160 x = Q * (U \ (L \ (P * b))); 0161 end 0162 case 5 %% 5 output args: UMFPACK w/row scaling 0163 if vec 0164 [L, U, p, q, R] = lu(A, 'vector'); 0165 x = zeros(size(A, 1), 1); 0166 x(q) = U \ ( L \ (R(:, p) \ b)); 0167 else 0168 [L, U, P, Q, R] = lu(A); 0169 x = Q * (U \ (L \ (P * (R \ b)))); 0170 end 0171 end 0172 case {'PARDISO'} 0173 %% set default options 0174 verbose = false; 0175 mtype = 11; 0176 pardiso_solver = 0; 0177 0178 %% override if provided via opt 0179 if ~isempty(opt) && isfield(opt, 'pardiso') 0180 if isfield(opt.pardiso, 'verbose') && opt.pardiso.verbose 0181 verbose = true; 0182 end 0183 if isfield(opt.pardiso, 'mtype') 0184 mtype = opt.pardiso.mtype; 0185 end 0186 if isfield(opt.pardiso, 'solver') 0187 pardiso_solver = opt.pardiso.solver; 0188 end 0189 end 0190 0191 %% begin setup and solve 0192 v6 = have_feature('pardiso_object'); 0193 if v6 %% PARDISO v6+ 0194 id = 1; 0195 p = pardiso(id, mtype, pardiso_solver); 0196 if verbose 0197 p.verbose(); 0198 end 0199 else %% PARDISO v5 0200 p = pardisoinit(mtype, pardiso_solver); 0201 end 0202 if ~isempty(opt) && isfield(opt, 'pardiso') 0203 if isfield(opt.pardiso, 'iparm') && ~isempty(opt.pardiso.iparm) 0204 p.iparm(opt.pardiso.iparm(:, 1)) = opt.pardiso.iparm(:, 2); 0205 end 0206 if isfield(opt.pardiso, 'dparm') && ~isempty(opt.pardiso.dparm) 0207 p.dparm(opt.pardiso.dparm(:, 1)) = opt.pardiso.dparm(:, 2); 0208 end 0209 end 0210 if v6 || abs(mtype) == 2 || mtype == 6 %% need non-zero diagonal 0211 nx = size(A, 1); 0212 if abs(mtype) == 2 || mtype == 6 %% symmetric 0213 myeps = 1e-14; 0214 A = tril(A); 0215 else %% non-symmetric 0216 myeps = 1e-8; 0217 end 0218 A = A + myeps * speye(nx, nx); 0219 end 0220 if v6 0221 p.factorize(id, A); 0222 x = p.solve(id, A, b); 0223 p.free(id); 0224 p.clear(); 0225 else 0226 p = pardisoreorder(A, p, verbose); 0227 p = pardisofactor(A, p, verbose); 0228 [x, p] = pardisosolve(A, b, p, verbose); 0229 pardisofree(p); 0230 end 0231 otherwise 0232 warning('mplinsolve: ''%s'' is not a valid value for SOLVER, using default.', solver); 0233 x = A \ b; 0234 end 0235 0236 0237 function num = vstr2num_(vstr) 0238 % Converts version string to numerical value suitable for < or > comparisons 0239 % E.g. '3.11.4' --> 3.011004 0240 pat = '\.?(\d+)'; 0241 [s,e,tE,m,t] = regexp(vstr, pat); 0242 b = 1; 0243 num = 0; 0244 for k = 1:length(t) 0245 num = num + b * str2num(t{k}{1}); 0246 b = b / 1000; 0247 end