Home > matpower7.1 > 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         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

Generated on Fri 09-Oct-2020 11:21:31 by m2html © 2005