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