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_MAKETRIANGULAR_INC
12 #include ELEM_UPDATEDIAGONAL_INC
13
14 #include ELEM_QR_INC
15 #include ELEM_FROBENIUSNORM_INC
16
17 #include ELEM_PERMUTECOLS_INC
18
19 #include ELEM_IDENTITY_INC
20 #include ELEM_UNIFORM_INC
21 using namespace std;
22 using namespace elem;
23
24 typedef double Real;
25 typedef Complex<Real> C;
26
27 int
main(int argc,char * argv[])28 main( int argc, char* argv[] )
29 {
30 Initialize( argc, argv );
31 const Int worldRank = mpi::WorldRank();
32
33 try
34 {
35 const Int m = Input("--height","height of matrix",100);
36 const Int n = Input("--width","width of matrix",100);
37 const bool alwaysRecompute = Input("--always","no norm updates?",false);
38 const bool blockedUnpiv =
39 Input("--blockUnpiv","blocked unpivoted QR?",false);
40 const bool display = Input("--display","display matrices?",false);
41 const bool print = Input("--print","print matrices?",false);
42 ProcessInput();
43 PrintInputReport();
44
45 DistMatrix<C> A;
46 Uniform( A, m, n );
47 const Real frobA = FrobeniusNorm( A );
48 if( display )
49 Display( A, "A" );
50 if( print )
51 Print( A, "A" );
52
53 // Compute the pivoted QR decomposition of A, but do not overwrite A
54 auto QRPiv( A );
55 DistMatrix<C,MD,STAR> tPiv;
56 DistMatrix<Real,MD,STAR> dPiv;
57 DistMatrix<Int,VR,STAR> perm;
58 qr::BusingerGolub( QRPiv, tPiv, dPiv, perm, alwaysRecompute );
59 if( display )
60 {
61 Display( QRPiv, "QRPiv" );
62 Display( tPiv, "tPiv" );
63 Display( dPiv, "dPiv" );
64 Display( perm, "perm" );
65 }
66 if( print )
67 {
68 Print( QRPiv, "QRPiv" );
69 Print( tPiv, "tPiv" );
70 Print( dPiv, "dPiv" );
71 Print( perm, "perm" );
72 }
73
74 // Compute the standard QR decomposition of A
75 auto QRNoPiv( A );
76 DistMatrix<C,MD,STAR> tNoPiv;
77 DistMatrix<Real,MD,STAR> dNoPiv;
78 if( blockedUnpiv )
79 QR( QRNoPiv, tNoPiv, dNoPiv );
80 else
81 qr::PanelHouseholder( QRNoPiv, tNoPiv, dNoPiv );
82 if( display )
83 {
84 Display( QRNoPiv, "QRNoPiv" );
85 Display( tNoPiv, "tNoPiv" );
86 Display( dNoPiv, "dNoPiv" );
87 }
88 if( print )
89 {
90 Print( QRNoPiv, "QRNoPiv" );
91 Print( tNoPiv, "tNoPiv" );
92 Print( dNoPiv, "dNoPiv" );
93 }
94
95 // Check the error in the pivoted QR factorization,
96 // || A P - Q R ||_F / || A ||_F
97 auto E( QRPiv );
98 MakeTriangular( UPPER, E );
99 qr::ApplyQ( LEFT, NORMAL, QRPiv, tPiv, dPiv, E );
100 PermuteCols( E, perm );
101 Axpy( C(-1), A, E );
102 const Real frobQRPiv = FrobeniusNorm( E );
103 if( display )
104 Display( E, "A P - Q R" );
105 if( print )
106 Print( E, "A P - Q R" );
107
108 // Check the error in the standard QR factorization,
109 // || A - Q R ||_F / || A ||_F
110 E = QRNoPiv;
111 MakeTriangular( UPPER, E );
112 qr::ApplyQ( LEFT, NORMAL, QRNoPiv, tNoPiv, dNoPiv, E );
113 Axpy( C(-1), A, E );
114 const Real frobQRNoPiv = FrobeniusNorm( E );
115 if( display )
116 Display( E, "A - Q R" );
117 if( print )
118 Print( E, "A - Q R" );
119
120 // Check orthogonality of pivoted Q, || I - Q^H Q ||_F / || A ||_F
121 Identity( E, m, n );
122 qr::ApplyQ( LEFT, NORMAL, QRPiv, tPiv, dPiv, E );
123 qr::ApplyQ( LEFT, ADJOINT, QRPiv, tPiv, dPiv, E );
124 const Int k = std::min(m,n);
125 auto EUpper = View( E, 0, 0, k, k );
126 UpdateDiagonal( EUpper, C(-1) );
127 const Real frobOrthogPiv = FrobeniusNorm( EUpper );
128 if( display )
129 Display( E, "pivoted I - Q^H Q" );
130 if( print )
131 Print( E, "pivoted I - Q^H Q" );
132
133 // Check orthogonality of unpivoted Q, || I - Q^H Q ||_F / || A ||_F
134 Identity( E, m, n );
135 qr::ApplyQ( LEFT, NORMAL, QRPiv, tPiv, dPiv, E );
136 qr::ApplyQ( LEFT, ADJOINT, QRPiv, tPiv, dPiv, E );
137 EUpper = View( E, 0, 0, k, k );
138 UpdateDiagonal( EUpper, C(-1) );
139 const Real frobOrthogNoPiv = FrobeniusNorm( EUpper );
140 if( display )
141 Display( E, "unpivoted I - Q^H Q" );
142 if( print )
143 Print( E, "unpivoted I - Q^H Q" );
144
145 if( worldRank == 0 )
146 {
147 std::cout << "|| A ||_F = " << frobA << "\n\n"
148 << "With pivoting: \n"
149 << " || A P - Q R ||_F / || A ||_F = "
150 << frobQRPiv/frobA << "\n"
151 << " || I - Q^H Q ||_F / || A ||_F = "
152 << frobOrthogPiv/frobA << "\n\n"
153 << "Without pivoting: \n"
154 << " || A - Q R ||_F / || A ||_F = "
155 << frobQRNoPiv/frobA << "\n"
156 << " || I - Q^H Q ||_F / || A ||_F = "
157 << frobOrthogNoPiv/frobA << "\n"
158 << std::endl;
159 }
160 }
161 catch( exception& e ) { ReportException(e); }
162
163 Finalize();
164 return 0;
165 }
166