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_HYPERBOLICREFLECTOR_HPP
11 #define ELEM_HYPERBOLICREFLECTOR_HPP
12 
13 #include ELEM_CONJUGATE_INC
14 #include ELEM_NRM2_INC
15 #include ELEM_SCALE_INC
16 #include ELEM_ZERO_INC
17 
18 #include "./HyperbolicReflector/Col.hpp"
19 #include "./HyperbolicReflector/Row.hpp"
20 
21 namespace elem {
22 
23 //
24 // (I - 1/tau w w^H Sigma) x = -lambda e_0
25 //
26 // where Sigma is diag(+1,-1,...,-1)
27 
28 template<typename F>
29 inline F
LeftHyperbolicReflector(Matrix<F> & chi,Matrix<F> & x)30 LeftHyperbolicReflector( Matrix<F>& chi, Matrix<F>& x )
31 {
32     DEBUG_ONLY(
33         CallStackEntry cse("LeftHyperbolicReflector");
34         if( chi.Height() != 1 || chi.Width() != 1 )
35             LogicError("chi must be a scalar");
36         if( x.Height() != 1 && x.Width() != 1 )
37             LogicError("x must be a vector");
38         if( ImagPart(chi.Get(0,0)) != Base<F>(0) )
39             LogicError("chi is assumed to be real");
40     )
41 
42     // Compute lambda = sgn(chi) sqrt([chi;x]^H Sigma [chi;x])
43     //                = sgn(chi) sqrt(chi^2 - x^H x)
44     typedef Base<F> Real;
45     const Real alpha = chi.GetRealPart(0,0);
46     const Real xNrm = Nrm2( x );
47     const Real delta = alpha*alpha - xNrm*xNrm;
48     if( delta < Real(0) )
49         LogicError("Attempted to square-root a negative number");
50     const Real lambda = ( alpha>=0 ? Sqrt(delta) : -Sqrt(delta) );
51     chi.Set(0,0,-lambda);
52 
53     // Implicitly define
54     //     w := [chi;x] + lambda e_0, and
55     //     kappa = chi + lambda,
56     // so that
57     //     tau := (w^H Sigma w) / 2 = (delta + lambda^2 + 2 chi lambda) / 2
58     //          = delta + chi lambda
59     // then normalize w so that its first entry is one
60     // TODO: Introduce a threshold instead of the approach from
61     //       van de Geijn and van Zee's "High-performance up-and-downdating
62     //       via Householder-like transformations"
63     const Real kappa = alpha + lambda;
64     if( kappa == Real(0) )
65     {
66         Zero( x );
67         return Real(1);
68     }
69     else
70     {
71         Scale( Real(1)/kappa, x );
72         return (delta+alpha*lambda)/(kappa*kappa);
73     }
74 }
75 
76 template<typename F>
77 inline F
LeftHyperbolicReflector(F & chi,Matrix<F> & x)78 LeftHyperbolicReflector( F& chi, Matrix<F>& x )
79 {
80     DEBUG_ONLY(
81         CallStackEntry cse("LeftHyperbolicReflector");
82         if( x.Height() != 1 && x.Width() != 1 )
83             LogicError("x must be a vector");
84         if( ImagPart(chi) != Base<F>(0) )
85             LogicError("chi is assumed to be real");
86     )
87 
88     // Compute lambda = sgn(chi) sqrt([chi;x]^H Sigma [chi;x])
89     //                = sgn(chi) sqrt(chi^2 - x^H x)
90     typedef Base<F> Real;
91     const Real alpha = RealPart(chi);
92     const Real xNrm = Nrm2( x );
93     const Real delta = alpha*alpha - xNrm*xNrm;
94     if( delta < Real(0) )
95         LogicError("Attempted to square-root a negative number");
96     const Real lambda = ( alpha>=0 ? Sqrt(delta) : -Sqrt(delta) );
97     chi = -lambda;
98 
99     // Implicitly define
100     //     w := [chi;x] + lambda e_0, and
101     //     kappa = chi + lambda,
102     // so that
103     //     tau := (w^H Sigma w) / 2 = (delta + lambda^2 + 2 chi lambda) / 2
104     //          = delta + chi lambda
105     // then normalize w so that its first entry is one
106     // TODO: Introduce a threshold instead of the approach from
107     //       van de Geijn and van Zee's "High-performance up-and-downdating
108     //       via Householder-like transformations"
109     const Real kappa = alpha + lambda;
110     if( kappa == Real(0) )
111     {
112         Zero( x );
113         return Real(1);
114     }
115     else
116     {
117         Scale( Real(1)/kappa, x );
118         return (delta+alpha*lambda)/(kappa*kappa);
119     }
120 }
121 
122 template<typename F,Dist U,Dist V>
123 inline F
LeftHyperbolicReflector(DistMatrix<F,U,V> & chi,DistMatrix<F,U,V> & x)124 LeftHyperbolicReflector( DistMatrix<F,U,V>& chi, DistMatrix<F,U,V>& x )
125 {
126     DEBUG_ONLY(
127         CallStackEntry cse("LeftHyperbolicReflector");
128         if( chi.Grid() != x.Grid() )
129             LogicError("chi and x must be distributed over the same grid");
130         if( chi.Height() != 1 || chi.Width() != 1 )
131             LogicError("chi must be a scalar");
132         if( x.Width() != 1 )
133             LogicError("x must be a column vector");
134     )
135     F tau;
136     if( x.RowRank() == x.RowAlign() )
137         tau = hyp_reflector::Col( chi, x );
138     mpi::Broadcast( tau, x.RowAlign(), x.RowComm() );
139     return tau;
140 }
141 
142 template<typename F,Dist U,Dist V>
143 inline F
LeftHyperbolicReflector(F & chi,DistMatrix<F,U,V> & x)144 LeftHyperbolicReflector( F& chi, DistMatrix<F,U,V>& x )
145 {
146     DEBUG_ONLY(
147         CallStackEntry cse("LeftHyperbolicReflector");
148         if( x.Width() != 1 )
149             LogicError("x must be a column vector");
150     )
151     F tau;
152     if( x.RowRank() == x.RowAlign() )
153         tau = hyp_reflector::Col( chi, x );
154     mpi::Broadcast( tau, x.RowAlign(), x.RowComm() );
155     return tau;
156 }
157 
158 //
159 // x (I - 1/tau Sigma w^H w) = -lambda e_0^T
160 //
161 // where Sigma is diag(+1,-1,...,-1)
162 
163 template<typename F>
164 inline F
RightHyperbolicReflector(Matrix<F> & chi,Matrix<F> & x)165 RightHyperbolicReflector( Matrix<F>& chi, Matrix<F>& x )
166 {
167     DEBUG_ONLY(CallStackEntry cse("RightHyperbolicReflector"))
168     const F tau = LeftHyperbolicReflector( chi, x );
169     Conjugate( x );
170     return tau;
171 }
172 
173 template<typename F>
174 inline F
RightHyperbolicReflector(F & chi,Matrix<F> & x)175 RightHyperbolicReflector( F& chi, Matrix<F>& x )
176 {
177     DEBUG_ONLY(CallStackEntry cse("RightHyperbolicReflector"))
178     const F tau = LeftHyperbolicReflector( chi, x );
179     Conjugate( x );
180     return tau;
181 }
182 
183 template<typename F,Dist U,Dist V>
184 inline F
RightHyperbolicReflector(DistMatrix<F,U,V> & chi,DistMatrix<F,U,V> & x)185 RightHyperbolicReflector( DistMatrix<F,U,V>& chi, DistMatrix<F,U,V>& x )
186 {
187     DEBUG_ONLY(
188         CallStackEntry cse("RightHyperbolicReflector");
189         if( chi.Grid() != x.Grid() )
190             LogicError("chi and x must be distributed over the same grid");
191         if( chi.Height() != 1 || chi.Width() != 1 )
192             LogicError("chi must be a scalar");
193         if( x.Height() != 1 )
194             LogicError("x must be a row vector");
195     )
196     F tau;
197     if( x.ColRank() == x.ColAlign() )
198         tau = hyp_reflector::Row( chi, x );
199     mpi::Broadcast( tau, x.ColAlign(), x.ColComm() );
200     return tau;
201 }
202 
203 template<typename F,Dist U,Dist V>
204 inline F
RightHyperbolicReflector(F & chi,DistMatrix<F,U,V> & x)205 RightHyperbolicReflector( F& chi, DistMatrix<F,U,V>& x )
206 {
207     DEBUG_ONLY(
208         CallStackEntry cse("RightHyperbolicReflector");
209         if( x.Height() != 1 )
210             LogicError("x must be a row vector");
211     )
212     F tau;
213     if( x.ColRank() == x.ColAlign() )
214         tau = hyp_reflector::Row( chi, x );
215     mpi::Broadcast( tau, x.ColAlign(), x.ColComm() );
216     return tau;
217 }
218 
219 } // namespace elem
220 
221 #endif // ifndef ELEM_HYPERBOLICREFLECTOR_HPP
222