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_MAKETRAPEZOIDAL_INC
12 #include ELEM_MAKETRIANGULAR_INC
13 #include ELEM_BIDIAG_INC
14 #include ELEM_INFINITYNORM_INC
15 #include ELEM_FROBENIUSNORM_INC
16 #include ELEM_IDENTITY_INC
17 #include ELEM_UNIFORM_INC
18 using namespace std;
19 using namespace elem;
20
21 template<typename F>
TestCorrectness(const DistMatrix<F> & A,const DistMatrix<F,STAR,STAR> & tP,const DistMatrix<F,STAR,STAR> & tQ,DistMatrix<F> & AOrig,bool print,bool display)22 void TestCorrectness
23 ( const DistMatrix<F>& A,
24 const DistMatrix<F,STAR,STAR>& tP,
25 const DistMatrix<F,STAR,STAR>& tQ,
26 DistMatrix<F>& AOrig,
27 bool print, bool display )
28 {
29 typedef Base<F> Real;
30 const Grid& g = A.Grid();
31 const Int m = AOrig.Height();
32 const Int n = AOrig.Width();
33 const Real infNormAOrig = InfinityNorm( AOrig );
34 const Real frobNormAOrig = FrobeniusNorm( AOrig );
35 if( g.Rank() == 0 )
36 cout << "Testing error..." << endl;
37
38 // Grab the diagonal and superdiagonal of the bidiagonal matrix
39 auto d = A.GetDiagonal( 0 );
40 auto e = A.GetDiagonal( (m>=n ? 1 : -1) );
41
42 // Zero B and then fill its bidiagonal
43 DistMatrix<F> B(g);
44 B.AlignWith( A );
45 Zeros( B, m, n );
46 B.SetDiagonal( d, 0 );
47 B.SetDiagonal( e, (m>=n ? 1 : -1) );
48 if( print )
49 Print( B, "Bidiagonal" );
50 if( display )
51 Display( B, "Bidiagonal" );
52
53 if( print || display )
54 {
55 DistMatrix<F> Q(g), P(g);
56 Identity( Q, m, m );
57 Identity( P, n, n );
58 bidiag::ApplyQ( LEFT, NORMAL, A, tQ, Q );
59 bidiag::ApplyP( RIGHT, NORMAL, A, tP, P );
60 if( print )
61 {
62 Print( Q, "Q" );
63 Print( P, "P" );
64 }
65 if( display )
66 {
67 Display( Q, "Q" );
68 Display( P, "P" );
69 }
70 }
71
72 // Reverse the accumulated Householder transforms
73 bidiag::ApplyQ( LEFT, ADJOINT, A, tQ, AOrig );
74 bidiag::ApplyP( RIGHT, NORMAL, A, tP, AOrig );
75 if( print )
76 Print( AOrig, "Manual bidiagonal" );
77 if( display )
78 Display( AOrig, "Manual bidiagonal" );
79
80 // Compare the appropriate portion of AOrig and B
81 if( m >= n )
82 {
83 MakeTriangular( UPPER, AOrig );
84 MakeTrapezoidal( LOWER, AOrig, 1 );
85 }
86 else
87 {
88 MakeTriangular( LOWER, AOrig );
89 MakeTrapezoidal( UPPER, AOrig, -1 );
90 }
91 Axpy( F(-1), AOrig, B );
92 if( print )
93 Print( B, "Error in rotated bidiagonal" );
94 if( display )
95 Display( B, "Error in rotated bidiagonal" );
96 const Real infNormError = InfinityNorm( B );
97 const Real frobNormError = FrobeniusNorm( B );
98
99 if( g.Rank() == 0 )
100 {
101 cout << " ||A||_oo = " << infNormAOrig << "\n"
102 << " ||A||_F = " << frobNormAOrig << "\n"
103 << " ||B - Q^H A P||_oo = " << infNormError << "\n"
104 << " ||B - Q^H A P||_F = " << frobNormError << endl;
105 }
106 }
107
108 template<typename F>
TestBidiag(Int m,Int n,const Grid & g,bool testCorrectness,bool print,bool display)109 void TestBidiag
110 ( Int m, Int n, const Grid& g, bool testCorrectness, bool print, bool display )
111 {
112 DistMatrix<F> A(g), AOrig(g);
113 DistMatrix<F,STAR,STAR> tP(g), tQ(g);
114
115 Uniform( A, m, n );
116 if( testCorrectness )
117 {
118 if( g.Rank() == 0 )
119 {
120 cout << " Making copy of original matrix...";
121 cout.flush();
122 }
123 AOrig = A;
124 if( g.Rank() == 0 )
125 cout << "DONE" << endl;
126 }
127 if( print )
128 Print( A, "A" );
129 if( display )
130 Display( A, "A" );
131
132 if( g.Rank() == 0 )
133 {
134 cout << " Starting bidiagonalization...";
135 cout.flush();
136 }
137 mpi::Barrier( g.Comm() );
138 const double startTime = mpi::Time();
139 Bidiag( A, tP, tQ );
140 mpi::Barrier( g.Comm() );
141 const double runTime = mpi::Time() - startTime;
142 // TODO: Flop calculation
143 if( g.Rank() == 0 )
144 {
145 cout << "DONE. " << endl
146 << " Time = " << runTime << " seconds." << std::endl;
147 }
148 if( print )
149 {
150 Print( A, "A after Bidiag" );
151 Print( tP, "tP after Bidiag" );
152 Print( tQ, "tQ after Bidiag" );
153 }
154 if( display )
155 {
156 Display( A, "A after Bidiag" );
157 Display( tP, "tP after Bidiag" );
158 Display( tQ, "tQ after Bidiag" );
159 }
160 if( testCorrectness )
161 TestCorrectness( A, tP, tQ, AOrig, print, display );
162 }
163
164 int
main(int argc,char * argv[])165 main( int argc, char* argv[] )
166 {
167 Initialize( argc, argv );
168 mpi::Comm comm = mpi::COMM_WORLD;
169 const Int commRank = mpi::Rank( comm );
170 const Int commSize = mpi::Size( comm );
171
172 try
173 {
174 Int r = Input("--gridHeight","height of process grid",0);
175 const bool colMajor = Input("--colMajor","column-major ordering?",true);
176 const Int m = Input("--height","height of matrix",100);
177 const Int n = Input("--width","width of matrix",100);
178 const Int nb = Input("--nb","algorithmic blocksize",96);
179 const bool testCorrectness = Input
180 ("--correctness","test correctness?",true);
181 const bool print = Input("--print","print matrices?",false);
182 const bool display = Input("--display","display matrices?",false);
183 ProcessInput();
184 PrintInputReport();
185
186 if( r == 0 )
187 r = Grid::FindFactor( commSize );
188 const GridOrder order = ( colMajor ? COLUMN_MAJOR : ROW_MAJOR );
189 const Grid g( comm, r, order );
190 SetBlocksize( nb );
191 ComplainIfDebug();
192
193 if( commRank == 0 )
194 cout << "Double-precision:" << endl;
195 TestBidiag<double>( m, n, g, testCorrectness, print, display );
196
197 if( commRank == 0 )
198 cout << "Double-precision complex:" << endl;
199 TestBidiag<Complex<double>>( m, n, g, testCorrectness, print, display );
200 }
201 catch( exception& e ) { ReportException(e); }
202
203 Finalize();
204 return 0;
205 }
206