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_UNIFORMHELMHOLTZGREENS_HPP
11 #define ELEM_UNIFORMHELMHOLTZGREENS_HPP
12
13 namespace elem {
14
15 // Generate the "random Green's matrix" from
16 // A. Goetschy and S. E. Skipetrov's "Non-Hermitian Euclidean random matrix
17 // theory".
18 // It is essentially the 3D Helmholtz Green's function with source and target
19 // points chosen uniformly from the unit ball. The behaviour of the spectrum
20 // should change dramatically dependinging upon the number of points per
21 // wavelength-cubed.
22
23 template<typename Real>
24 inline void
UniformHelmholtzGreens(Matrix<Complex<Real>> & A,Int n,Real lambda)25 UniformHelmholtzGreens( Matrix<Complex<Real>>& A, Int n, Real lambda )
26 {
27 DEBUG_ONLY(CallStackEntry cse("UniformHelmholtzGreens"))
28 typedef Complex<Real> C;
29 const Real pi = 4*Atan( Real(1) );
30 const Real k0 = 2*pi/lambda;
31
32 // Generate a list of n uniform samples from the 3D unit ball
33 Matrix<Real> X(3,n);
34 for( Int j=0; j<n; ++j )
35 {
36 Real x0, x1, x2;
37 // Sample uniformly from [-1,+1]^3 until a point is drawn from the ball
38 while( true )
39 {
40 x0 = SampleUniform( Real(-1), Real(1) );
41 x1 = SampleUniform( Real(-1), Real(1) );
42 x2 = SampleUniform( Real(-1), Real(1) );
43 const Real radiusSq = x0*x0 + x1*x1 + x2*x2;
44 if( radiusSq > 0 && radiusSq <= 1 )
45 break;
46 }
47 X.Set( 0, j, x0 );
48 X.Set( 1, j, x1 );
49 X.Set( 2, j, x2 );
50 }
51
52 A.Resize( n, n );
53 for( Int j=0; j<n; ++j )
54 {
55 const Real xj0 = X.Get(0,j);
56 const Real xj1 = X.Get(1,j);
57 const Real xj2 = X.Get(2,j);
58 for( Int i=0; i<n; ++i )
59 {
60 if( i == j )
61 {
62 A.Set( i, j, 0 );
63 }
64 else
65 {
66 const Real d0 = X.Get(0,i)-xj0;
67 const Real d1 = X.Get(1,i)-xj1;
68 const Real d2 = X.Get(2,i)-xj2;
69 const Real gamma = k0*Sqrt(d0*d0+d1*d1+d2*d2);
70 const Real realPart = Cos(gamma)/gamma;
71 const Real imagPart = Sin(gamma)/gamma;
72 A.Set( i, j, C(realPart,imagPart) );
73 }
74 }
75 }
76 }
77
78 template<typename Real,Dist U,Dist V>
79 inline void
UniformHelmholtzGreens(DistMatrix<Complex<Real>,U,V> & A,Int n,Real lambda)80 UniformHelmholtzGreens( DistMatrix<Complex<Real>,U,V>& A, Int n, Real lambda )
81 {
82 DEBUG_ONLY(CallStackEntry cse("UniformHelmholtzGreens"))
83 typedef Complex<Real> C;
84 const Real pi = 4*Atan( Real(1) );
85 const Real k0 = 2*pi/lambda;
86 const Grid& g = A.Grid();
87
88 // Generate a list of n uniform samples from the 3D unit ball
89 DistMatrix<Real,STAR,VR> X_STAR_VR(3,n,g);
90 for( Int jLoc=0; jLoc<X_STAR_VR.LocalWidth(); ++jLoc )
91 {
92 Real x0, x1, x2;
93 // Sample uniformly from [-1,+1]^3 until a point is drawn from the ball
94 while( true )
95 {
96 x0 = SampleUniform( Real(-1), Real(1) );
97 x1 = SampleUniform( Real(-1), Real(1) );
98 x2 = SampleUniform( Real(-1), Real(1) );
99 const Real radiusSq = x0*x0 + x1*x1 + x2*x2;
100 if( radiusSq > 0 && radiusSq <= 1 )
101 break;
102 }
103 X_STAR_VR.SetLocal( 0, jLoc, x0 );
104 X_STAR_VR.SetLocal( 1, jLoc, x1 );
105 X_STAR_VR.SetLocal( 2, jLoc, x2 );
106 }
107 DistMatrix<Real,STAR,STAR> X_STAR_STAR( X_STAR_VR );
108
109 A.Resize( n, n );
110 for( Int jLoc=0; jLoc<A.LocalWidth(); ++jLoc )
111 {
112 const Int j = A.GlobalCol(jLoc);
113 const Real xj0 = X_STAR_STAR.GetLocal(0,j);
114 const Real xj1 = X_STAR_STAR.GetLocal(1,j);
115 const Real xj2 = X_STAR_STAR.GetLocal(2,j);
116 for( Int iLoc=0; iLoc<A.LocalHeight(); ++iLoc )
117 {
118 const Int i = A.GlobalRow(iLoc);
119 if( i == j )
120 {
121 A.SetLocal( iLoc, jLoc, 0 );
122 }
123 else
124 {
125 const Real d0 = X_STAR_STAR.GetLocal(0,i)-xj0;
126 const Real d1 = X_STAR_STAR.GetLocal(1,i)-xj1;
127 const Real d2 = X_STAR_STAR.GetLocal(2,i)-xj2;
128 const Real gamma = k0*Sqrt(d0*d0+d1*d1+d2*d2);
129 const Real realPart = Cos(gamma)/gamma;
130 const Real imagPart = Sin(gamma)/gamma;
131 A.SetLocal( iLoc, jLoc, C(realPart,imagPart) );
132 }
133 }
134 }
135 }
136
137 template<typename Real,Dist U,Dist V>
138 inline void
UniformHelmholtzGreens(BlockDistMatrix<Complex<Real>,U,V> & A,Int n,Real lambda)139 UniformHelmholtzGreens
140 ( BlockDistMatrix<Complex<Real>,U,V>& A, Int n, Real lambda )
141 {
142 DEBUG_ONLY(CallStackEntry cse("UniformHelmholtzGreens"))
143 typedef Complex<Real> C;
144 const Real pi = 4*Atan( Real(1) );
145 const Real k0 = 2*pi/lambda;
146 const Grid& g = A.Grid();
147
148 // Generate a list of n uniform samples from the 3D unit ball
149 DistMatrix<Real,STAR,VR> X_STAR_VR(3,n,g);
150 for( Int jLoc=0; jLoc<X_STAR_VR.LocalWidth(); ++jLoc )
151 {
152 Real x0, x1, x2;
153 // Sample uniformly from [-1,+1]^3 until a point is drawn from the ball
154 while( true )
155 {
156 x0 = SampleUniform( Real(-1), Real(1) );
157 x1 = SampleUniform( Real(-1), Real(1) );
158 x2 = SampleUniform( Real(-1), Real(1) );
159 const Real radiusSq = x0*x0 + x1*x1 + x2*x2;
160 if( radiusSq > 0 && radiusSq <= 1 )
161 break;
162 }
163 X_STAR_VR.SetLocal( 0, jLoc, x0 );
164 X_STAR_VR.SetLocal( 1, jLoc, x1 );
165 X_STAR_VR.SetLocal( 2, jLoc, x2 );
166 }
167 DistMatrix<Real,STAR,STAR> X_STAR_STAR( X_STAR_VR );
168
169 A.Resize( n, n );
170 for( Int jLoc=0; jLoc<A.LocalWidth(); ++jLoc )
171 {
172 const Int j = A.GlobalCol(jLoc);
173 const Real xj0 = X_STAR_STAR.GetLocal(0,j);
174 const Real xj1 = X_STAR_STAR.GetLocal(1,j);
175 const Real xj2 = X_STAR_STAR.GetLocal(2,j);
176 for( Int iLoc=0; iLoc<A.LocalHeight(); ++iLoc )
177 {
178 const Int i = A.GlobalRow(iLoc);
179 if( i == j )
180 {
181 A.SetLocal( iLoc, jLoc, 0 );
182 }
183 else
184 {
185 const Real d0 = X_STAR_STAR.GetLocal(0,i)-xj0;
186 const Real d1 = X_STAR_STAR.GetLocal(1,i)-xj1;
187 const Real d2 = X_STAR_STAR.GetLocal(2,i)-xj2;
188 const Real gamma = k0*Sqrt(d0*d0+d1*d1+d2*d2);
189 const Real realPart = Cos(gamma)/gamma;
190 const Real imagPart = Sin(gamma)/gamma;
191 A.SetLocal( iLoc, jLoc, C(realPart,imagPart) );
192 }
193 }
194 }
195 }
196
197 } // namespace elem
198
199 #endif // ifndef ELEM_UNIFORMHELMHOLTZGREENS_HPP
200