本文整理汇总了C++中Axpy函数的典型用法代码示例。如果您正苦于以下问题：C++ Axpy函数的具体用法？C++ Axpy怎么用？C++ Axpy使用的例子？那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了Axpy函数的20个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于我们的系统推荐出更棒的C++代码示例。
示例1: MakeExplicitlyHermitian
void MakeExplicitlyHermitian( UpperOrLower uplo, DistMatrix<F,MC,MR>& A )
{
const Grid& g = A.Grid();
DistMatrix<F,MC,MR> ATL(g), ATR(g), A00(g), A01(g), A02(g),
ABL(g), ABR(g), A10(g), A11(g), A12(g),
A20(g), A21(g), A22(g);
DistMatrix<F,MC,MR> A11Adj(g);
DistMatrix<F,MR,MC> A11_MR_MC(g);
DistMatrix<F,MR,MC> A21_MR_MC(g);
DistMatrix<F,MR,MC> A12_MR_MC(g);
PartitionDownDiagonal
( A, ATL, ATR,
ABL, ABR, 0 );
while( ATL.Height() < A.Height() )
{
RepartitionDownDiagonal
( ATL, /**/ ATR, A00, /**/ A01, A02,
/*************/ /******************/
/**/ A10, /**/ A11, A12,
ABL, /**/ ABR, A20, /**/ A21, A22 );
A11Adj.AlignWith( A11 );
A11_MR_MC.AlignWith( A11 );
A12_MR_MC.AlignWith( A21 );
A21_MR_MC.AlignWith( A12 );
//--------------------------------------------------------------------//
A11_MR_MC = A11;
A11Adj.ResizeTo( A11.Height(), A11.Width() );
Adjoint( A11_MR_MC.LocalMatrix(), A11Adj.LocalMatrix() );
if( uplo == LOWER )
{
MakeTrapezoidal( LEFT, UPPER, 1, A11Adj );
Axpy( (F)1, A11Adj, A11 );
A21_MR_MC = A21;
Adjoint( A21_MR_MC.LocalMatrix(), A12.LocalMatrix() );
}
else
{
MakeTrapezoidal( LEFT, LOWER, -1, A11Adj );
Axpy( (F)1, A11Adj, A11 );
A12_MR_MC = A12;
Adjoint( A12_MR_MC.LocalMatrix(), A21.LocalMatrix() );
}
//--------------------------------------------------------------------//
A21_MR_MC.FreeAlignments();
A12_MR_MC.FreeAlignments();
A11_MR_MC.FreeAlignments();
A11Adj.FreeAlignments();
SlidePartitionDownDiagonal
( ATL, /**/ ATR, A00, A01, /**/ A02,
/**/ A10, A11, /**/ A12,
/*************/ /******************/
ABL, /**/ ABR, A20, A21, /**/ A22 );
}
}
开发者ID:ahmadia，项目名称:elemental，代码行数:60，代码来源:HPSDCholesky.hpp
示例2: a1_like_a2
void Transform2x2
( const Matrix<T>& G,
AbstractDistMatrix<T>& a1,
AbstractDistMatrix<T>& a2 )
{
DEBUG_CSE
typedef unique_ptr<AbstractDistMatrix<T>> ADMPtr;
// TODO: Optimize by attempting SendRecv when possible
ADMPtr a1_like_a2( a2.Construct( a2.Grid(), a2.Root() ) );
a1_like_a2->AlignWith( DistData(a2) );
Copy( a1, *a1_like_a2 );
ADMPtr a2_like_a1( a1.Construct( a1.Grid(), a1.Root() ) );
a2_like_a1->AlignWith( DistData(a1) );
Copy( a2, *a2_like_a1 );
// TODO: Generalized axpy?
Scale( G(0,0), a1 );
Axpy( G(0,1), *a2_like_a1, a1 );
// TODO: Generalized axpy?
Scale( G(1,1), a2 );
Axpy( G(1,0), *a1_like_a2, a2 );
}
开发者ID:jeffhammond，项目名称:Elemental，代码行数:26，代码来源:Transform2x2.cpp
示例3: Copy
/* Prototype implementation for specialized functions */
void Vector::AddTwoVectorsImpl(Number a, const Vector& v1,
Number b, const Vector& v2, Number c)
{
if (c==0.) {
if (a==1.) {
Copy(v1);
if (b!=0.) {
Axpy(b, v2);
}
}
else if (a==0.) {
if (b==0.) {
Set(0.);
}
else {
Copy(v2);
if (b!=1.) {
Scal(b);
}
}
}
else {
if (b==1.) {
Copy(v2);
Axpy(a, v1);
}
else if (b==0.) {
Copy(v1);
Scal(a);
}
else {
Copy(v1);
Scal(a);
Axpy(b, v2);
}
}
}
else { /* c==0. */
if (c!=1.) {
Scal(c);
}
if (a!=0.) {
Axpy(a, v1);
}
if (b!=0.) {
Axpy(b, v2);
}
}
}
开发者ID:Gjacquenot，项目名称:simbody，代码行数:50，代码来源:IpVector.cpp
示例4: entry
inline void
NewtonStep
( const DistMatrix<F>& X, DistMatrix<F>& XNew, Scaling scaling=FROB_NORM )
{
#ifndef RELEASE
CallStackEntry entry("sign::NewtonStep");
#endif
typedef BASE(F) Real;
// Calculate mu while forming B := inv(X)
Real mu;
DistMatrix<Int,VC,STAR> p( X.Grid() );
XNew = X;
LU( XNew, p );
if( scaling == DETERMINANT )
{
SafeProduct<F> det = determinant::AfterLUPartialPiv( XNew, p );
mu = Real(1)/Exp(det.kappa);
}
inverse::AfterLUPartialPiv( XNew, p );
if( scaling == FROB_NORM )
mu = Sqrt( FrobeniusNorm(XNew)/FrobeniusNorm(X) );
else if( scaling == NONE )
mu = 1;
else
LogicError("Scaling case not handled");
// Overwrite XNew with the new iterate
const Real halfMu = mu/Real(2);
const Real halfMuInv = Real(1)/(2*mu);
Scale( halfMuInv, XNew );
Axpy( halfMu, X, XNew );
}
开发者ID:khalid-hasanov，项目名称:Elemental，代码行数:33，代码来源:Sign.hpp
示例5: P
void
NewtonStep
( const DistMatrix<Field>& X,
DistMatrix<Field>& XNew,
SignScaling scaling=SIGN_SCALE_FROB )
{
EL_DEBUG_CSE
typedef Base<Field> Real;
// Calculate mu while forming B := inv(X)
Real mu=1;
DistPermutation P( X.Grid() );
XNew = X;
LU( XNew, P );
if( scaling == SIGN_SCALE_DET )
{
SafeProduct<Field> det = det::AfterLUPartialPiv( XNew, P );
mu = Real(1)/Exp(det.kappa);
}
inverse::AfterLUPartialPiv( XNew, P );
if( scaling == SIGN_SCALE_FROB )
mu = Sqrt( FrobeniusNorm(XNew)/FrobeniusNorm(X) );
// Overwrite XNew with the new iterate
const Real halfMu = mu/Real(2);
const Real halfMuInv = Real(1)/(2*mu);
XNew *= halfMuInv;
Axpy( halfMu, X, XNew );
}
开发者ID:elemental，项目名称:Elemental，代码行数:29，代码来源:Sign.cpp
示例6: DEBUG_ONLY
const BlockCyclicMatrix<T>&
BlockCyclicMatrix<T>::operator-=( const BlockCyclicMatrix<T>& A )
{
DEBUG_ONLY(CSE cse("BCM::operator-="))
Axpy( T(-1), A, *this );
return *this;
}
开发者ID:bluehope，项目名称:Elemental，代码行数:7，代码来源:BlockCyclic.cpp
示例7: MakeSymmetric
inline void
MakeSymmetric( UpperOrLower uplo, DistMatrix<T>& A )
{
#ifndef RELEASE
PushCallStack("MakeSymmetric");
#endif
if( A.Height() != A.Width() )
throw std::logic_error("Cannot make non-square matrix symmetric");
const Grid& g = A.Grid();
DistMatrix<T,MD,STAR> d(g);
A.GetDiagonal( d );
if( uplo == LOWER )
MakeTrapezoidal( LEFT, LOWER, -1, A );
else
MakeTrapezoidal( LEFT, UPPER, +1, A );
DistMatrix<T> ATrans(g);
Transpose( A, ATrans );
Axpy( T(1), ATrans, A );
A.SetDiagonal( d );
#ifndef RELEASE
PopCallStack();
#endif
}
开发者ID:jimgoo，项目名称:Elemental，代码行数:26，代码来源:MakeSymmetric.hpp
示例8: HermitianTridiagU
inline void HermitianTridiagU( Matrix<R>& A )
{
#ifndef RELEASE
PushCallStack("HermitianTridiagU");
if( A.Height() != A.Width() )
throw std::logic_error( "A must be square." );
#endif
// Matrix views
Matrix<R>
ATL, ATR, A00, a01, A02, a01T,
ABL, ABR, a10, alpha11, a12, alpha01B,
A20, a21, A22;
// Temporary matrices
Matrix<R> w01;
PushBlocksizeStack( 1 );
PartitionUpDiagonal
( A, ATL, ATR,
ABL, ABR, 0 );
while( ABR.Height()+1 < A.Height() )
{
RepartitionUpDiagonal
( ATL, /**/ ATR, A00, a01, /**/ A02,
/**/ a10, alpha11, /**/ a12,
/*************/ /**********************/
ABL, /**/ ABR, A20, a21, /**/ A22 );
PartitionUp
( a01, a01T,
alpha01B, 1 );
w01.ResizeTo( a01.Height(), 1 );
//--------------------------------------------------------------------//
const R tau = Reflector( alpha01B, a01T );
const R epsilon1 = alpha01B.Get(0,0);
alpha01B.Set(0,0,R(1));
Symv( UPPER, tau, A00, a01, R(0), w01 );
const R alpha = -tau*Dot( w01, a01 )/R(2);
Axpy( alpha, a01, w01 );
Syr2( UPPER, R(-1), a01, w01, A00 );
alpha01B.Set(0,0,epsilon1);
//--------------------------------------------------------------------//
SlidePartitionUpDiagonal
( ATL, /**/ ATR, A00, /**/ a01, A02,
/*************/ /**********************/
/**/ a10, /**/ alpha11, a12,
ABL, /**/ ABR, A20, /**/ a21, A22 );
}
PopBlocksizeStack();
#ifndef RELEASE
PopCallStack();
#endif
}
开发者ID:certik，项目名称:Elemental，代码行数:56，代码来源:Local.hpp
示例9: L
void L( Matrix<F>& A, Matrix<F>& t )
{
#ifndef RELEASE
CallStackEntry entry("hermitian_tridiag::L");
if( A.Height() != A.Width() )
LogicError("A must be square");
#endif
typedef BASE(F) R;
const Int tHeight = Max(A.Height()-1,0);
t.ResizeTo( tHeight, 1 );
// Matrix views
Matrix<F>
ATL, ATR, A00, a01, A02, alpha21T,
ABL, ABR, a10, alpha11, a12, a21B,
A20, a21, A22;
// Temporary matrices
Matrix<F> w21;
PartitionDownDiagonal
( A, ATL, ATR,
ABL, ABR, 0 );
while( ATL.Height()+1 < A.Height() )
{
RepartitionDownDiagonal
( ATL, /**/ ATR, A00, /**/ a01, A02,
/*************/ /**********************/
/**/ a10, /**/ alpha11, a12,
ABL, /**/ ABR, A20, /**/ a21, A22, 1 );
PartitionDown
( a21, alpha21T,
a21B, 1 );
//--------------------------------------------------------------------//
const F tau = Reflector( alpha21T, a21B );
const R epsilon1 = alpha21T.GetRealPart(0,0);
t.Set(A00.Height(),0,tau);
alpha21T.Set(0,0,F(1));
Zeros( w21, a21.Height(), 1 );
Hemv( LOWER, tau, A22, a21, F(0), w21 );
const F alpha = -tau*Dot( w21, a21 )/F(2);
Axpy( alpha, a21, w21 );
Her2( LOWER, F(-1), a21, w21, A22 );
alpha21T.Set(0,0,epsilon1);
//--------------------------------------------------------------------//
SlidePartitionDownDiagonal
( ATL, /**/ ATR, A00, a01, /**/ A02,
/**/ a10, alpha11, /**/ a12,
/*************/ /**********************/
ABL, /**/ ABR, A20, a21, /**/ A22 );
}
}
开发者ID:khalid-hasanov，项目名称:Elemental，代码行数:56，代码来源:L.hpp
示例10: RunRoutine
// Describes how to run the CLBlast routine
static StatusCode RunRoutine(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
auto queue_plain = queue();
auto event = cl_event{};
auto status = Axpy(args.n, args.alpha,
buffers.x_vec(), args.x_offset, args.x_inc,
buffers.y_vec(), args.y_offset, args.y_inc,
&queue_plain, &event);
if (status == StatusCode::kSuccess) { clWaitForEvents(1, &event); clReleaseEvent(event); }
return status;
}
开发者ID:dividiti，项目名称:CLBlast，代码行数:11，代码来源:xaxpy.hpp
示例11: entry
inline void
NewtonStep
( const Matrix<F>& A, const Matrix<F>& X, Matrix<F>& XNew, Matrix<F>& XTmp )
{
#ifndef RELEASE
CallStackEntry entry("square_root::NewtonStep");
#endif
// XNew := inv(X) A
XTmp = X;
Matrix<Int> p;
LU( XTmp, p );
XNew = A;
lu::SolveAfter( NORMAL, XTmp, p, XNew );
// XNew := 1/2 ( X + XNew )
typedef BASE(F) R;
Axpy( R(1)/R(2), X, XNew );
}
开发者ID:khalid-hasanov，项目名称:Elemental，代码行数:18，代码来源:SquareRoot.hpp
示例12: Blocksize
void SUMMA_NTB
( Orientation orientB,
T alpha,
const AbstractDistMatrix<T>& APre,
const AbstractDistMatrix<T>& BPre,
AbstractDistMatrix<T>& CPre )
{
EL_DEBUG_CSE
const Int m = CPre.Height();
const Int bsize = Blocksize();
const Grid& g = APre.Grid();
DistMatrixReadProxy<T,T,MC,MR> AProx( APre );
DistMatrixReadProxy<T,T,MC,MR> BProx( BPre );
DistMatrixReadWriteProxy<T,T,MC,MR> CProx( CPre );
auto& A = AProx.GetLocked();
auto& B = BProx.GetLocked();
auto& C = CProx.Get();
// Temporary distributions
DistMatrix<T,MR,STAR> A1Trans_MR_STAR(g);
DistMatrix<T,STAR,MC> D1_STAR_MC(g);
DistMatrix<T,MR,MC> D1_MR_MC(g);
A1Trans_MR_STAR.AlignWith( B );
D1_STAR_MC.AlignWith( B );
for( Int k=0; k<m; k+=bsize )
{
const Int nb = Min(bsize,m-k);
auto A1 = A( IR(k,k+nb), ALL );
auto C1 = C( IR(k,k+nb), ALL );
// D1[*,MC] := alpha A1[*,MR] (B[MC,MR])^T
// = alpha (A1^T)[MR,*] (B^T)[MR,MC]
Transpose( A1, A1Trans_MR_STAR );
LocalGemm( TRANSPOSE, orientB, alpha, A1Trans_MR_STAR, B, D1_STAR_MC );
// C1[MC,MR] += scattered & transposed D1[*,MC] summed over grid rows
Contract( D1_STAR_MC, D1_MR_MC );
Axpy( T(1), D1_MR_MC, C1 );
}
}
开发者ID:elemental，项目名称:Elemental，代码行数:43，代码来源:NT.hpp
示例13: Blocksize
void SUMMA_TNA
( Orientation orientA,
T alpha,
const AbstractDistMatrix<T>& APre,
const AbstractDistMatrix<T>& BPre,
AbstractDistMatrix<T>& CPre )
{
DEBUG_CSE
const Int n = CPre.Width();
const Int bsize = Blocksize();
const Grid& g = APre.Grid();
DistMatrixReadProxy<T,T,MC,MR> AProx( APre );
DistMatrixReadProxy<T,T,MC,MR> BProx( BPre );
DistMatrixReadWriteProxy<T,T,MC,MR> CProx( CPre );
auto& A = AProx.GetLocked();
auto& B = BProx.GetLocked();
auto& C = CProx.Get();
// Temporary distributions
DistMatrix<T,MC,STAR> B1_MC_STAR(g);
DistMatrix<T,MR,STAR> D1_MR_STAR(g);
DistMatrix<T,MR,MC > D1_MR_MC(g);
B1_MC_STAR.AlignWith( A );
D1_MR_STAR.AlignWith( A );
for( Int k=0; k<n; k+=bsize )
{
const Int nb = Min(bsize,n-k);
auto B1 = B( ALL, IR(k,k+nb) );
auto C1 = C( ALL, IR(k,k+nb) );
// D1[MR,*] := alpha (A1[MC,MR])^T B1[MC,*]
// = alpha (A1^T)[MR,MC] B1[MC,*]
B1_MC_STAR = B1;
LocalGemm( orientA, NORMAL, alpha, A, B1_MC_STAR, D1_MR_STAR );
// C1[MC,MR] += scattered & transposed D1[MR,*] summed over grid cols
Contract( D1_MR_STAR, D1_MR_MC );
Axpy( T(1), D1_MR_MC, C1 );
}
}
开发者ID:YingzhouLi，项目名称:Elemental，代码行数:43，代码来源:TN.hpp
示例14: Newton
Int
Newton( DistMatrix<Field>& A, const SignCtrl<Base<Field>>& ctrl )
{
EL_DEBUG_CSE
typedef Base<Field> Real;
Real tol = ctrl.tol;
if( tol == Real(0) )
tol = A.Height()*limits::Epsilon<Real>();
Int numIts=0;
DistMatrix<Field> B( A.Grid() );
DistMatrix<Field> *X=&A, *XNew=&B;
while( numIts < ctrl.maxIts )
{
// Overwrite XNew with the new iterate
NewtonStep( *X, *XNew, ctrl.scaling );
// Use the difference in the iterates to test for convergence
Axpy( Real(-1), *XNew, *X );
const Real oneDiff = OneNorm( *X );
const Real oneNew = OneNorm( *XNew );
// Ensure that X holds the current iterate and break if possible
++numIts;
std::swap( X, XNew );
if( ctrl.progress && A.Grid().Rank() == 0 )
cout << "after " << numIts << " Newton iter's: "
<< "oneDiff=" << oneDiff << ", oneNew=" << oneNew
<< ", oneDiff/oneNew=" << oneDiff/oneNew << ", tol="
<< tol << endl;
if( oneDiff/oneNew <= Pow(oneNew,ctrl.power)*tol )
break;
}
if( X != &A )
A = *X;
return numIts;
}
开发者ID:elemental，项目名称:Elemental，代码行数:37，代码来源:Sign.cpp
示例15: entry
inline void
Symv
( UpperOrLower uplo,
T alpha, const DistMatrix<T>& A,
const DistMatrix<T>& x,
T beta, DistMatrix<T>& y,
bool conjugate=false )
{
#ifndef RELEASE
CallStackEntry entry("Symv");
if( A.Grid() != x.Grid() || x.Grid() != y.Grid() )
throw std::logic_error
("{A,x,y} must be distributed over the same grid");
if( A.Height() != A.Width() )
throw std::logic_error("A must be square");
if( ( x.Width() != 1 && x.Height() != 1 ) ||
( y.Width() != 1 && y.Height() != 1 ) )
throw std::logic_error("x and y are assumed to be vectors");
const int xLength = ( x.Width()==1 ? x.Height() : x.Width() );
const int yLength = ( y.Width()==1 ? y.Height() : y.Width() );
if( A.Height() != xLength || A.Height() != yLength )
{
std::ostringstream msg;
msg << "Nonconformal Symv: \n"
<< " A ~ " << A.Height() << " x " << A.Width() << "\n"
<< " x ~ " << x.Height() << " x " << x.Width() << "\n"
<< " y ~ " << y.Height() << " x " << y.Width() << "\n";
throw std::logic_error( msg.str() );
}
#endif
const Grid& g = A.Grid();
if( x.Width() == 1 && y.Width() == 1 )
{
// Temporary distributions
DistMatrix<T,MC,STAR> x_MC_STAR(g), z_MC_STAR(g);
DistMatrix<T,MR,STAR> x_MR_STAR(g), z_MR_STAR(g);
DistMatrix<T,MR,MC > z_MR_MC(g);
DistMatrix<T> z(g);
// Begin the algoritm
Scale( beta, y );
x_MC_STAR.AlignWith( A );
x_MR_STAR.AlignWith( A );
z_MC_STAR.AlignWith( A );
z_MR_STAR.AlignWith( A );
z.AlignWith( y );
Zeros( z_MC_STAR, y.Height(), 1 );
Zeros( z_MR_STAR, y.Height(), 1 );
//--------------------------------------------------------------------//
x_MC_STAR = x;
x_MR_STAR = x_MC_STAR;
if( uplo == LOWER )
{
internal::LocalSymvColAccumulateL
( alpha, A, x_MC_STAR, x_MR_STAR, z_MC_STAR, z_MR_STAR, conjugate );
}
else
{
internal::LocalSymvColAccumulateU
( alpha, A, x_MC_STAR, x_MR_STAR, z_MC_STAR, z_MR_STAR, conjugate );
}
z_MR_MC.SumScatterFrom( z_MR_STAR );
z = z_MR_MC;
z.SumScatterUpdate( T(1), z_MC_STAR );
Axpy( T(1), z, y );
//--------------------------------------------------------------------//
x_MC_STAR.FreeAlignments();
x_MR_STAR.FreeAlignments();
z_MC_STAR.FreeAlignments();
z_MR_STAR.FreeAlignments();
z.FreeAlignments();
}
else if( x.Width() == 1 )
{
// Temporary distributions
DistMatrix<T,MC,STAR> x_MC_STAR(g), z_MC_STAR(g);
DistMatrix<T,MR,STAR> x_MR_STAR(g), z_MR_STAR(g);
DistMatrix<T,MR,MC > z_MR_MC(g);
DistMatrix<T> z(g), zTrans(g);
// Begin the algoritm
Scale( beta, y );
x_MC_STAR.AlignWith( A );
x_MR_STAR.AlignWith( A );
z_MC_STAR.AlignWith( A );
z_MR_STAR.AlignWith( A );
z.AlignWith( y );
z_MR_MC.AlignWith( y );
Zeros( z_MC_STAR, y.Width(), 1 );
Zeros( z_MR_STAR, y.Width(), 1 );
//--------------------------------------------------------------------//
x_MC_STAR = x;
x_MR_STAR = x_MC_STAR;
if( uplo == LOWER )
{
internal::LocalSymvColAccumulateL
( alpha, A, x_MC_STAR, x_MR_STAR, z_MC_STAR, z_MR_STAR, conjugate );
}
//.........这里部分代码省略.........
开发者ID:ahmadia，项目名称:Elemental-1，代码行数:101，代码来源:Symv.hpp
示例16: Min
void LUMod
( Matrix<F>& A,
Permutation& P,
const Matrix<F>& u,
const Matrix<F>& v,
bool conjugate,
Base<F> tau )
{
DEBUG_CSE
typedef Base<F> Real;
const Int m = A.Height();
const Int n = A.Width();
const Int minDim = Min(m,n);
if( minDim != m )
LogicError("It is assumed that height(A) <= width(A)");
if( u.Height() != m || u.Width() != 1 )
LogicError("u is expected to be a conforming column vector");
if( v.Height() != n || v.Width() != 1 )
LogicError("v is expected to be a conforming column vector");
// w := inv(L) P u
auto w( u );
P.PermuteRows( w );
Trsv( LOWER, NORMAL, UNIT, A, w );
// Maintain an external vector for the temporary subdiagonal of U
Matrix<F> uSub;
Zeros( uSub, minDim-1, 1 );
// Reduce w to a multiple of e0
for( Int i=minDim-2; i>=0; --i )
{
// Decide if we should pivot the i'th and i+1'th rows of w
const F lambdaSub = A(i+1,i);
const F ups_ii = A(i,i);
const F omega_i = w(i);
const F omega_ip1 = w(i+1);
const Real rightTerm = Abs(lambdaSub*omega_i+omega_ip1);
const bool pivot = ( Abs(omega_i) < tau*rightTerm );
const Range<Int> indi( i, i+1 ),
indip1( i+1, i+2 ),
indB( i+2, m ),
indR( i+1, n );
auto lBi = A( indB, indi );
auto lBip1 = A( indB, indip1 );
auto uiR = A( indi, indR );
auto uip1R = A( indip1, indR );
if( pivot )
{
// P := P_i P
P.Swap( i, i+1 );
// Simultaneously perform
// U := P_i U and
// L := P_i L P_i^T
//
// Then update
// L := L T_{i,L}^{-1},
// U := T_{i,L} U,
// w := T_{i,L} P_i w,
// where T_{i,L} is the Gauss transform which zeros (P_i w)_{i+1}.
//
// More succinctly,
// gamma := w(i) / w(i+1),
// w(i) := w(i+1),
// w(i+1) := 0,
// L(:,i) += gamma L(:,i+1),
// U(i+1,:) -= gamma U(i,:).
const F gamma = omega_i / omega_ip1;
const F lambda_ii = F(1) + gamma*lambdaSub;
A(i, i) = gamma;
A(i+1,i) = 0;
auto lBiCopy = lBi;
Swap( NORMAL, lBi, lBip1 );
Axpy( gamma, lBiCopy, lBi );
auto uip1RCopy = uip1R;
RowSwap( A, i, i+1 );
Axpy( -gamma, uip1RCopy, uip1R );
// Force L back to *unit* lower-triangular form via the transform
// L := L T_{i,U}^{-1} D^{-1},
// where D is diagonal and responsible for forcing L(i,i) and
// L(i+1,i+1) back to 1. The effect on L is:
// eta := L(i,i+1)/L(i,i),
// L(:,i+1) -= eta L(:,i),
// delta_i := L(i,i),
// delta_ip1 := L(i+1,i+1),
// L(:,i) /= delta_i,
// L(:,i+1) /= delta_ip1,
// while the effect on U is
// U(i,:) += eta U(i+1,:)
// U(i,:) *= delta_i,
// U(i+1,:) *= delta_{i+1},
// and the effect on w is
// w(i) *= delta_i.
//.........这里部分代码省略.........
开发者ID:YingzhouLi，项目名称:Elemental，代码行数:101，代码来源:Mod.hpp
示例17: TwoSidedTrmmUVar5
inline void
TwoSidedTrmmUVar5( UnitOrNonUnit diag, Matrix<F>& A, const Matrix<F>& U )
{
#ifndef RELEASE
PushCallStack("internal::TwoSidedTrmmUVar5");
if( A.Height() != A.Width() )
throw std::logic_error("A must be square");
if( U.Height() != U.Width() )
throw std::logic_error("Triangular matrices must be square");
if( A.Height() != U.Height() )
throw std::logic_error("A and U must be the same size");
#endif
// Matrix views
Matrix<F>
ATL, ATR, A00, A01, A02,
ABL, ABR, A10, A11, A12,
A20, A21, A22;
Matrix<F>
UTL, UTR, U00, U01, U02,
UBL, UBR, U10, U11, U12,
U20, U21, U22;
// Temporary products
Matrix<F> Y01;
PartitionDownDiagonal
( A, ATL, ATR,
ABL, ABR, 0 );
LockedPartitionDownDiagonal
( U, UTL, UTR,
UBL, UBR, 0 );
while( ATL.Height() < A.Height() )
{
RepartitionDownDiagonal
( ATL, /**/ ATR, A00, /**/ A01, A02,
/*************/ /******************/
/**/ A10, /**/ A11, A12,
ABL, /**/ ABR, A20, /**/ A21, A22 );
LockedRepartitionDownDiagonal
( UTL, /**/ UTR, U00, /**/ U01, U02,
/*************/ /******************/
/**/ U10, /**/ U11, U12,
UBL, /**/ UBR, U20, /**/ U21, U22 );
//--------------------------------------------------------------------//
// Y01 := U01 A11
Zeros( A01.Height(), A01.Width(), Y01 );
Hemm( RIGHT, UPPER, F(1), A11, U01, F(0), Y01 );
// A01 := U00 A01
Trmm( LEFT, UPPER, NORMAL, diag, F(1), U00, A01 );
// A01 := A01 + 1/2 Y01
Axpy( F(1)/F(2), Y01, A01 );
// A00 := A00 + (U01 A01' + A01 U01')
Her2k( UPPER, NORMAL, F(1), U01, A01, F(1), A00 );
// A01 := A01 + 1/2 Y01
Axpy( F(1)/F(2), Y01, A01 );
// A01 := A01 U11'
Trmm( RIGHT, UPPER, ADJOINT, diag, F(1), U11, A01 );
// A11 := U11 A11 U11'
TwoSidedTrmmUUnb( diag, A11, U11 );
//--------------------------------------------------------------------//
SlidePartitionDownDiagonal
( ATL, /**/ ATR, A00, A01, /**/ A02,
/**/ A10, A11, /**/ A12,
/*************/ /******************/
ABL, /**/ ABR, A20, A21, /**/ A22 );
SlideLockedPartitionDownDiagonal
( UTL, /**/ UTR, U00, U01, /**/ U02,
/**/ U10, U11, /**/ U12,
/*************/ /******************/
UBL, /**/ UBR, U20, U21, /**/ U22 );
}
#ifndef RELEASE
PopCallStack();
#endif
}
开发者ID:jimgoo，项目名称:Elemental，代码行数:85，代码来源:UVar5.hpp
示例18: PushCallStack
inline void
TwoSidedTrmmUVar5
( UnitOrNonUnit diag, DistMatrix<F>& A, const DistMatrix<F>& U )
{
#ifndef RELEASE
PushCallStack("internal::TwoSidedTrmmUVar5");
if( A.Height() != A.Width() )
throw std::logic_error("A must be square");
if( U.Height() != U.Width() )
throw std::logic_error("Triangular matrices must be square");
if( A.Height() != U.Height() )
throw std::logic_error("A and U must be the same size");
#endif
const Grid& g = A.Grid();
// Matrix views
DistMatrix<F>
ATL(g), ATR(g), A00(g), A01(g), A02(g),
ABL(g), ABR(g), A10(g), A11(g), A12(g),
A20(g), A21(g), A22(g);
DistMatrix<F>
UTL(g), UTR(g), U00(g), U01(g), U02(g),
UBL(g), UBR(g), U10(g), U11(g), U12(g),
U20(g), U21(g), U22(g);
// Temporary distributions
DistMatrix<F,STAR,STAR> A11_STAR_STAR(g);
DistMatrix<F,MC, STAR> A01_MC_STAR(g);
DistMatrix<F,MR, STAR> A01_MR_STAR(g);
DistMatrix<F,VC, STAR> A01_VC_STAR(g);
DistMatrix<F,STAR,STAR> U11_STAR_STAR(g);
DistMatrix<F,MC, STAR> U01_MC_STAR(g);
DistMatrix<F,MR, STAR> U01_MR_STAR(g);
DistMatrix<F,VC, STAR> U01_VC_STAR(g);
DistMatrix<F,VC, STAR> Y01_VC_STAR(g);
DistMatrix<F> Y01(g);
PartitionDownDiagonal
( A, ATL, ATR,
ABL, ABR, 0 );
LockedPartitionDownDiagonal
( U, UTL, UTR,
UBL, UBR, 0 );
while( ATL.Height() < A.Height() )
{
RepartitionDownDiagonal
( ATL, /**/ ATR, A00, /**/ A01, A02,
/*************/ /******************/
/**/ A10, /**/ A11, A12,
ABL, /**/ ABR, A20, /**/ A21, A22 );
LockedRepartitionDownDiagonal
( UTL, /**/ UTR, U00, /**/ U01, U02,
/*************/ /******************/
/**/ U10, /**/ U11, U12,
UBL, /**/ UBR, U20, /**/ U21, U22 );
A01_MC_STAR.AlignWith( A00 );
A01_MR_STAR.AlignWith( A00 );
A01_VC_STAR.AlignWith( A00 );
U01_MC_STAR.AlignWith( A00 );
U01_MR_STAR.AlignWith( A00 );
U01_VC_STAR.AlignWith( A00 );
Y01.AlignWith( A01 );
Y01_VC_STAR.AlignWith( A01 );
//--------------------------------------------------------------------//
// Y01 := U01 A11
A11_STAR_STAR = A11;
U01_VC_STAR = U01;
Y01_VC_STAR.ResizeTo( A01.Height(), A01.Width() );
Hemm
( RIGHT, UPPER,
F(1), A11_STAR_STAR.LocalMatrix(), U01_VC_STAR.LocalMatrix(),
F(0), Y01_VC_STAR.LocalMatrix() );
Y01 = Y01_VC_STAR;
// A01 := U00 A01
Trmm( LEFT, UPPER, NORMAL, diag, F(1), U00, A01 );
// A01 := A01 + 1/2 Y01
Axpy( F(1)/F(2), Y01, A01 );
// A00 := A00 + (U01 A01' + A01 U01')
A01_MC_STAR = A01;
U01_MC_STAR = U01;
A01_VC_STAR = A01_MC_STAR;
A01_MR_STAR = A01_VC_STAR;
U01_MR_STAR = U01_MC_STAR;
LocalTrr2k
( UPPER, ADJOINT, ADJOINT,
F(1), U01_MC_STAR, A01_MR_STAR,
A01_MC_STAR, U01_MR_STAR,
F(1), A00 );
// A01 := A01 + 1/2 Y01
Axpy( F(1)/F(2), Y01_VC_STAR, A01_VC_STAR );
// A01 := A01 U11'
U11_STAR_STAR = U11;
LocalTrmm
//.........这里部分代码省略.........
开发者ID:jimgoo，项目名称:Elemental，代码行数:101，代码来源:UVar5.hpp
示例19: PushCallStack
inline void
SymmLLA
( T alpha, const DistMatrix<T>& A, const DistMatrix<T>& B,
T beta, DistMatrix<T>& C )
{
#ifndef RELEASE
PushCallStack("internal::SymmLLA");
if( A.Grid() != B.Grid() || B.Grid() != C.Grid() )
throw std::logic_error
("{A,B,C} must be distributed over the same grid");
#endif
const Grid& g = A.Grid();
DistMatrix<T>
BL(g), BR(g),
B0(g), B1(g), B2(g);
DistMatrix<T>
CL(g), CR(g),
C0(g), C1(g), C2(g);
DistMatrix<T,MC,STAR> B1_MC_STAR(g);
DistMatrix<T,VR,STAR> B1_VR_STAR(g);
DistMatrix<T,STAR,MR> B1Trans_STAR_MR(g);
DistMatrix<T> Z1(g);
DistMatrix<T,MC,STAR> Z1_MC_STAR(g);
DistMatrix<T,MR,STAR> Z1_MR_STAR(g);
DistMatrix<T,MR,MC > Z1_MR_MC(g);
B1_MC_STAR.AlignWith( A );
B1_VR_STAR.AlignWith( A );
B1Trans_STAR_MR.AlignWith( A );
Z1_MC_STAR.AlignWith( A );
Z1_MR_STAR.AlignWith( A );
Scale( beta, C );
LockedPartitionRight
( B, BL, BR, 0 );
PartitionRight
( C, CL, CR, 0 );
while( CL.Width() < C.Width() )
{
LockedRepartitionRight
( BL, /**/ BR,
B0, /**/ B1, B2 );
RepartitionRight
( CL, /**/ CR,
C0, /**/ C1, C2 );
Z1.AlignWith( C1 );
Zeros( C1.Height(), C1.Width(), Z1_MC_STAR );
Zeros( C1.Height(), C1.Width(), Z1_MR_STAR );
//--------------------------------------------------------------------//
B1_MC_STAR = B1;
B1_VR_STAR = B1_MC_STAR;
B1Trans_STAR_MR.TransposeFrom( B1_VR_STAR );
LocalSymmetricAccumulateLL
( TRANSPOSE,
alpha, A, B1_MC_STAR, B1Trans_STAR_MR, Z1_MC_STAR, Z1_MR_STAR );
Z1_MR_MC.SumScatterFrom( Z1_MR_STAR );
Z1 = Z1_MR_MC;
Z1.SumScatterUpdate( T(1), Z1_MC_STAR );
Axpy( T(1), Z1, C1 );
//--------------------------------------------------------------------//
Z1.FreeAlignments();
SlideLockedPartitionRight
( BL, /**/ BR,
B0, B1, /**/ B2 );
SlidePartitionRight
( CL, /**/ CR,
C0, C1, /**/ C2 );
}
#ifndef RELEASE
PopCallStack();
#endif
}
开发者ID:jimgoo，项目名称:Elemental，代码行数:80，代码来源:LL.hpp
示例20: entry
inline void
TwoSidedTrsmUVar1
( UnitOrNonUnit diag, DistMatrix<F>& A, const DistMatrix<F>& U )
{
#ifndef RELEASE
CallStackEntry entry("internal::TwoSidedTrsmUVar1");
if( A.Height() != A.Width() )
LogicError("A must be square");
if( U.Height() != U.Width() )
LogicError("Triangular matrices must be square");
if( A.Height() != U.Height() )
LogicError("A and U must be the same size");
#endif
const Grid& g = A.Grid();
// Matrix views
DistMatrix<F>
ATL(g), ATR(g), A00(g), A01(g), A02(g),
ABL(g), ABR(g), A10(g), A11(g), A12(g),
A20(g), A21(g), A22(g);
DistMatrix<F>
UTL(g), UTR(g), U00(g), U01(g), U02(g),
UBL(g), UBR(g), U10(g), U11(g), U12(g),
U20(g), U21(g), U22(g);
// Temporary distributions
DistMatrix<F,STAR,STAR> A11_STAR_STAR(g);
DistMatrix<F,VC, STAR> A01_VC_STAR(g);
DistMatrix<F,STAR,STAR> U11_STAR_STAR(g);
DistMatrix<F,MC, STAR> U01_MC_STAR(g);
DistMatrix<F,VC, STAR> U01_VC_STAR(g);
DistMatrix<F,VR, STAR> U01_VR_STAR(g);
DistMatrix<F,STAR,MR > U01Adj_STAR_MR(g);
DistMatrix<F,STAR,STAR> X11_STAR_STAR(g);
DistMatrix<F,MR, MC > Z01_MR_MC(g);
DistMatrix<F,MC, STAR> Z01_MC_STAR(g);
DistMatrix<F,MR, STAR> Z01_MR_STAR(g);
DistMatrix<F> Y01(g);
PartitionDownDiagonal
( A, ATL, ATR,
ABL, ABR, 0 );
LockedPartitionDownDiagonal
( U, UTL, UTR,
UBL, UBR, 0 );
while( ATL.Height() < A.Height() )
{
RepartitionDownDiagonal
( ATL, /**/ ATR, A00, /**/ A01, A02,
/*************/ /******************/
/**/ A10, /**/ A11, A12,
ABL, /**/ ABR, A20, /**/ A21, A22 );
LockedRepartitionDownDiagonal
( UTL, /**/ UTR, U00, /**/ U01, U02,
/*************/ /******************/
/**/ U10, /**/ U11, U12,
UBL, /**/ UBR, U20, /**/ U21, U22 );
A01_VC_STAR.AlignWith( A01 );
U01_MC_STAR.AlignWith( A00 );
U01_VR_STAR.AlignWith( A00 );
U01_VC_STAR.AlignWith( A00 );
U01Adj_STAR_MR.AlignWith( A00 );
Y01.AlignWith( A01 );
Z01_MR_MC.AlignWith( A01 );
Z01_MC_STAR.AlignWith( A00 );
Z01_MR_STAR.AlignWith( A00 );
//--------------------------------------------------------------------//
// Y01 := A00 U01
U01_MC_STAR = U01;
U01_VR_STAR = U01_MC_STAR;
U01Adj_STAR_MR.AdjointFrom( U01_VR_STAR );
Zeros( Z01_MC_STAR, A01.Height(), A01.Width() );
Zeros( Z01_MR_STAR, A01.Height(), A01.Width() );
LocalSymmetricAccumulateLU
( ADJOINT,
F(1), A00, U01_MC_STAR, U01Adj_STAR_MR, Z01_MC_STAR, Z01_MR_STAR );
Z01_MR_MC.SumScatterFrom( Z01_MR_STAR );
Y01 = Z01_MR_MC;
Y01.SumScatterUpdate( F(1), Z01_MC_STAR );
// A01 := inv(U00)' A01
//
// This is the bottleneck because A01 only has blocksize columns
Trsm( LEFT, UPPER, ADJOINT, diag, F(1), U00, A01 );
// A01 := A01 - 1/2 Y01
Axpy( F(-1)/F(2), Y01, A01 );
// A11 := A11 - (U01' A01 + A01' U01)
A01_VC_STAR = A01;
U01_VC_STAR = U01_MC_STAR;
Zeros( X11_STAR_STAR, A11.Height(), A11.Width() );
Her2k
( UPPER, ADJOINT,
F(-1), A01_VC_STAR.Matrix(), U01_VC_STAR.Matrix(),
F(0), X11_STAR_STAR.Matrix() );
A11.SumScatterUpdate( F(1), X11_STAR_STAR );
//.........这里部分代码省略.........
开发者ID:khalid-hasanov，项目名称:Elemental，代码行数:101，代码来源:UVar1.hpp
注：本文中的Axpy函数示例整理自Github/MSDocs等源码及文档管理平台，相关代码片段筛选自各路编程大神贡献的开源项目，源码版权归原作者所有，传播和使用请参考对应项目的License；未经允许，请勿转载。 |
客服电话
APP下载
官方微信
请发表评论