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_FROBENIUSNORM_INC
12 #include ELEM_GLM_INC
13 #include ELEM_UNIFORM_INC
14 using namespace elem;
15
16 typedef double Real;
17 typedef Complex<Real> F;
18
19 int
main(int argc,char * argv[])20 main( int argc, char* argv[] )
21 {
22 Initialize( argc, argv );
23 mpi::Comm comm = mpi::COMM_WORLD;
24 const Int commRank = mpi::Rank( comm );
25 const Int commSize = mpi::Size( comm );
26
27 try
28 {
29 const Int m = Input("--m","height of A",100);
30 const Int n = Input("--n","width of A",100);
31 const Int p = Input("--p","width of B",20);
32 const Int numRhs = Input("--numRhs","# of right-hand sides",5);
33 const Int blocksize = Input("--blocksize","algorithmic blocksize",64);
34 const bool print = Input("--print","print matrices?",false);
35 Int gridHeight = Input("--gridHeight","grid height",0);
36 ProcessInput();
37 PrintInputReport();
38
39 // Set the algorithmic blocksize
40 SetBlocksize( blocksize );
41
42 // If the grid height wasn't specified, then we should attempt to build
43 // a nearly-square process grid
44 if( gridHeight == 0 )
45 gridHeight = Grid::FindFactor( commSize );
46 Grid g( comm, gridHeight );
47
48 DistMatrix<F> A(g), B(g), D(g), Y(g);
49 Uniform( A, m, n );
50 Uniform( B, m, p );
51 Uniform( D, m, numRhs );
52 auto ACopy( A );
53 auto BCopy( B );
54 auto DCopy( D );
55
56 if( print )
57 {
58 Print( A, "A" );
59 Print( B, "B" );
60 Print( D, "D" );
61 }
62
63 GLM( A, B, D, Y );
64
65 if( print )
66 {
67 Print( D, "X" );
68 Print( Y, "Y" );
69 }
70
71 const Real DFrob = FrobeniusNorm( DCopy );
72 Gemm( NORMAL, NORMAL, F(-1), ACopy, D, F(1), DCopy );
73 Gemm( NORMAL, NORMAL, F(-1), BCopy, Y, F(1), DCopy );
74 const Real EFrob = FrobeniusNorm( DCopy );
75 if( print )
76 Print( DCopy, "D - A X - B Y" );
77 if( commRank == 0 )
78 std::cout << "|| D ||_F = " << DFrob << "\n"
79 << "|| A X + B Y - D ||_F = " << EFrob << "\n"
80 << std::endl;
81 }
82 catch( std::exception& e ) { ReportException(e); }
83
84 Finalize();
85 return 0;
86 }
87