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_CAUCHY_HPP
11 #define ELEM_CAUCHY_HPP
12
13 namespace elem {
14
15 template<typename F1,typename F2>
16 inline void
Cauchy(Matrix<F1> & A,const std::vector<F2> & x,const std::vector<F2> & y)17 Cauchy( Matrix<F1>& A, const std::vector<F2>& x, const std::vector<F2>& y )
18 {
19 DEBUG_ONLY(CallStackEntry cse("Cauchy"))
20 const Int m = x.size();
21 const Int n = y.size();
22 A.Resize( m, n );
23
24 const F1 one = F1(1);
25 for( Int j=0; j<n; ++j )
26 {
27 for( Int i=0; i<m; ++i )
28 {
29 DEBUG_ONLY(
30 // TODO: Use tolerance instead?
31 if( x[i] == y[j] )
32 LogicError
33 ( "x[", i, "] = y[", j, "] (", x[i],
34 ") is not allowed for Cauchy matrices" );
35 )
36 A.Set( i, j, one/(x[i]-y[j]) );
37 }
38 }
39 }
40
41 template<typename F1,typename F2,Dist U,Dist V>
42 inline void
Cauchy(DistMatrix<F1,U,V> & A,const std::vector<F2> & x,const std::vector<F2> & y)43 Cauchy
44 ( DistMatrix<F1,U,V>& A, const std::vector<F2>& x, const std::vector<F2>& y )
45 {
46 DEBUG_ONLY(CallStackEntry cse("Cauchy"))
47 const Int m = x.size();
48 const Int n = y.size();
49 A.Resize( m, n );
50
51 const F1 one = F1(1);
52 const Int localHeight = A.LocalHeight();
53 const Int localWidth = A.LocalWidth();
54 for( Int jLoc=0; jLoc<localWidth; ++jLoc )
55 {
56 const Int j = A.GlobalCol(jLoc);
57 for( Int iLoc=0; iLoc<localHeight; ++iLoc )
58 {
59 const Int i = A.GlobalRow(iLoc);
60 DEBUG_ONLY(
61 // TODO: Use tolerance instead?
62 if( x[i] == y[j] )
63 LogicError
64 ( "x[", i, "] = y[", j, "] (", x[i],
65 ") is not allowed for Cauchy matrices" );
66 )
67 A.SetLocal( iLoc, jLoc, one/(x[i]-y[j]) );
68 }
69 }
70 }
71
72 template<typename F1,typename F2,Dist U,Dist V>
73 inline void
Cauchy(BlockDistMatrix<F1,U,V> & A,const std::vector<F2> & x,const std::vector<F2> & y)74 Cauchy
75 ( BlockDistMatrix<F1,U,V>& A, const std::vector<F2>& x, const std::vector<F2>& y )
76 {
77 DEBUG_ONLY(CallStackEntry cse("Cauchy"))
78 const Int m = x.size();
79 const Int n = y.size();
80 A.Resize( m, n );
81
82 const F1 one = F1(1);
83 const Int localHeight = A.LocalHeight();
84 const Int localWidth = A.LocalWidth();
85 for( Int jLoc=0; jLoc<localWidth; ++jLoc )
86 {
87 const Int j = A.GlobalCol(jLoc);
88 for( Int iLoc=0; iLoc<localHeight; ++iLoc )
89 {
90 const Int i = A.GlobalRow(iLoc);
91 DEBUG_ONLY(
92 // TODO: Use tolerance instead?
93 if( x[i] == y[j] )
94 LogicError
95 ( "x[", i, "] = y[", j, "] (", x[i],
96 ") is not allowed for Cauchy matrices" );
97 )
98 A.SetLocal( iLoc, jLoc, one/(x[i]-y[j]) );
99 }
100 }
101 }
102
103 } // namespace elem
104
105 #endif // ifndef ELEM_CAUCHY_HPP
106