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_HERMITIANUNIFORMSPECTRUM_HPP
11 #define ELEM_HERMITIANUNIFORMSPECTRUM_HPP
12
13 #include "./Diagonal.hpp"
14 #include "./Haar.hpp"
15 #include "./Uniform.hpp"
16
17 namespace elem {
18
19 // Draw the spectrum from the specified half-open interval on the real line,
20 // then rotate with a Haar matrix
21
22 template<typename F>
23 inline void
MakeHermitianUniformSpectrum(Matrix<F> & A,Base<F> lower=0,Base<F> upper=1)24 MakeHermitianUniformSpectrum( Matrix<F>& A, Base<F> lower=0, Base<F> upper=1 )
25 {
26 DEBUG_ONLY(CallStackEntry cse("MakeHermitianUniformSpectrum"))
27 if( A.Height() != A.Width() )
28 LogicError("Cannot make a non-square matrix Hermitian");
29 typedef Base<F> Real;
30 const bool isComplex = IsComplex<F>::val;
31
32 // Form d and D
33 const Int n = A.Height();
34 std::vector<F> d( n );
35 for( Int j=0; j<n; ++j )
36 d[j] = SampleUniform<Real>( lower, upper );
37 Diagonal( A, d );
38
39 // Apply a Haar matrix from both sides
40 Matrix<F> Q, t;
41 Matrix<Real> s;
42 ImplicitHaar( Q, t, s, n );
43 qr::ApplyQ( LEFT, NORMAL, Q, t, s, A );
44 qr::ApplyQ( RIGHT, ADJOINT, Q, t, s, A );
45
46 if( isComplex )
47 {
48 const Int height = A.Height();
49 for( Int j=0; j<height; ++j )
50 A.SetImagPart( j, j, Real(0) );
51 }
52 }
53
54 template<typename F,Dist U,Dist V>
55 inline void
MakeHermitianUniformSpectrum(DistMatrix<F,U,V> & A,Base<F> lower=0,Base<F> upper=1)56 MakeHermitianUniformSpectrum
57 ( DistMatrix<F,U,V>& A, Base<F> lower=0, Base<F> upper=1 )
58 {
59 DEBUG_ONLY(CallStackEntry cse("MakeHermitianUniformSpectrum"))
60 if( A.Height() != A.Width() )
61 LogicError("Cannot make a non-square matrix Hermitian");
62 const Grid& grid = A.Grid();
63 typedef Base<F> Real;
64 const bool isComplex = IsComplex<F>::val;
65 const bool standardDist = ( U == MC && V == MR );
66
67 // Form d and D
68 const Int n = A.Height();
69 std::vector<F> d( n );
70 if( grid.Rank() == 0 )
71 for( Int j=0; j<n; ++j )
72 d[j] = SampleUniform<Real>( lower, upper );
73 mpi::Broadcast( d.data(), n, 0, grid.Comm() );
74 DistMatrix<F> ABackup( grid );
75 if( standardDist )
76 Diagonal( A, d );
77 else
78 {
79 ABackup.AlignWith( A );
80 Diagonal( ABackup, d );
81 }
82
83 // Apply a Haar matrix from both sides
84 DistMatrix<F> Q(grid);
85 DistMatrix<F,MD,STAR> t(grid);
86 DistMatrix<Real,MD,STAR> s(grid);
87 ImplicitHaar( Q, t, s, n );
88
89 // Copy the result into the correct distribution
90 if( standardDist )
91 {
92 qr::ApplyQ( LEFT, NORMAL, Q, t, s, A );
93 qr::ApplyQ( RIGHT, ADJOINT, Q, t, s, A );
94 }
95 else
96 {
97 qr::ApplyQ( LEFT, NORMAL, Q, t, s, ABackup );
98 qr::ApplyQ( RIGHT, ADJOINT, Q, t, s, ABackup );
99 A = ABackup;
100 }
101
102 // Force the diagonal to be real-valued
103 if( isComplex )
104 {
105 const Int localHeight = A.LocalHeight();
106 const Int localWidth = A.LocalWidth();
107 for( Int jLoc=0; jLoc<localWidth; ++jLoc )
108 {
109 const Int j = A.GlobalCol(jLoc);
110 for( Int iLoc=0; iLoc<localHeight; ++iLoc )
111 {
112 const Int i = A.GlobalRow(iLoc);
113 if( i == j )
114 A.SetLocalImagPart( iLoc, jLoc, Real(0) );
115 }
116 }
117 }
118 }
119
120 template<typename F>
121 inline void
HermitianUniformSpectrum(Matrix<F> & A,Int n,Base<F> lower=0,Base<F> upper=1)122 HermitianUniformSpectrum
123 ( Matrix<F>& A, Int n, Base<F> lower=0, Base<F> upper=1 )
124 {
125 DEBUG_ONLY(CallStackEntry cse("HermitianUniformSpectrum"))
126 A.Resize( n, n );
127 MakeHermitianUniformSpectrum( A, lower, upper );
128 }
129
130 template<typename F,Dist U,Dist V>
131 inline void
HermitianUniformSpectrum(DistMatrix<F,U,V> & A,Int n,Base<F> lower=0,Base<F> upper=1)132 HermitianUniformSpectrum
133 ( DistMatrix<F,U,V>& A, Int n, Base<F> lower=0, Base<F> upper=1 )
134 {
135 DEBUG_ONLY(CallStackEntry cse("HermitianUniformSpectrum"))
136 A.Resize( n, n );
137 MakeHermitianUniformSpectrum( A, lower, upper );
138 }
139
140 } // namespace elem
141
142 #endif // ifndef ELEM_HERMITIANUNIFORMSPECTRUM_HPP
143