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_SYLVESTER_HPP
11 #define ELEM_SYLVESTER_HPP
12
13 #include ELEM_AXPY_INC
14 #include ELEM_SCALE_INC
15 #include ELEM_UPDATEDIAGONAL_INC
16
17 #include ELEM_FROBENIUSNORM_INC
18 #include ELEM_SIGN_INC
19
20 #include ELEM_ZEROS_INC
21
22 namespace elem {
23
24 // W = | A -C |, where A is m x m, B is n x n, and both are assumed to have
25 // | 0 -B | all of their eigenvalues in the open right-half plane.
26 //
27 // The solution, X, to the equation
28 // A X + X B = C
29 // is returned, as well as the number of Newton iterations for computing sgn(W).
30 //
31 // See Chapter 2 of Nicholas J. Higham's "Functions of Matrices"
32
33 template<typename F>
34 inline void
Sylvester(Int m,Matrix<F> & W,Matrix<F> & X,SignCtrl<Base<F>> signCtrl=SignCtrl<Base<F>> ())35 Sylvester
36 ( Int m, Matrix<F>& W, Matrix<F>& X,
37 SignCtrl<Base<F>> signCtrl=SignCtrl<Base<F>>() )
38 {
39 DEBUG_ONLY(CallStackEntry cse("Sylvester"))
40 Sign( W, signCtrl );
41 Matrix<F> WTL, WTR,
42 WBL, WBR;
43 PartitionDownDiagonal
44 ( W, WTL, WTR,
45 WBL, WBR, m );
46 // WTL and WBR should be the positive and negative identity, WBL should be
47 // zero, and WTR should be -2 X
48 X = WTR;
49 Scale( -F(1)/F(2), X );
50
51 // TODO: Think of how to probe for checks on other quadrants.
52 /*
53 typedef Base<F> Real;
54 UpdateDiagonal( WTL, F(-1) );
55 const Real errorWTL = FrobeniusNorm( WTL );
56 const Int n = W.Height() - m;
57 UpdateDiagonal( WBR, F(1) );
58 const Real errorWBR = FrobeniusNorm( WBR );
59 const Real errorWBL = FrobeniusNorm( WBL );
60 */
61 }
62
63 template<typename F>
64 inline void
Sylvester(Int m,DistMatrix<F> & W,DistMatrix<F> & X,SignCtrl<Base<F>> signCtrl=SignCtrl<Base<F>> ())65 Sylvester
66 ( Int m, DistMatrix<F>& W, DistMatrix<F>& X,
67 SignCtrl<Base<F>> signCtrl=SignCtrl<Base<F>>() )
68 {
69 DEBUG_ONLY(CallStackEntry cse("Sylvester"))
70 const Grid& g = W.Grid();
71 Sign( W, signCtrl );
72 DistMatrix<F> WTL(g), WTR(g),
73 WBL(g), WBR(g);
74 PartitionDownDiagonal
75 ( W, WTL, WTR,
76 WBL, WBR, m );
77 // WTL and WBR should be the positive and negative identity, WBL should be
78 // zero, and WTR should be -2 X
79 X = WTR;
80 Scale( -F(1)/F(2), X );
81
82 // TODO: Think of how to probe for checks on other quadrants.
83 // Add UpdateDiagonal routine to avoid explicit identity Axpy?
84 /*
85 typedef Base<F> Real;
86 UpdateDiagonal( WTL, F(-1) );
87 const Real errorWTL = FrobeniusNorm( WTL );
88 const Int n = W.Height() - m;
89 UpdateDiagonal( WBR, F(1) );
90 const Real errorWBR = FrobeniusNorm( WBR );
91 const Real errorWBL = FrobeniusNorm( WBL );
92 */
93 }
94
95 template<typename F>
96 inline void
Sylvester(const Matrix<F> & A,const Matrix<F> & B,const Matrix<F> & C,Matrix<F> & X,SignCtrl<Base<F>> signCtrl=SignCtrl<Base<F>> ())97 Sylvester
98 ( const Matrix<F>& A, const Matrix<F>& B, const Matrix<F>& C, Matrix<F>& X,
99 SignCtrl<Base<F>> signCtrl=SignCtrl<Base<F>>() )
100 {
101 DEBUG_ONLY(
102 CallStackEntry cse("Sylvester");
103 if( A.Height() != A.Width() )
104 LogicError("A must be square");
105 if( B.Height() != B.Width() )
106 LogicError("B must be square");
107 if( C.Height() != A.Height() || C.Width() != B.Height() )
108 LogicError("C must conform with A and B");
109 )
110 const Int m = C.Height();
111 const Int n = C.Width();
112 Matrix<F> W, WTL, WTR,
113 WBL, WBR;
114 Zeros( W, m+n, m+n );
115 PartitionDownDiagonal
116 ( W, WTL, WTR,
117 WBL, WBR, m );
118 WTL = A;
119 WBR = B; Scale( F(-1), WBR );
120 WTR = C; Scale( F(-1), WTR );
121 Sylvester( m, W, X, signCtrl );
122 }
123
124 template<typename F>
125 inline void
Sylvester(const DistMatrix<F> & A,const DistMatrix<F> & B,const DistMatrix<F> & C,DistMatrix<F> & X,SignCtrl<Base<F>> signCtrl=SignCtrl<Base<F>> ())126 Sylvester
127 ( const DistMatrix<F>& A, const DistMatrix<F>& B, const DistMatrix<F>& C,
128 DistMatrix<F>& X, SignCtrl<Base<F>> signCtrl=SignCtrl<Base<F>>() )
129 {
130 DEBUG_ONLY(
131 CallStackEntry cse("Sylvester");
132 if( A.Height() != A.Width() )
133 LogicError("A must be square");
134 if( B.Height() != B.Width() )
135 LogicError("B must be square");
136 if( C.Height() != A.Height() || C.Width() != B.Height() )
137 LogicError("C must conform with A and B");
138 if( A.Grid() != B.Grid() || B.Grid() != C.Grid() )
139 LogicError("A, B, and C must have the same grid");
140 )
141 const Int m = C.Height();
142 const Int n = C.Width();
143 const Grid& g = A.Grid();
144 DistMatrix<F> W(g), WTL(g), WTR(g),
145 WBL(g), WBR(g);
146 Zeros( W, m+n, m+n );
147 PartitionDownDiagonal
148 ( W, WTL, WTR,
149 WBL, WBR, m );
150 WTL = A;
151 WBR = B; Scale( F(-1), WBR );
152 WTR = C; Scale( F(-1), WTR );
153 Sylvester( m, W, X, signCtrl );
154 }
155
156 } // namespace elem
157
158 #endif // ifndef ELEM_SYLVESTER_HPP
159