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