Home > matpower7.0 > mips > lib > mplinsolve.m

mplinsolve

PURPOSE ^

MPLINSOLVE Solves A * x = b using specified solver

SYNOPSIS ^

function [x, info] = mplinsolve(A, b, solver, opt)

DESCRIPTION ^

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 respectively

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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 %           '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

Generated on Mon 24-Jun-2019 15:58:45 by m2html © 2005