1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64
|
function [p,q,r,s,cc,rr] = cs_dmperm (A,seed) %#ok
%CS_DMPERM maximum matching or Dulmage-Mendelsohn permutation.
% p = cs_dmperm(A) finds a maximum matching p such that p(j) = i if column j
% is matched to row i, or 0 if column j is unmatched. If A is square and
% full structural rank, p is a row permutation and A(p,:) has a zero-free
% diagonal. The structural rank of A is sprank(A) = sum(p>0).
%
% [p,q,r,s,cc,rr] = cs_dmperm(A) finds the Dulmage-Mendelsohn decomposition
% of A. p and q are permutation vectors. cc and rr are vectors of length 5.
% C = A(p,q) is split into a 4-by-4 set of coarse blocks:
%
% A11 A12 A13 A14
% 0 0 A23 A24
% 0 0 0 A34
% 0 0 0 A44
%
% where A12, A23, and A34 are square with zero-free diagonals. The columns of
% A11 are the unmatched columns, and the rows of A44 are the unmatched rows.
% Any of these blocks can be empty. In the "coarse" decomposition, the
% (i,j)th block is C(rr(i):rr(i+1)-1,cc(j):cc(j+1)-1). In terms of a linear
% system, [A11 A12] is the underdetermined part of the system (it is always
% rectangular and with more columns and rows, or 0-by-0), A23 is the well-
% determined part of the system (it is always square), and [A34 ; A44] is
% the over-determined part of the system (it is always rectangular with more
% rows than columns, or 0-by-0).
%
% The structural rank of A is sprank(A) = rr(4)-1, which is an upper bound on
% the numerical rank of A. sprank(A) = rank(full(sprand(A))) with probability
% 1 in exact arithmetic.
%
% The A23 submatrix is further subdivided into block upper triangular form
% via the "fine" decomposition (the strongly-connected components of A23).
% If A is square and structurally non-singular, A23 is the entire matrix.
%
% C(r(i):r(i+1)-1,s(j):s(j+1)-1) is the (i,j)th block of the fine
% decomposition. The (1,1) block is the rectangular block [A11 A12], unless
% this block is 0-by-0. The (b,b) block is the rectangular block [A34 ; A44],
% unless this block is 0-by-0, where b = length(r)-1. All other blocks of the
% form C(r(i):r(i+1)-1,s(i):s(i+1)-1) are diagonal blocks of A23, and are
% square with a zero-free diagonal.
%
% The matching algorithm used in cs_dmperm can take a very long time
% in rare cases. This can be avoided by exploiting a randomized algorithm,
% with cs_dmperm(A,seed). If seed=0, the non-randomized algorithm is used
% (columns are considered in order 1:n). If seed=-1, columns are considered
% in reverse order. Otherwise, the columns are considered in a random order,
% using seed as the random number generator seed. Try cs_dmpmerm(A,1) or
% cs_dmperm(A,rand), for a randomized order, for example. Seed defaults to 0.
%
%
% Example:
% Prob = ssget ('HB/west0479') ; A = Prob.A ; cspy (A) ;
% p = cs_dmperm (A) ;
% cspy (A (p,:)) ;
% [p q r s cc rr] = cs_dmperm (A) ;
% cspy (A (p,q)) ;
% cs_dmspy (A) ;
%
% See also CS_DMSPY, CS_DMSOL, DMPERM, SPRANK, CS_RANDPERM, RAND
% CXSparse, Copyright (c) 2006-2022, Timothy A. Davis. All Rights Reserved.
% SPDX-License-Identifier: LGPL-2.1+
error ('cs_dmperm mexFunction not found') ;
|