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