本文整理匯總了C++中ABL函數的典型用法代碼示例。如果您正苦於以下問題：C++ ABL函數的具體用法？C++ ABL怎麽用？C++ ABL使用的例子？那麽, 這裏精選的函數代碼示例或許可以為您提供幫助。
在下文中一共展示了ABL函數的15個代碼示例，這些例子默認根據受歡迎程度排序。您可以為喜歡或者感覺有用的代碼點讚，您的評價將有助於係統推薦出更棒的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 );
}
}示例2: LQ
inline void
LQ( DistMatrix<R,MC,MR>& A )
{
#ifndef RELEASE
PushCallStack("LQ");
#endif
if( IsComplex<R>::val )
throw std::logic_error("Called real routine with complex datatype");
const Grid& g = A.Grid();
// Matrix views
DistMatrix<R,MC,MR>
ATL(g), ATR(g), A00(g), A01(g), A02(g), ATopPan(g), ABottomPan(g),
ABL(g), ABR(g), A10(g), A11(g), A12(g),
A20(g), A21(g), A22(g);
PartitionDownLeftDiagonal
( A, ATL, ATR,
ABL, ABR, 0 );
while( ATL.Height() < A.Height() && ATL.Width() < A.Width() )
{
RepartitionDownDiagonal
( ATL, /**/ ATR, A00, /**/ A01, A02,
/*************/ /******************/
/**/ A10, /**/ A11, A12,
ABL, /**/ ABR, A20, /**/ A21, A22 );
ATopPan.View1x2( A11, A12 );
ABottomPan.View1x2( A21, A22 );
//--------------------------------------------------------------------//
internal::PanelLQ( ATopPan );
ApplyPackedReflectors
( RIGHT, UPPER, HORIZONTAL, FORWARD, 0, ATopPan, ABottomPan );
//--------------------------------------------------------------------//
SlidePartitionDownDiagonal
( ATL, /**/ ATR, A00, A01, /**/ A02,
/**/ A10, A11, /**/ A12,
/*************/ /******************/
ABL, /**/ ABR, A20, A21, /**/ A22 );
}
#ifndef RELEASE
PopCallStack();
#endif
}示例3: PushCallStack
inline void
HemmRUC
( T alpha, const DistMatrix<T>& A,
const DistMatrix<T>& B,
T beta, DistMatrix<T>& C )
{
#ifndef RELEASE
PushCallStack("internal::HemmRUC");
if( A.Grid() != B.Grid() || B.Grid() != C.Grid() )
throw std::logic_error("{A,B,C} must be distributed on the same grid");
#endif
const Grid& g = A.Grid();
// Matrix views
DistMatrix<T>
ATL(g), ATR(g), A00(g), A01(g), A02(g), AColPan(g),
ABL(g), ABR(g), A10(g), A11(g), A12(g), ARowPan(g),
A20(g), A21(g), A22(g);
DistMatrix<T> BL(g), BR(g),
B0(g), B1(g), B2(g);
DistMatrix<T> CL(g), CR(g),
C0(g), C1(g), C2(g),
CLeft(g), CRight(g);
// Temporary distributions
DistMatrix<T,MC,STAR> B1_MC_STAR(g);
DistMatrix<T,VR, STAR> AColPan_VR_STAR(g);
DistMatrix<T,STAR,MR > AColPanAdj_STAR_MR(g);
DistMatrix<T,MR, STAR> ARowPanAdj_MR_STAR(g);
B1_MC_STAR.AlignWith( C );
// Start the algorithm
Scale( beta, C );
LockedPartitionDownDiagonal
( A, ATL, ATR,
ABL, ABR, 0 );
LockedPartitionRight( B, BL, BR, 0 );
PartitionRight( C, CL, CR, 0 );
while( CR.Width() > 0 )
{
LockedRepartitionDownDiagonal
( ATL, /**/ ATR, A00, /**/ A01, A02,
/*************/ /******************/
/**/ A10, /**/ A11, A12,
ABL, /**/ ABR, A20, /**/ A21, A22 );
LockedRepartitionRight
( BL, /**/ BR,
B0, /**/ B1, B2 );
RepartitionRight
( CL, /**/ CR,
C0, /**/ C1, C2 );
ARowPan.LockedView1x2( A11, A12 );
AColPan.LockedView2x1
( A01,
A11 );
CLeft.View1x2( C0, C1 );
CRight.View1x2( C1, C2 );
AColPan_VR_STAR.AlignWith( CLeft );
AColPanAdj_STAR_MR.AlignWith( CLeft );
ARowPanAdj_MR_STAR.AlignWith( CRight );
//--------------------------------------------------------------------//
B1_MC_STAR = B1;
AColPan_VR_STAR = AColPan;
AColPanAdj_STAR_MR.AdjointFrom( AColPan_VR_STAR );
ARowPanAdj_MR_STAR.AdjointFrom( ARowPan );
MakeTrapezoidal( LEFT, LOWER, 0, ARowPanAdj_MR_STAR );
MakeTrapezoidal( RIGHT, LOWER, -1, AColPanAdj_STAR_MR );
LocalGemm
( NORMAL, ADJOINT,
alpha, B1_MC_STAR, ARowPanAdj_MR_STAR, T(1), CRight );
LocalGemm
( NORMAL, NORMAL,
alpha, B1_MC_STAR, AColPanAdj_STAR_MR, T(1), CLeft );
//--------------------------------------------------------------------//
AColPan_VR_STAR.FreeAlignments();
AColPanAdj_STAR_MR.FreeAlignments();
ARowPanAdj_MR_STAR.FreeAlignments();
SlideLockedPartitionDownDiagonal
( ATL, /**/ ATR, A00, A01, /**/ A02,
/**/ A10, A11, /**/ A12,
/*************/ /******************/
ABL, /**/ ABR, A20, A21, /**/ A22 );
SlideLockedPartitionRight
( BL, /**/ BR,
B0, B1, /**/ B2 );
SlidePartitionRight
( CL, /**/ CR,
C0, C1, /**/ C2 );
//.........這裏部分代碼省略.........
示例4: LSquare
void LSquare( DistMatrix<R>& A )
{
#ifndef RELEASE
CallStackEntry entry("hermitian_tridiag::LSquare");
if( A.Height() != A.Width() )
throw std::logic_error("A must be square");
if( A.Grid().Height() != A.Grid().Width() )
throw std::logic_error("The process grid must be square");
#endif
const Grid& g = A.Grid();
// Matrix views
DistMatrix<R>
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);
// Temporary distributions
DistMatrix<R> WPan(g);
DistMatrix<R,STAR,STAR> A11_STAR_STAR(g);
DistMatrix<R,MC, STAR> APan_MC_STAR(g), A11_MC_STAR(g),
A21_MC_STAR(g);
DistMatrix<R,MR, STAR> APan_MR_STAR(g), A11_MR_STAR(g),
A21_MR_STAR(g);
DistMatrix<R,MC, STAR> WPan_MC_STAR(g), W11_MC_STAR(g),
W21_MC_STAR(g);
DistMatrix<R,MR, STAR> WPan_MR_STAR(g), W11_MR_STAR(g),
W21_MR_STAR(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 );
if( A22.Height() > 0 )
{
WPan.AlignWith( A11 );
APan_MC_STAR.AlignWith( A11 );
WPan_MC_STAR.AlignWith( A11 );
APan_MR_STAR.AlignWith( A11 );
WPan_MR_STAR.AlignWith( A11 );
//----------------------------------------------------------------//
WPan.ResizeTo( ABR.Height(), A11.Width() );
APan_MC_STAR.ResizeTo( ABR.Height(), A11.Width() );
WPan_MC_STAR.ResizeTo( ABR.Height(), A11.Width() );
APan_MR_STAR.ResizeTo( ABR.Height(), A11.Width() );
WPan_MR_STAR.ResizeTo( ABR.Height(), A11.Width() );
hermitian_tridiag::PanelLSquare
( ABR, WPan,
APan_MC_STAR, APan_MR_STAR, WPan_MC_STAR, WPan_MR_STAR );
PartitionDown
( APan_MC_STAR, A11_MC_STAR,
A21_MC_STAR, A11.Height() );
PartitionDown
( APan_MR_STAR, A11_MR_STAR,
A21_MR_STAR, A11.Height() );
PartitionDown
( WPan_MC_STAR, W11_MC_STAR,
W21_MC_STAR, A11.Height() );
PartitionDown
( WPan_MR_STAR, W11_MR_STAR,
W21_MR_STAR, A11.Height() );
LocalTrr2k
( LOWER, TRANSPOSE, TRANSPOSE,
R(-1), A21_MC_STAR, W21_MR_STAR,
W21_MC_STAR, A21_MR_STAR,
R(1), A22 );
//----------------------------------------------------------------//
WPan_MR_STAR.FreeAlignments();
APan_MR_STAR.FreeAlignments();
WPan_MC_STAR.FreeAlignments();
APan_MC_STAR.FreeAlignments();
WPan.FreeAlignments();
}
else
{
A11_STAR_STAR = A11;
HermitianTridiag( LOWER, A11_STAR_STAR.Matrix() );
A11 = A11_STAR_STAR;
}
SlidePartitionDownDiagonal
( ATL, /**/ ATR, A00, A01, /**/ A02,
/**/ A10, A11, /**/ A12,
/*************/ /******************/
ABL, /**/ ABR, A20, A21, /**/ A22 );
}
}示例5: PushCallStack
inline void
SymmLLC
( T alpha, const DistMatrix<T>& A, const DistMatrix<T>& B,
T beta, DistMatrix<T>& C )
{
#ifndef RELEASE
PushCallStack("internal::SymmLLC");
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();
// Matrix views
DistMatrix<T>
ATL(g), ATR(g), A00(g), A01(g), A02(g), AColPan(g),
ABL(g), ABR(g), A10(g), A11(g), A12(g), ARowPan(g),
A20(g), A21(g), A22(g);
DistMatrix<T>
BT(g), B0(g),
BB(g), B1(g),
B2(g);
DistMatrix<T>
CT(g), C0(g), CAbove(g),
CB(g), C1(g), CBelow(g),
C2(g);
// Temporary distributions
DistMatrix<T,MC, STAR> AColPan_MC_STAR(g);
DistMatrix<T,STAR,MC > ARowPan_STAR_MC(g);
DistMatrix<T,MR, STAR> B1Trans_MR_STAR(g);
B1Trans_MR_STAR.AlignWith( C );
// Start the algorithm
Scale( beta, C );
LockedPartitionDownDiagonal
( A, ATL, ATR,
ABL, ABR, 0 );
LockedPartitionDown
( B, BT,
BB, 0 );
PartitionDown
( C, CT,
CB, 0 );
while( CB.Height() > 0 )
{
LockedRepartitionDownDiagonal
( ATL, /**/ ATR, A00, /**/ A01, A02,
/*************/ /******************/
/**/ A10, /**/ A11, A12,
ABL, /**/ ABR, A20, /**/ A21, A22 );
LockedRepartitionDown
( BT, B0,
/**/ /**/
B1,
BB, B2 );
RepartitionDown
( CT, C0,
/**/ /**/
C1,
CB, C2 );
LockedView1x2( ARowPan, A10, A11 );
LockedView2x1
( AColPan, A11,
A21 );
View2x1
( CAbove, C0,
C1 );
View2x1
( CBelow, C1,
C2 );
AColPan_MC_STAR.AlignWith( CBelow );
ARowPan_STAR_MC.AlignWith( CAbove );
//--------------------------------------------------------------------//
AColPan_MC_STAR = AColPan;
ARowPan_STAR_MC = ARowPan;
MakeTrapezoidal( LEFT, LOWER, 0, AColPan_MC_STAR );
MakeTrapezoidal( RIGHT, LOWER, -1, ARowPan_STAR_MC );
B1Trans_MR_STAR.TransposeFrom( B1 );
LocalGemm
( NORMAL, TRANSPOSE,
alpha, AColPan_MC_STAR, B1Trans_MR_STAR, T(1), CBelow );
LocalGemm
( TRANSPOSE, TRANSPOSE,
alpha, ARowPan_STAR_MC, B1Trans_MR_STAR, T(1), CAbove );
//--------------------------------------------------------------------//
AColPan_MC_STAR.FreeAlignments();
ARowPan_STAR_MC.FreeAlignments();
SlideLockedPartitionDownDiagonal
( ATL, /**/ ATR, A00, A01, /**/ A02,
//.........這裏部分代碼省略.........
示例6: 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 );
//.........這裏部分代碼省略.........
示例7: PanelHouseholder
inline void
PanelHouseholder( DistMatrix<F>& A, DistMatrix<F,MD,STAR>& t )
{
#ifndef RELEASE
CallStackEntry entry("lq::PanelHouseholder");
if( A.Grid() != t.Grid() )
LogicError("{A,t} must be distributed over the same grid");
if( t.Height() != Min(A.Height(),A.Width()) || t.Width() != 1 )
LogicError
("t must be a vector of height equal to the minimum dimension of A");
if( !t.AlignedWithDiagonal( A, 0 ) )
LogicError("t must be aligned with A's main diagonal");
#endif
const Grid& g = A.Grid();
// Matrix views
DistMatrix<F>
ATL(g), ATR(g), A00(g), a01(g), A02(g), aTopRow(g), ABottomPan(g),
ABL(g), ABR(g), a10(g), alpha11(g), a12(g),
A20(g), a21(g), A22(g);
DistMatrix<F,MD,STAR>
tT(g), t0(g),
tB(g), tau1(g),
t2(g);
// Temporary distributions
DistMatrix<F> aTopRowConj(g);
DistMatrix<F,STAR,MR > aTopRowConj_STAR_MR(g);
DistMatrix<F,MC, STAR> z_MC_STAR(g);
PartitionDownDiagonal
( A, ATL, ATR,
ABL, ABR, 0 );
PartitionDown
( t, tT,
tB, 0 );
while( ATL.Height() < A.Height() && ATL.Width() < A.Width() )
{
RepartitionDownDiagonal
( ATL, /**/ ATR, A00, /**/ a01, A02,
/*************/ /**********************/
/**/ a10, /**/ alpha11, a12,
ABL, /**/ ABR, A20, /**/ a21, A22, 1 );
RepartitionDown
( tT, t0,
/**/ /****/
tau1,
tB, t2, 1 );
View1x2( aTopRow, alpha11, a12 );
View1x2( ABottomPan, a21, A22 );
aTopRowConj_STAR_MR.AlignWith( ABottomPan );
z_MC_STAR.AlignWith( ABottomPan );
//--------------------------------------------------------------------//
// Compute the Householder reflector
const F tau = Reflector( alpha11, a12 );
tau1.Set( 0, 0, tau );
// Apply the Householder reflector
const bool myDiagonalEntry = ( g.Row() == alpha11.ColAlignment() &&
g.Col() == alpha11.RowAlignment() );
F alpha = 0;
if( myDiagonalEntry )
{
alpha = alpha11.GetLocal(0,0);
alpha11.SetLocal(0,0,1);
}
Conjugate( aTopRow, aTopRowConj );
aTopRowConj_STAR_MR = aTopRowConj;
Zeros( z_MC_STAR, ABottomPan.Height(), 1 );
LocalGemv
( NORMAL, F(1), ABottomPan, aTopRowConj_STAR_MR, F(0), z_MC_STAR );
z_MC_STAR.SumOverRow();
Ger
( -Conj(tau),
z_MC_STAR.LockedMatrix(),
aTopRowConj_STAR_MR.LockedMatrix(),
ABottomPan.Matrix() );
if( myDiagonalEntry )
alpha11.SetLocal(0,0,alpha);
//--------------------------------------------------------------------//
SlidePartitionDown
( tT, t0,
tau1,
/**/ /****/
tB, t2 );
SlidePartitionDownDiagonal
( ATL, /**/ ATR, A00, a01, /**/ A02,
/**/ a10, alpha11, /**/ a12,
/*************/ /**********************/
ABL, /**/ ABR, A20, a21, /**/ A22 );
}
}示例8: Householder
inline void
Householder( DistMatrix<F>& A, DistMatrix<F,MD,STAR>& t )
{
#ifndef RELEASE
CallStackEntry entry("qr::Householder");
if( A.Grid() != t.Grid() )
LogicError("{A,s} must be distributed over the same grid");
#endif
const Grid& g = A.Grid();
if( t.Viewing() )
{
if( !t.AlignedWithDiagonal( A ) )
LogicError("t was not aligned with A");
}
else
{
t.AlignWithDiagonal( A );
}
t.ResizeTo( Min(A.Height(),A.Width()), 1 );
// Matrix views
DistMatrix<F>
ATL(g), ATR(g), A00(g), A01(g), A02(g), ALeftPan(g), ARightPan(g),
ABL(g), ABR(g), A10(g), A11(g), A12(g),
A20(g), A21(g), A22(g);
DistMatrix<F,MD,STAR>
tT(g), t0(g),
tB(g), t1(g),
t2(g);
PartitionDownDiagonal
( A, ATL, ATR,
ABL, ABR, 0 );
PartitionDown
( t, tT,
tB, 0 );
while( ATL.Height() < A.Height() && ATL.Width() < A.Width() )
{
RepartitionDownDiagonal
( ATL, /**/ ATR, A00, /**/ A01, A02,
/*************/ /******************/
/**/ A10, /**/ A11, A12,
ABL, /**/ ABR, A20, /**/ A21, A22 );
RepartitionDown
( tT, t0,
/**/ /**/
t1,
tB, t2 );
View2x1
( ALeftPan, A11,
A21 );
View2x1
( ARightPan, A12,
A22 );
//--------------------------------------------------------------------//
PanelHouseholder( ALeftPan, t1 );
ApplyQ( LEFT, ADJOINT, ALeftPan, t1, ARightPan );
//--------------------------------------------------------------------//
SlidePartitionDown
( tT, t0,
t1,
/**/ /**/
tB, t2 );
SlidePartitionDownDiagonal
( ATL, /**/ ATR, A00, A01, /**/ A02,
/**/ A10, A11, /**/ A12,
/*************/ /******************/
ABL, /**/ ABR, A20, A21, /**/ A22 );
}
}示例9: PushCallStack
inline void
internal::HermitianTridiagU
( DistMatrix<Complex<R>,MC, MR >& A,
DistMatrix<Complex<R>,STAR,STAR>& t )
{
#ifndef RELEASE
PushCallStack("internal::HermitianTridiagU");
if( A.Grid() != t.Grid() )
throw std::logic_error("{A,t} must be distributed over the same grid");
if( A.Height() != A.Width() )
throw std::logic_error("A must be square");
if( t.Viewing() )
throw std::logic_error("t must not be a view");
#endif
typedef Complex<R> C;
const Grid& g = A.Grid();
DistMatrix<C,MD,STAR> tDiag(g);
tDiag.AlignWithDiagonal( A, 1 );
tDiag.ResizeTo( A.Height()-1, 1 );
if( g.InGrid() )
{
// Matrix views
DistMatrix<C,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<C,MD,STAR> tT(g), t0(g),
tB(g), t1(g),
t2(g);
// Temporary distributions
DistMatrix<C,MC, MR > WPan(g);
DistMatrix<C,STAR,STAR> t1_STAR_STAR(g);
DistMatrix<C,STAR,STAR> A11_STAR_STAR(g);
DistMatrix<C,MC, STAR> APan_MC_STAR(g), A01_MC_STAR(g),
A11_MC_STAR(g);
DistMatrix<C,MR, STAR> APan_MR_STAR(g), A01_MR_STAR(g),
A11_MR_STAR(g);
DistMatrix<C,MC, STAR> WPan_MC_STAR(g), W01_MC_STAR(g),
W11_MC_STAR(g);
DistMatrix<C,MR, STAR> WPan_MR_STAR(g), W01_MR_STAR(g),
W11_MR_STAR(g);
PartitionUpDiagonal
( A, ATL, ATR,
ABL, ABR, 0 );
PartitionUp
( tDiag, tT,
tB, 0 );
while( ABR.Height() < A.Height() )
{
RepartitionUpDiagonal
( ATL, /**/ ATR, A00, A01, /**/ A02,
/**/ A10, A11, /**/ A12,
/*************/ /******************/
ABL, /**/ ABR, A20, A21, /**/ A22 );
RepartitionUp
( tT, t0,
t1,
/**/ /**/
tB, t2 );
if( A00.Height() > 0 )
{
WPan.AlignWith( A01 );
APan_MC_STAR.AlignWith( A00 );
WPan_MC_STAR.AlignWith( A00 );
APan_MR_STAR.AlignWith( A00 );
WPan_MR_STAR.AlignWith( A00 );
//------------------------------------------------------------//
WPan.ResizeTo( ATL.Height(), A11.Width() );
APan_MC_STAR.ResizeTo( ATL.Height(), A11.Width() );
WPan_MC_STAR.ResizeTo( ATL.Height(), A11.Width() );
APan_MR_STAR.ResizeTo( ATL.Height(), A11.Width() );
WPan_MR_STAR.ResizeTo( ATL.Height(), A11.Width() );
internal::HermitianPanelTridiagU
( ATL, WPan, t1,
APan_MC_STAR, APan_MR_STAR, WPan_MC_STAR, WPan_MR_STAR );
PartitionUp
( APan_MC_STAR, A01_MC_STAR,
A11_MC_STAR, A11.Height() );
PartitionUp
( APan_MR_STAR, A01_MR_STAR,
A11_MR_STAR, A11.Height() );
PartitionUp
( WPan_MC_STAR, W01_MC_STAR,
W11_MC_STAR, A11.Height() );
PartitionUp
( WPan_MR_STAR, W01_MR_STAR,
W11_MR_STAR, A11.Height() );
internal::LocalTrr2k
( UPPER, ADJOINT, ADJOINT,
(C)-1, A01_MC_STAR, W01_MR_STAR,
W01_MC_STAR, A01_MR_STAR,
//.........這裏部分代碼省略.........
示例10: PushCallStack
inline void UnblockedBidiagU
( DistMatrix<Complex<R> >& A,
DistMatrix<Complex<R>,MD,STAR>& tP,
DistMatrix<Complex<R>,MD,STAR>& tQ )
{
#ifndef RELEASE
PushCallStack("BidiagU");
#endif
const int tPHeight = std::max(A.Width()-1,0);
const int tQHeight = A.Width();
#ifndef RELEASE
if( A.Grid() != tP.Grid() || tP.Grid() != tQ.Grid() )
throw std::logic_error("Process grids do not match");
if( A.Height() < A.Width() )
throw std::logic_error("A must be at least as tall as it is wide");
if( tP.Viewing() && (tP.Height() != tPHeight || tP.Width() != 1) )
throw std::logic_error("tP is the wrong height");
if( tQ.Viewing() && (tQ.Height() != tQHeight || tQ.Width() != 1) )
throw std::logic_error("tQ is the wrong height");
#endif
typedef Complex<R> C;
const Grid& g = A.Grid();
if( !tP.Viewing() )
tP.ResizeTo( tPHeight, 1 );
if( !tQ.Viewing() )
tQ.ResizeTo( tQHeight, 1 );
// Matrix views
DistMatrix<C>
ATL(g), ATR(g), A00(g), a01(g), A02(g), alpha12L(g), a12R(g),
ABL(g), ABR(g), a10(g), alpha11(g), a12(g), aB1(g), AB2(g),
A20(g), a21(g), A22(g);
// Temporary matrices
DistMatrix<C,STAR,MR > a12_STAR_MR(g);
DistMatrix<C,MC, STAR> aB1_MC_STAR(g);
DistMatrix<C,MR, STAR> x12Adj_MR_STAR(g);
DistMatrix<C,MC, STAR> w21_MC_STAR(g);
PushBlocksizeStack( 1 );
PartitionDownDiagonal
( A, ATL, ATR,
ABL, ABR, 0 );
while( ATL.Width() < A.Width() )
{
RepartitionDownDiagonal
( ATL, /**/ ATR, A00, /**/ a01, A02,
/*************/ /**********************/
/**/ a10, /**/ alpha11, a12,
ABL, /**/ ABR, A20, /**/ a21, A22 );
View2x1
( aB1, alpha11,
a21 );
View2x1
( AB2, a12,
A22 );
aB1_MC_STAR.AlignWith( aB1 );
a12_STAR_MR.AlignWith( a12 );
x12Adj_MR_STAR.AlignWith( AB2 );
w21_MC_STAR.AlignWith( A22 );
Zeros( a12.Width(), 1, x12Adj_MR_STAR );
Zeros( a21.Height(), 1, w21_MC_STAR );
const bool thisIsMyRow = ( g.Row() == alpha11.ColAlignment() );
const bool thisIsMyCol = ( g.Col() == alpha11.RowAlignment() );
const bool nextIsMyCol = ( g.Col() == a12.RowAlignment() );
//--------------------------------------------------------------------//
// Find tauQ, u, and epsilonQ such that
// I - conj(tauQ) | 1 | | 1, u^H | | alpha11 | = | epsilonQ |
// | u | | a21 | | 0 |
const C tauQ = Reflector( alpha11, a21 );
tQ.Set(A00.Height(),0,tauQ );
C epsilonQ=0;
if( thisIsMyCol && thisIsMyRow )
epsilonQ = alpha11.GetLocal(0,0);
// Set aB1 = | 1 | and form x12^H := (aB1^H AB2)^H = AB2^H aB1
// | u |
alpha11.Set(0,0,C(1));
aB1_MC_STAR = aB1;
internal::LocalGemv
( ADJOINT, C(1), AB2, aB1_MC_STAR, C(0), x12Adj_MR_STAR );
x12Adj_MR_STAR.SumOverCol();
// Update AB2 := AB2 - conj(tauQ) aB1 x12
// = AB2 - conj(tauQ) aB1 aB1^H AB2
// = (I - conj(tauQ) aB1 aB1^H) AB2
internal::LocalGer( -Conj(tauQ), aB1_MC_STAR, x12Adj_MR_STAR, AB2 );
// Put epsilonQ back instead of the temporary value, 1
if( thisIsMyCol && thisIsMyRow )
alpha11.SetLocal(0,0,epsilonQ);
if( A22.Width() != 0 )
{
// Due to the deficiencies in the BLAS ?gemv routines, this section
// is easier if we temporarily conjugate a12
//.........這裏部分代碼省略.........
示例11: UnblockedBidiagU
inline void UnblockedBidiagU( DistMatrix<R>& A )
{
#ifndef RELEASE
PushCallStack("bidiag::UnblockedBidiagU");
if( A.Height() < A.Width() )
throw std::logic_error("A must be at least as tall as it is wide");
#endif
const Grid& g = A.Grid();
// Matrix views
DistMatrix<R>
ATL(g), ATR(g), A00(g), a01(g), A02(g), alpha12L(g), a12R(g),
ABL(g), ABR(g), a10(g), alpha11(g), a12(g), aB1(g), AB2(g),
A20(g), a21(g), A22(g);
// Temporary matrices
DistMatrix<R,STAR,MR > a12_STAR_MR(g);
DistMatrix<R,MC, STAR> aB1_MC_STAR(g);
DistMatrix<R,MR, STAR> x12Trans_MR_STAR(g);
DistMatrix<R,MC, STAR> w21_MC_STAR(g);
PushBlocksizeStack( 1 );
PartitionDownDiagonal
( A, ATL, ATR,
ABL, ABR, 0 );
while( ATL.Width() < A.Width() )
{
RepartitionDownDiagonal
( ATL, /**/ ATR, A00, /**/ a01, A02,
/*************/ /**********************/
/**/ a10, /**/ alpha11, a12,
ABL, /**/ ABR, A20, /**/ a21, A22 );
View2x1
( aB1, alpha11,
a21 );
View2x1
( AB2, a12,
A22 );
aB1_MC_STAR.AlignWith( aB1 );
a12_STAR_MR.AlignWith( a12 );
x12Trans_MR_STAR.AlignWith( AB2 );
w21_MC_STAR.AlignWith( A22 );
Zeros( a12.Width(), 1, x12Trans_MR_STAR );
Zeros( a21.Height(), 1, w21_MC_STAR );
const bool thisIsMyRow = ( g.Row() == alpha11.ColAlignment() );
const bool thisIsMyCol = ( g.Col() == alpha11.RowAlignment() );
const bool nextIsMyCol = ( g.Col() == a12.RowAlignment() );
//--------------------------------------------------------------------//
// Find tauQ, u, and epsilonQ such that
// I - tauQ | 1 | | 1, u^T | | alpha11 | = | epsilonQ |
// | u | | a21 | = | 0 |
const R tauQ = Reflector( alpha11, a21 );
R epsilonQ=0;
if( thisIsMyCol && thisIsMyRow )
epsilonQ = alpha11.GetLocal(0,0);
// Set aB1 = | 1 | and form x12^T := (aB1^T AB2)^T = AB2^T aB1
// | u |
alpha11.Set(0,0,R(1));
aB1_MC_STAR = aB1;
internal::LocalGemv
( TRANSPOSE, R(1), AB2, aB1_MC_STAR, R(0), x12Trans_MR_STAR );
x12Trans_MR_STAR.SumOverCol();
// Update AB2 := AB2 - tauQ aB1 x12
// = AB2 - tauQ aB1 aB1^T AB2
// = (I - tauQ aB1 aB1^T) AB2
internal::LocalGer( -tauQ, aB1_MC_STAR, x12Trans_MR_STAR, AB2 );
// Put epsilonQ back instead of the temporary value, 1
if( thisIsMyCol && thisIsMyRow )
alpha11.SetLocal(0,0,epsilonQ);
if( A22.Width() != 0 )
{
// Expose the subvector we seek to zero, a12R
PartitionRight( a12, alpha12L, a12R );
// Find tauP, v, and epsilonP such that
// I - tauP | 1 | | 1, v^T | | alpha12L | = | epsilonP |
// | v | | a12R^T | = | 0 |
const R tauP = Reflector( alpha12L, a12R );
R epsilonP=0;
if( nextIsMyCol && thisIsMyRow )
epsilonP = alpha12L.GetLocal(0,0);
// Set a12^T = | 1 | and form w21 := A22 a12^T = A22 | 1 |
// | v | | v |
alpha12L.Set(0,0,R(1));
a12_STAR_MR = a12;
internal::LocalGemv
( NORMAL, R(1), A22, a12_STAR_MR, R(0), w21_MC_STAR );
w21_MC_STAR.SumOverRow();
// A22 := A22 - tauP w21 a12
// = A22 - tauP A22 a12^T a12
// = A22 (I - tauP a12^T a12)
//.........這裏部分代碼省略.........
示例12: PushCallStack
inline void
PanelLU
( DistMatrix<F, STAR,STAR>& A,
DistMatrix<F, MC, STAR>& B,
DistMatrix<int,STAR,STAR>& p,
int pivotOffset )
{
#ifndef RELEASE
PushCallStack("internal::PanelLU");
if( A.Grid() != p.Grid() || p.Grid() != B.Grid() )
throw std::logic_error
("Matrices must be distributed over the same grid");
if( A.Width() != B.Width() )
throw std::logic_error("A and B must be the same width");
if( A.Height() != p.Height() || p.Width() != 1 )
throw std::logic_error("p must be a vector that conforms with A");
#endif
const Grid& g = A.Grid();
const int r = g.Height();
const int colShift = B.ColShift();
const int colAlignment = B.ColAlignment();
// Matrix views
DistMatrix<F,STAR,STAR>
ATL(g), ATR(g), A00(g), a01(g), A02(g),
ABL(g), ABR(g), a10(g), alpha11(g), a12(g),
A20(g), a21(g), A22(g);
DistMatrix<F,MC,STAR>
BL(g), BR(g),
B0(g), b1(g), B2(g);
const int width = A.Width();
const int numBytes = (width+1)*sizeof(F)+sizeof(int);
std::vector<byte> sendData(numBytes);
std::vector<byte> recvData(numBytes);
// Extract pointers to send and recv data
// TODO: Think of how to make this safer with respect to alignment issues
F* sendBufFloat = (F*)&sendData[0];
F* recvBufFloat = (F*)&recvData[0];
int* sendBufInt = (int*)&sendData[(width+1)*sizeof(F)];
int* recvBufInt = (int*)&recvData[(width+1)*sizeof(F)];
// Start the algorithm
PushBlocksizeStack( 1 );
PartitionDownDiagonal
( A, ATL, ATR,
ABL, ABR, 0 );
PartitionRight( B, BL, BR, 0 );
while( ATL.Height() < A.Height() )
{
RepartitionDownDiagonal
( ATL, /**/ ATR, A00, /**/ a01, A02,
/*************/ /**********************/
/**/ a10, /**/ alpha11, a12,
ABL, /**/ ABR, A20, /**/ a21, A22 );
RepartitionRight
( BL, /**/ BR,
B0, /**/ b1, B2 );
//--------------------------------------------------------------------//
const int currentRow = a01.Height();
// Store the index/value of the pivot candidate in A
F pivot = alpha11.GetLocal(0,0);
int pivotRow = currentRow;
for( int i=0; i<a21.Height(); ++i )
{
F value = a21.GetLocal(i,0);
if( FastAbs(value) > FastAbs(pivot) )
{
pivot = value;
pivotRow = currentRow + i + 1;
}
}
// Update the pivot candidate to include local data from B
for( int i=0; i<B.LocalHeight(); ++i )
{
F value = b1.GetLocal(i,0);
if( FastAbs(value) > FastAbs(pivot) )
{
pivot = value;
pivotRow = A.Height() + colShift + i*r;
}
}
// Fill the send buffer with:
// [ pivotValue | pivot row data | pivotRow ]
if( pivotRow < A.Height() )
{
sendBufFloat[0] = A.GetLocal(pivotRow,a10.Width());
const int ALDim = A.LocalLDim();
const F* ABuffer = A.LocalBuffer(pivotRow,0);
for( int j=0; j<width; ++j )
sendBufFloat[j+1] = ABuffer[j*ALDim];
}
//.........這裏部分代碼省略.........
示例13: PushCallStack
inline void
TwoSidedTrsmLVar2
( UnitOrNonUnit diag, DistMatrix<F>& A, const DistMatrix<F>& L )
{
#ifndef RELEASE
PushCallStack("internal::TwoSidedTrsmLVar2");
if( A.Height() != A.Width() )
throw std::logic_error("A must be square");
if( L.Height() != L.Width() )
throw std::logic_error("Triangular matrices must be square");
if( A.Height() != L.Height() )
throw std::logic_error("A and L 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>
LTL(g), LTR(g), L00(g), L01(g), L02(g),
LBL(g), LBR(g), L10(g), L11(g), L12(g),
L20(g), L21(g), L22(g);
// Temporary distributions
DistMatrix<F,MR, STAR> A10Adj_MR_STAR(g);
DistMatrix<F,STAR,VR > A10_STAR_VR(g);
DistMatrix<F,STAR,STAR> A11_STAR_STAR(g);
DistMatrix<F,VC, STAR> A21_VC_STAR(g);
DistMatrix<F,MR, STAR> F10Adj_MR_STAR(g);
DistMatrix<F,MR, STAR> L10Adj_MR_STAR(g);
DistMatrix<F,VC, STAR> L10Adj_VC_STAR(g);
DistMatrix<F,STAR,MC > L10_STAR_MC(g);
DistMatrix<F,STAR,STAR> L11_STAR_STAR(g);
DistMatrix<F,MC, STAR> X11_MC_STAR(g);
DistMatrix<F,MC, STAR> X21_MC_STAR(g);
DistMatrix<F,MC, STAR> Y10Adj_MC_STAR(g);
DistMatrix<F,MR, MC > Y10Adj_MR_MC(g);
DistMatrix<F> X11(g);
DistMatrix<F> Y10Adj(g);
Matrix<F> Y10Local;
PartitionDownDiagonal
( A, ATL, ATR,
ABL, ABR, 0 );
LockedPartitionDownDiagonal
( L, LTL, LTR,
LBL, LBR, 0 );
while( ATL.Height() < A.Height() )
{
RepartitionDownDiagonal
( ATL, /**/ ATR, A00, /**/ A01, A02,
/*************/ /******************/
/**/ A10, /**/ A11, A12,
ABL, /**/ ABR, A20, /**/ A21, A22 );
LockedRepartitionDownDiagonal
( LTL, /**/ LTR, L00, /**/ L01, L02,
/*************/ /******************/
/**/ L10, /**/ L11, L12,
LBL, /**/ LBR, L20, /**/ L21, L22 );
A10Adj_MR_STAR.AlignWith( L10 );
F10Adj_MR_STAR.AlignWith( A00 );
L10Adj_MR_STAR.AlignWith( A00 );
L10Adj_VC_STAR.AlignWith( A00 );
L10_STAR_MC.AlignWith( A00 );
X11.AlignWith( A11 );
X11_MC_STAR.AlignWith( L10 );
X21_MC_STAR.AlignWith( A20 );
Y10Adj_MC_STAR.AlignWith( A00 );
Y10Adj_MR_MC.AlignWith( A10 );
//--------------------------------------------------------------------//
// Y10 := L10 A00
L10Adj_MR_STAR.AdjointFrom( L10 );
L10Adj_VC_STAR = L10Adj_MR_STAR;
L10_STAR_MC.AdjointFrom( L10Adj_VC_STAR );
Y10Adj_MC_STAR.ResizeTo( A10.Width(), A10.Height() );
F10Adj_MR_STAR.ResizeTo( A10.Width(), A10.Height() );
Zero( Y10Adj_MC_STAR );
Zero( F10Adj_MR_STAR );
LocalSymmetricAccumulateRL
( ADJOINT,
F(1), A00, L10_STAR_MC, L10Adj_MR_STAR,
Y10Adj_MC_STAR, F10Adj_MR_STAR );
Y10Adj.SumScatterFrom( Y10Adj_MC_STAR );
Y10Adj_MR_MC = Y10Adj;
Y10Adj_MR_MC.SumScatterUpdate( F(1), F10Adj_MR_STAR );
Adjoint( Y10Adj_MR_MC.LockedLocalMatrix(), Y10Local );
// X11 := A10 L10'
X11_MC_STAR.ResizeTo( A11.Height(), A11.Width() );
LocalGemm
( NORMAL, NORMAL, F(1), A10, L10Adj_MR_STAR, F(0), X11_MC_STAR );
// A10 := A10 - Y10
Axpy( F(-1), Y10Local, A10.LocalMatrix() );
//.........這裏部分代碼省略.........
示例14: entry
inline void
TwoSidedTrsmLVar5
( UnitOrNonUnit diag, DistMatrix<F>& A, const DistMatrix<F>& L )
{
#ifndef RELEASE
CallStackEntry entry("internal::TwoSidedTrsmLVar5");
if( A.Height() != A.Width() )
LogicError("A must be square");
if( L.Height() != L.Width() )
LogicError("Triangular matrices must be square");
if( A.Height() != L.Height() )
LogicError("A and L 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>
LTL(g), LTR(g), L00(g), L01(g), L02(g),
LBL(g), LBR(g), L10(g), L11(g), L12(g),
L20(g), L21(g), L22(g);
// Temporary distributions
DistMatrix<F,STAR,STAR> A11_STAR_STAR(g);
DistMatrix<F,MC, STAR> A21_MC_STAR(g);
DistMatrix<F,VC, STAR> A21_VC_STAR(g);
DistMatrix<F,VR, STAR> A21_VR_STAR(g);
DistMatrix<F,STAR,MR > A21Adj_STAR_MR(g);
DistMatrix<F,STAR,STAR> L11_STAR_STAR(g);
DistMatrix<F,MC, STAR> L21_MC_STAR(g);
DistMatrix<F,VC, STAR> L21_VC_STAR(g);
DistMatrix<F,VR, STAR> L21_VR_STAR(g);
DistMatrix<F,STAR,MR > L21Adj_STAR_MR(g);
DistMatrix<F,VC, STAR> Y21_VC_STAR(g);
DistMatrix<F> Y21(g);
PartitionDownDiagonal
( A, ATL, ATR,
ABL, ABR, 0 );
LockedPartitionDownDiagonal
( L, LTL, LTR,
LBL, LBR, 0 );
while( ATL.Height() < A.Height() )
{
RepartitionDownDiagonal
( ATL, /**/ ATR, A00, /**/ A01, A02,
/*************/ /******************/
/**/ A10, /**/ A11, A12,
ABL, /**/ ABR, A20, /**/ A21, A22 );
LockedRepartitionDownDiagonal
( LTL, /**/ LTR, L00, /**/ L01, L02,
/*************/ /******************/
/**/ L10, /**/ L11, L12,
LBL, /**/ LBR, L20, /**/ L21, L22 );
A21_MC_STAR.AlignWith( A22 );
A21_VC_STAR.AlignWith( A22 );
A21_VR_STAR.AlignWith( A22 );
A21Adj_STAR_MR.AlignWith( A22 );
L21_MC_STAR.AlignWith( A22 );
L21_VC_STAR.AlignWith( A22 );
L21_VR_STAR.AlignWith( A22 );
L21Adj_STAR_MR.AlignWith( A22 );
Y21.AlignWith( A21 );
Y21_VC_STAR.AlignWith( A22 );
//--------------------------------------------------------------------//
// A11 := inv(L11) A11 inv(L11)'
L11_STAR_STAR = L11;
A11_STAR_STAR = A11;
LocalTwoSidedTrsm( LOWER, diag, A11_STAR_STAR, L11_STAR_STAR );
A11 = A11_STAR_STAR;
// Y21 := L21 A11
L21_VC_STAR = L21;
Zeros( Y21_VC_STAR, A21.Height(), A21.Width() );
Hemm
( RIGHT, LOWER,
F(1), A11_STAR_STAR.Matrix(), L21_VC_STAR.Matrix(),
F(0), Y21_VC_STAR.Matrix() );
Y21 = Y21_VC_STAR;
// A21 := A21 inv(L11)'
A21_VC_STAR = A21;
LocalTrsm
( RIGHT, LOWER, ADJOINT, diag, F(1), L11_STAR_STAR, A21_VC_STAR );
A21 = A21_VC_STAR;
// A21 := A21 - 1/2 Y21
Axpy( F(-1)/F(2), Y21, A21 );
// A22 := A22 - (L21 A21' + A21 L21')
A21_MC_STAR = A21;
L21_MC_STAR = L21;
A21_VC_STAR = A21_MC_STAR;
A21_VR_STAR = A21_VC_STAR;
L21_VR_STAR = L21_VC_STAR;
//.........這裏部分代碼省略.........
示例15: PushCallStack
inline void
internal::HegstRLVar3( DistMatrix<F,MC,MR>& A, const DistMatrix<F,MC,MR>& L )
{
#ifndef RELEASE
PushCallStack("internal::HegstRLVar4");
if( A.Height() != A.Width() )
throw std::logic_error("A must be square");
if( L.Height() != L.Width() )
throw std::logic_error("Triangular matrices must be square");
if( A.Height() != L.Height() )
throw std::logic_error("A and L must be the same size");
#endif
const Grid& g = A.Grid();
// Matrix views
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>
YTL(g), YTR(g), Y00(g), Y01(g), Y02(g),
YBL(g), YBR(g), Y10(g), Y11(g), Y12(g),
Y20(g), Y21(g), Y22(g);
DistMatrix<F,MC,MR>
LTL(g), LTR(g), L00(g), L01(g), L02(g),
LBL(g), LBR(g), L10(g), L11(g), L12(g),
L20(g), L21(g), L22(g);
// Temporary distributions
DistMatrix<F,STAR,MR > A11_STAR_MR(g);
DistMatrix<F,STAR,STAR> A11_STAR_STAR(g);
DistMatrix<F,VC, STAR> A21_VC_STAR(g);
DistMatrix<F,STAR,VR > A10_STAR_VR(g);
DistMatrix<F,STAR,MR > A10_STAR_MR(g);
DistMatrix<F,STAR,STAR> L11_STAR_STAR(g);
DistMatrix<F,STAR,VR > L10_STAR_VR(g);
DistMatrix<F,STAR,MR > L10_STAR_MR(g);
DistMatrix<F,MC, STAR> L21_MC_STAR(g);
DistMatrix<F,STAR,STAR> X11_STAR_STAR(g);
DistMatrix<F,MC, STAR> X21_MC_STAR(g);
DistMatrix<F,MC, STAR> Z21_MC_STAR(g);
// We will use an entire extra matrix as temporary storage. If this is not
// acceptable, use HegstRLVar4 instead.
DistMatrix<F,MC,MR> Y(g);
Y.AlignWith( A );
Y.ResizeTo( A.Height(), A.Width() );
Zero( Y );
PartitionDownDiagonal
( A, ATL, ATR,
ABL, ABR, 0 );
PartitionDownDiagonal
( Y, YTL, YTR,
YBL, YBR, 0 );
LockedPartitionDownDiagonal
( L, LTL, LTR,
LBL, LBR, 0 );
while( ATL.Height() < A.Height() )
{
RepartitionDownDiagonal
( ATL, /**/ ATR, A00, /**/ A01, A02,
/*************/ /******************/
/**/ A10, /**/ A11, A12,
ABL, /**/ ABR, A20, /**/ A21, A22 );
RepartitionDownDiagonal
( YTL, /**/ YTR, Y00, /**/ Y01, Y02,
/*************/ /******************/
/**/ Y10, /**/ Y11, Y12,
YBL, /**/ YBR, Y20, /**/ Y21, Y22 );
LockedRepartitionDownDiagonal
( LTL, /**/ LTR, L00, /**/ L01, L02,
/*************/ /******************/
/**/ L10, /**/ L11, L12,
LBL, /**/ LBR, L20, /**/ L21, L22 );
A11_STAR_MR.AlignWith( Y21 );
A21_VC_STAR.AlignWith( A21 );
A10_STAR_VR.AlignWith( A10 );
A10_STAR_MR.AlignWith( A10 );
L10_STAR_VR.AlignWith( A10 );
L10_STAR_MR.AlignWith( A10 );
L21_MC_STAR.AlignWith( Y21 );
X21_MC_STAR.AlignWith( A20 );
Z21_MC_STAR.AlignWith( L20 );
//--------------------------------------------------------------------//
// A10 := A10 - 1/2 Y10
Axpy( (F)-0.5, Y10, A10 );
// A11 := A11 - (A10 L10' + L10 A10')
A10_STAR_VR = A10;
L10_STAR_VR = L10;
X11_STAR_STAR.ResizeTo( A11.Height(), A11.Width() );
Her2k
( LOWER, NORMAL,
(F)1, A10_STAR_VR.LocalMatrix(), L10_STAR_VR.LocalMatrix(),
(F)0, X11_STAR_STAR.LocalMatrix() );
MakeTrapezoidal( LEFT, LOWER, 0, X11_STAR_STAR );
//.........這裏部分代碼省略.........