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_DIAGONALSCALE_INC
12 #include ELEM_HEMM_INC
13 #include ELEM_HERK_INC
14 #include ELEM_FROBENIUSNORM_INC
15 #include ELEM_IDENTITY_INC
16 using namespace std;
17 using namespace elem;
18
19 // Typedef our real and complex types to 'Real' and 'C' for convenience
20 typedef double Real;
21 typedef Complex<Real> C;
22
23 int
main(int argc,char * argv[])24 main( int argc, char* argv[] )
25 {
26 // This detects whether or not you have already initialized MPI and
27 // does so if necessary. The full routine is elem::Initialize.
28 Initialize( argc, argv );
29
30 // Surround the Elemental calls with try/catch statements in order to
31 // safely handle any exceptions that were thrown during execution.
32 try
33 {
34 const Int n = Input("--size","size of matrix",100);
35 const bool print = Input("--print","print matrices?",false);
36 ProcessInput();
37 PrintInputReport();
38
39 // Create a 2d process grid from a communicator. In our case, it is
40 // MPI_COMM_WORLD. There is another constructor that allows you to
41 // specify the grid dimensions, Grid g( comm, r ), which creates a
42 // grid of height r.
43 Grid g( mpi::COMM_WORLD );
44
45 // Create an n x n complex distributed matrix,
46 // We distribute the matrix using grid 'g'.
47 //
48 // There are quite a few available constructors, including ones that
49 // allow you to pass in your own local buffer and to specify the
50 // distribution alignments (i.e., which process row and column owns the
51 // top-left element)
52 DistMatrix<C> H( n, n, g );
53
54 // Fill the matrix since we did not pass in a buffer.
55 //
56 // We will fill entry (i,j) with the complex value (i+j,i-j) so that
57 // the global matrix is Hermitian. However, only one triangle of the
58 // matrix actually needs to be filled, the symmetry can be implicit.
59 //
60 const Int localHeight = H.LocalHeight();
61 const Int localWidth = H.LocalWidth();
62 for( Int jLoc=0; jLoc<localWidth; ++jLoc )
63 {
64 // Our process owns the rows colShift:colStride:n,
65 // and the columns rowShift:rowStride:n
66 const Int j = H.GlobalCol(jLoc);
67 for( Int iLoc=0; iLoc<localHeight; ++iLoc )
68 {
69 const Int i = H.GlobalRow(iLoc);
70 H.SetLocal( iLoc, jLoc, C(i+j,i-j) );
71 }
72 }
73
74 // Make a backup of H before we overwrite it within the eigensolver
75 auto HCopy( H );
76
77 // Call the eigensolver. We first create an empty complex eigenvector
78 // matrix, X[MC,MR], and an eigenvalue column vector, w[VR,* ]
79 //
80 // Optional: set blocksizes and algorithmic choices here. See the
81 // 'Tuning' section of the README for details.
82 DistMatrix<Real,VR,STAR> w( g );
83 DistMatrix<C> X( g );
84 HermitianEig( LOWER, H, w, X, ASCENDING );
85
86 if( print )
87 {
88 Print( HCopy, "H" );
89 Print( X, "X" );
90 Print( w, "w" );
91 }
92
93 // Check the residual, || H X - Omega X ||_F
94 const Real frobH = HermitianFrobeniusNorm( LOWER, HCopy );
95 auto E( X );
96 DiagonalScale( RIGHT, NORMAL, w, E );
97 Hemm( LEFT, LOWER, C(-1), HCopy, X, C(1), E );
98 const Real frobResid = FrobeniusNorm( E );
99
100 // Check the orthogonality of X
101 Identity( E, n, n );
102 Herk( LOWER, NORMAL, C(-1), X, C(1), E );
103 const Real frobOrthog = HermitianFrobeniusNorm( LOWER, E );
104
105 if( mpi::WorldRank() == 0 )
106 {
107 std::cout << "|| H ||_F = " << frobH << "\n"
108 << "|| H X - X Omega ||_F / || A ||_F = "
109 << frobResid / frobH << "\n"
110 << "|| X X^H - I ||_F = " << frobOrthog / frobH
111 << "\n" << std::endl;
112 }
113 }
114 catch( exception& e ) { ReportException(e); }
115
116 Finalize();
117 return 0;
118 }
119