1 /*
2 Copyright (c) 2009-2014, Jack Poulson
3 All rights reserved.
4
5 This file is part of Elemental and is under the BSD 2-Clause License,
6 which can be found in the LICENSE file in the root directory, or at
7 http://opensource.org/licenses/BSD-2-Clause
8 */
9 #pragma once
10 #ifndef ELEM_PIVOTPARITY_HPP
11 #define ELEM_PIVOTPARITY_HPP
12
13 namespace elem {
14
15 // A permutation is even if and only if it is the product of an even number
16 // of transpositions, so we can decide this by simply checking how many
17 // nontrivial pivots were performed.
18
19 inline bool
PivotParity(const Matrix<Int> & p,Int pivotOffset=0)20 PivotParity( const Matrix<Int>& p, Int pivotOffset=0 )
21 {
22 DEBUG_ONLY(
23 CallStackEntry cse("PivotParity");
24 if( p.Width() != 1 )
25 LogicError("p must be a column vector");
26 if( pivotOffset < 0 )
27 LogicError("pivot offset cannot be negative");
28 )
29 const Int n = p.Height();
30 bool isOdd = false;
31 for( Int i=0; i<n; ++i )
32 if( p.Get(i,0) != i+pivotOffset )
33 isOdd = !isOdd;
34 return isOdd;
35 }
36
37 inline bool
PivotParity(const DistMatrix<Int,VC,STAR> & p,Int pivotOffset=0)38 PivotParity( const DistMatrix<Int,VC,STAR>& p, Int pivotOffset=0 )
39 {
40 DEBUG_ONLY(
41 CallStackEntry cse("PivotParity");
42 if( p.Width() != 1 )
43 LogicError("p must be a column vector");
44 if( pivotOffset < 0 )
45 LogicError("pivot offset cannot be negative");
46 )
47 bool isLocallyOdd = false;
48 const Int mLocal = p.LocalHeight();
49 for( Int iLoc=0; iLoc<mLocal; ++iLoc )
50 {
51 const Int i = p.GlobalRow(iLoc);
52 if( p.GetLocal(iLoc,0) != i+pivotOffset )
53 isLocallyOdd = !isLocallyOdd;
54 }
55
56 Int localContribution = isLocallyOdd;
57 Int globalContribution;
58 mpi::AllReduce
59 ( &localContribution, &globalContribution, 1, MPI_SUM, p.ColComm() );
60 return globalContribution % 2;
61 }
62
63 } // namespace elem
64
65 #endif // ifndef ELEM_PIVOTPARITY_HPP
66