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_FOXLI_HPP
11 #define ELEM_FOXLI_HPP
12
13 #include ELEM_DIAGONALSCALE_INC
14 #include ELEM_HERMITIANTRIDIAGEIG_INC
15 #include ELEM_ZEROS_INC
16
17 namespace elem {
18
19 template<typename Real>
20 inline void
FoxLi(Matrix<Complex<Real>> & A,Int n,Real omega)21 FoxLi( Matrix<Complex<Real>>& A, Int n, Real omega )
22 {
23 DEBUG_ONLY(CallStackEntry cse("FoxLi"))
24 typedef Complex<Real> C;
25 const Real pi = 4*Atan( Real(1) );
26 const C phi = Sqrt( C(0,omega/pi) );
27
28 // Compute Gauss quadrature points and weights
29 Matrix<Real> d, e;
30 Zeros( d, n, 1 );
31 e.Resize( n-1, 1 );
32 for( Int j=0; j<n-1; ++j )
33 {
34 const Real betaInv = 2*Sqrt(1-Pow(j+Real(1),-2)/4);
35 e.Set( j, 0, 1/betaInv );
36 }
37 Matrix<Real> x, Z;
38 HermitianTridiagEig( d, e, x, Z, UNSORTED );
39 auto z = LockedView( Z, 0, 0, 1, n );
40 Matrix<Real> sqrtWeights( z );
41 for( Int j=0; j<n; ++j )
42 sqrtWeights.Set( 0, j, Sqrt(2)*Abs(sqrtWeights.Get(0,j)) );
43 herm_eig::Sort( x, sqrtWeights, ASCENDING );
44 Transpose( sqrtWeights );
45
46 // Form the integral operator
47 A.Resize( n, n );
48 for( Int j=0; j<n; ++j )
49 {
50 for( Int i=0; i<n; ++i )
51 {
52 const Real theta = -omega*Pow(x.Get(i,0)-x.Get(j,0),2);
53 const Real realPart = Cos(theta);
54 const Real imagPart = Sin(theta);
55 A.Set( i, j, phi*C(realPart,imagPart) );
56 }
57 }
58
59 // Apply the weighting
60 DiagonalScale( LEFT, NORMAL, sqrtWeights, A );
61 DiagonalScale( RIGHT, NORMAL, sqrtWeights, A );
62 }
63
64 template<typename Real,Dist U,Dist V>
65 inline void
FoxLi(DistMatrix<Complex<Real>,U,V> & A,Int n,Real omega)66 FoxLi( DistMatrix<Complex<Real>,U,V>& A, Int n, Real omega )
67 {
68 DEBUG_ONLY(CallStackEntry cse("FoxLi"))
69 typedef Complex<Real> C;
70 const Real pi = 4*Atan( Real(1) );
71 const C phi = Sqrt( C(0,omega/pi) );
72
73 // Compute Gauss quadrature points and weights
74 const Grid& g = A.Grid();
75 DistMatrix<Real,VR,STAR> d(g), e(g);
76 Zeros( d, n, 1 );
77 e.Resize( n-1, 1 );
78 for( Int iLoc=0; iLoc<e.LocalHeight(); ++iLoc )
79 {
80 const Int i = e.GlobalRow(iLoc);
81 const Real betaInv = 2*Sqrt(1-Pow(i+Real(1),-2)/4);
82 e.SetLocal( iLoc, 0, 1/betaInv );
83 }
84 DistMatrix<Real,VR,STAR> x(g);
85 DistMatrix<Real,STAR,VR> Z(g);
86 HermitianTridiagEig( d, e, x, Z, UNSORTED );
87 auto z = LockedView( Z, 0, 0, 1, n );
88 DistMatrix<Real,STAR,VR> sqrtWeights( z );
89 for( Int jLoc=0; jLoc<sqrtWeights.LocalWidth(); ++jLoc )
90 sqrtWeights.SetLocal
91 ( 0, jLoc, Sqrt(2)*Abs(sqrtWeights.GetLocal(0,jLoc)) );
92 herm_eig::Sort( x, sqrtWeights, ASCENDING );
93
94 // Form the integral operator
95 A.Resize( n, n );
96 DistMatrix<Real,U,STAR> x_U_STAR( x );
97 DistMatrix<Real,V,STAR> x_V_STAR( x );
98 for( Int jLoc=0; jLoc<A.LocalWidth(); ++jLoc )
99 {
100 for( Int iLoc=0; iLoc<A.LocalHeight(); ++iLoc )
101 {
102 const Real diff = x_U_STAR.GetLocal(iLoc,0)-
103 x_V_STAR.GetLocal(jLoc,0);
104 const Real theta = -omega*Pow(diff,2);
105 const Real realPart = Cos(theta);
106 const Real imagPart = Sin(theta);
107 A.SetLocal( iLoc, jLoc, phi*C(realPart,imagPart) );
108 }
109 }
110
111 // Apply the weighting
112 DistMatrix<Real,VR,STAR> sqrtWeightsTrans(g);
113 Transpose( sqrtWeights, sqrtWeightsTrans );
114 DiagonalScale( LEFT, NORMAL, sqrtWeightsTrans, A );
115 DiagonalScale( RIGHT, NORMAL, sqrtWeightsTrans, A );
116 }
117
118 } // namespace elem
119
120 #endif // ifndef ELEM_FOXLI_HPP
121