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 // NOTE: It is possible to simply include "elemental.hpp" instead
10 #include "elemental-lite.hpp"
11 #include ELEM_HERMITIANTRIDIAGEIG_INC
12 #include ELEM_ZEROS_INC
13 using namespace elem;
14
15 /*
16 This is a testbed for some of the formulas for
17 Prolate Spheroidal Wave Functions (PSFWs) found in Osipov and Rokhlin's
18 "On the evaluation of prolate spheroidal wave functions and associated
19 quadrature rules", arXiv:1301.1707.
20 */
21
22 int
main(int argc,char * argv[])23 main( int argc, char* argv[] )
24 {
25 Initialize( argc, argv );
26
27 try
28 {
29 const Int n = Input("--n","maximum k value",10);
30 const double c = Input("--c","coefficient for PSFW",100.);
31 const bool display = Input("--display","display matrix?",true);
32 const bool print = Input("--print","print matrix?",false);
33 ProcessInput();
34 PrintInputReport();
35
36 const Int nOdd = n/2;
37 const Int nEven = n - nOdd;
38 DistMatrix<double> AEven, AOdd;
39 Zeros( AEven, nEven, nEven );
40 Zeros( AOdd, nOdd, nOdd );
41
42 // Fill AEven and AOdd
43 for( Int k=0; k<n; ++k )
44 {
45 const Int kHalf = k/2;
46 const double kd = k;
47 const double diag =
48 kd*(kd+1) + (2*kd*(kd+1)-1)/((2*kd+3)*(2*kd-1))*c*c;
49 const double offDiag =
50 (kd+2)*(kd+1)/((2*kd+3)*Sqrt((2*kd+1)*(2*kd+5)))*c*c;
51 if( k % 2 == 0 )
52 {
53 // Fill relevant entries of AEven
54 AEven.Set(kHalf,kHalf,diag);
55 if( kHalf+1 < nEven )
56 {
57 AEven.Set(kHalf+1,kHalf, offDiag);
58 AEven.Set(kHalf, kHalf+1,offDiag);
59 }
60 }
61 else
62 {
63 // Fill relevant entries of AOdd
64 AOdd.Set(kHalf,kHalf,diag);
65 if( kHalf+1 < nOdd )
66 {
67 AOdd.Set(kHalf+1,kHalf, offDiag);
68 AOdd.Set(kHalf, kHalf+1,offDiag);
69 }
70 }
71 }
72 if( display )
73 {
74 Display( AEven, "Even matrix" );
75 Display( AOdd, "Odd matrix" );
76 }
77 if( print )
78 {
79 Print( AEven, "AEven" );
80 Print( AOdd, "AOdd" );
81 }
82
83 DistMatrix<double,VR, STAR> wEven, wOdd;
84 DistMatrix<double> XEven, XOdd;
85 HermitianTridiagEig
86 ( AEven.GetDiagonal(), AEven.GetDiagonal(-1), wEven, XEven, ASCENDING );
87 HermitianTridiagEig
88 ( AOdd.GetDiagonal(), AOdd.GetDiagonal(-1), wOdd, XOdd, ASCENDING );
89 if( display )
90 {
91 Display( XEven, "Even eigenvectors" );
92 Display( XOdd, "Odd eigenvectors" );
93 Display( wEven, "Even eigenvalues" );
94 Display( wOdd, "Odd eigenvalues" );
95 }
96 if( print )
97 {
98 Print( XEven, "XEven" );
99 Print( XOdd, "XOdd" );
100 Print( wEven, "wEven" );
101 Print( wOdd, "wOdd" );
102 }
103 }
104 catch( std::exception& e ) { ReportException(e); }
105
106 Finalize();
107 return 0;
108 }
109