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 #pragma once
10 #ifndef ELEM_RIEMANN_HPP
11 #define ELEM_RIEMANN_HPP
12
13 namespace elem {
14
15 template<typename T>
16 inline void
Riemann(Matrix<T> & R,Int n)17 Riemann( Matrix<T>& R, Int n )
18 {
19 DEBUG_ONLY(CallStackEntry cse("Riemann"))
20 R.Resize( n, n );
21 for( Int j=0; j<n; ++j )
22 {
23 for( Int i=0; i<n; ++i )
24 {
25 if( ((j+2)%(i+2))==0 )
26 R.Set( i, j, T(i+1) );
27 else
28 R.Set( i, j, T(-1) );
29 }
30 }
31 }
32
33 template<typename T,Dist U,Dist V>
34 inline void
Riemann(DistMatrix<T,U,V> & R,Int n)35 Riemann( DistMatrix<T,U,V>& R, Int n )
36 {
37 DEBUG_ONLY(CallStackEntry cse("Riemann"))
38 R.Resize( n, n );
39 const Int localHeight = R.LocalHeight();
40 const Int localWidth = R.LocalWidth();
41 for( Int jLoc=0; jLoc<localWidth; ++jLoc )
42 {
43 const Int j = R.GlobalCol(jLoc);
44 for( Int iLoc=0; iLoc<localHeight; ++iLoc )
45 {
46 const Int i = R.GlobalRow(iLoc);
47 if( ((j+2)%(i+2))==0 )
48 R.SetLocal( iLoc, jLoc, T(i+1) );
49 else
50 R.SetLocal( iLoc, jLoc, T(-1) );
51 }
52 }
53 }
54
55 template<typename T,Dist U,Dist V>
56 inline void
Riemann(BlockDistMatrix<T,U,V> & R,Int n)57 Riemann( BlockDistMatrix<T,U,V>& R, Int n )
58 {
59 DEBUG_ONLY(CallStackEntry cse("Riemann"))
60 R.Resize( n, n );
61 const Int localHeight = R.LocalHeight();
62 const Int localWidth = R.LocalWidth();
63 for( Int jLoc=0; jLoc<localWidth; ++jLoc )
64 {
65 const Int j = R.GlobalCol(jLoc);
66 for( Int iLoc=0; iLoc<localHeight; ++iLoc )
67 {
68 const Int i = R.GlobalRow(iLoc);
69 if( ((j+2)%(i+2))==0 )
70 R.SetLocal( iLoc, jLoc, T(i+1) );
71 else
72 R.SetLocal( iLoc, jLoc, T(-1) );
73 }
74 }
75 }
76
77 } // namespace elem
78
79 #endif // ifndef ELEM_RIEMANN_HPP
80