本文整理匯總了C++中A22函數的典型用法代碼示例。如果您正苦於以下問題：C++ A22函數的具體用法？C++ A22怎麽用？C++ A22使用的例子？那麽, 這裏精選的函數代碼示例或許可以為您提供幫助。
在下文中一共展示了A22函數的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: main
int main()
{
Rat a1(-1,2);
Rat a2(2,9);
Rat b2(3,2);
Rat c1(-1,9);
DoubleInterval A11(a2,b2);
DoubleInterval A12(a1,0);
DoubleInterval A21(-1,c1);
DoubleInterval A22(a2,b2);
DoubleMatrix A(2,2);
A(0,0) = A11;
A(0,1) = A12;
A(1,0) = A21;
A(1,1) = A22;
DoubleInterval B1(1,3);
DoubleInterval B2(3,4);
DoubleVector b(2, (DoubleInterval)0);
b[0] = B1;
b[1] = B2;
std::cout << A << std::endl;
std::cout << b << std::endl;
DoubleVector x(2, (DoubleInterval)0);
x = A * b;
std::cout << x << std::endl;
//Need to sort this out
DoubleInterval test(-1,2);
DoubleInterval testb(5,100);
DoubleInterval ans;
ans = testb/test;
std::cout << ans << std::endl;
DoubleInterval k00(0,2);
DoubleInterval k01(1,3);
DoubleInterval k10(3,5);
DoubleInterval k11(5,7);
std::cout << "INVERSE IS " << std::endl;
std::cout << boost::numeric::Doubleinterval_lib::multiplicative_inverse(k11) << std::endl;
std::cout << boost::numeric::norm(k11) << std::endl;
std::cout << k00 * k11 << std::endl;
std::cout << k01 * k10 << std::endl;
std::cout << (k00 * k11) - (k01 * k10) << std::endl;
return 0;
}示例3: 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
}示例4: entry
void LSquare
( DistMatrix<Complex<R> >& A,
DistMatrix<Complex<R>,STAR,STAR>& t )
{
#ifndef RELEASE
CallStackEntry entry("hermitian_tridiag::LSquare");
if( A.Grid() != t.Grid() )
throw std::logic_error("{A,t} must be distributed over the same grid");
#endif
const Grid& g = A.Grid();
#ifndef RELEASE
if( g.Height() != g.Width() )
throw std::logic_error("The process grid must be square");
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;
DistMatrix<C,MD,STAR> tDiag(g);
tDiag.AlignWithDiagonal( A, -1 );
tDiag.ResizeTo( A.Height()-1, 1 );
// Matrix views
DistMatrix<C>
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> 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), A11_MC_STAR(g),
A21_MC_STAR(g);
DistMatrix<C,MR, STAR> APan_MR_STAR(g), A11_MR_STAR(g),
A21_MR_STAR(g);
DistMatrix<C,MC, STAR> WPan_MC_STAR(g), W11_MC_STAR(g),
W21_MC_STAR(g);
DistMatrix<C,MR, STAR> WPan_MR_STAR(g), W11_MR_STAR(g),
W21_MR_STAR(g);
PartitionDownDiagonal
( A, ATL, ATR,
ABL, ABR, 0 );
PartitionDown
( tDiag, tT,
tB, 0 );
while( ATL.Height() < A.Height() )
{
RepartitionDownDiagonal
( ATL, /**/ ATR, A00, /**/ A01, A02,
/*************/ /******************/
/**/ A10, /**/ A11, A12,
ABL, /**/ ABR, A20, /**/ A21, A22 );
RepartitionDown
( tT, t0,
/**/ /**/
t1,
tB, t2 );
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, t1,
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, ADJOINT, ADJOINT,
C(-1), A21_MC_STAR, W21_MR_STAR,
//.........這裏部分代碼省略.........
示例5: 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 );
}
}示例6: 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,
//.........這裏部分代碼省略.........
示例7: 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 );
//.........這裏部分代碼省略.........
示例8: 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;
//.........這裏部分代碼省略.........
示例9: 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 );
//.........這裏部分代碼省略.........
示例10: 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];
}
//.........這裏部分代碼省略.........
示例11: PushCallStack
inline void
internal::HermitianTridiagU( DistMatrix<R,MC,MR>& A )
{
#ifndef RELEASE
PushCallStack("internal::HermitianTridiagU");
if( A.Height() != A.Width() )
throw std::logic_error( "A must be square." );
#endif
const Grid& g = A.Grid();
if( g.InGrid() )
{
// Matrix views
DistMatrix<R,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);
// Temporary distributions
DistMatrix<R,MC, MR > WPan(g);
DistMatrix<R,STAR,STAR> A11_STAR_STAR(g);
DistMatrix<R,MC, STAR> APan_MC_STAR(g), A01_MC_STAR(g),
A11_MC_STAR(g);
DistMatrix<R,MR, STAR> APan_MR_STAR(g), A01_MR_STAR(g),
A11_MR_STAR(g);
DistMatrix<R,MC, STAR> WPan_MC_STAR(g), W01_MC_STAR(g),
W11_MC_STAR(g);
DistMatrix<R,MR, STAR> WPan_MR_STAR(g), W01_MR_STAR(g),
W11_MR_STAR(g);
PartitionUpDiagonal
( A, ATL, ATR,
ABL, ABR, 0 );
while( ABR.Height() < A.Height() )
{
RepartitionUpDiagonal
( ATL, /**/ ATR, A00, A01, /**/ A02,
/**/ A10, A11, /**/ A12,
/*************/ /******************/
ABL, /**/ ABR, A20, A21, /**/ A22 );
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,
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, TRANSPOSE, TRANSPOSE,
(R)-1, A01_MC_STAR, W01_MR_STAR,
W01_MC_STAR, A01_MR_STAR,
(R)1, A00 );
//------------------------------------------------------------//
WPan_MR_STAR.FreeAlignments();
APan_MR_STAR.FreeAlignments();
WPan_MC_STAR.FreeAlignments();
APan_MC_STAR.FreeAlignments();
WPan.FreeAlignments();
}
else
{
A11_STAR_STAR = A11;
HermitianTridiag( UPPER, A11_STAR_STAR.LocalMatrix() );
A11 = A11_STAR_STAR;
}
SlidePartitionUpDiagonal
( ATL, /**/ ATR, A00, /**/ A01, A02,
/*************/ /******************/
/**/ A10, /**/ A11, A12,
ABL, /**/ ABR, A20, /**/ A21, A22 );
}
}
#ifndef RELEASE
PopCallStack();
//.........這裏部分代碼省略.........
示例12: PushCallStack
inline void
internal::CholeskyUVar3Square( DistMatrix<F,MC,MR>& A )
{
#ifndef RELEASE
PushCallStack("internal::CholeskyUVar3Square");
if( A.Height() != A.Width() )
throw std::logic_error
("Can only compute Cholesky factor of square matrices.");
if( A.Grid().Height() != A.Grid().Width() )
throw std::logic_error
("CholeskyUVar3Square assumes a square process grid.");
#endif
const Grid& g = A.Grid();
// Find the process holding our transposed data
const int r = g.Height();
int transposeRank;
{
const int colAlignment = A.ColAlignment();
const int rowAlignment = A.RowAlignment();
const int colShift = A.ColShift();
const int rowShift = A.RowShift();
const int transposeRow = (colAlignment+rowShift) % r;
const int transposeCol = (rowAlignment+colShift) % r;
transposeRank = transposeRow + r*transposeCol;
}
const bool onDiagonal = ( transposeRank == g.VCRank() );
// 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);
// Temporary matrix distributions
DistMatrix<F,STAR,STAR> A11_STAR_STAR(g);
DistMatrix<F,STAR,VR > A12_STAR_VR(g);
DistMatrix<F,STAR,MC > A12_STAR_MC(g);
DistMatrix<F,STAR,MR > A12_STAR_MR(g);
// Start the algorithm
PartitionDownDiagonal
( A, ATL, ATR,
ABL, ABR, 0 );
while( ABR.Height() > 0 )
{
RepartitionDownDiagonal
( ATL, /**/ ATR, A00, /**/ A01, A02,
/*************/ /******************/
/**/ A10, /**/ A11, A12,
ABL, /**/ ABR, A20, /**/ A21, A22 );
A12_STAR_MC.AlignWith( A22 );
A12_STAR_MR.AlignWith( A22 );
A12_STAR_VR.AlignWith( A22 );
//--------------------------------------------------------------------//
A11_STAR_STAR = A11;
internal::LocalCholesky( UPPER, A11_STAR_STAR );
A11 = A11_STAR_STAR;
A12_STAR_VR = A12;
internal::LocalTrsm
( LEFT, UPPER, ADJOINT, NON_UNIT, (F)1, A11_STAR_STAR, A12_STAR_VR );
A12_STAR_MR = A12_STAR_VR;
// SendRecv to form A12[* ,MC] from A12[* ,MR]
A12_STAR_MC.ResizeTo( A12.Height(), A12.Width() );
{
if( onDiagonal )
{
const int size = A11.Height()*A22.LocalWidth();
MemCopy
( A12_STAR_MC.LocalBuffer(),
A12_STAR_MR.LocalBuffer(), size );
}
else
{
const int sendSize = A11.Height()*A22.LocalWidth();
const int recvSize = A11.Width()*A22.LocalHeight();
// We know that the ldim is the height since we have manually
// created both temporary matrices.
mpi::SendRecv
( A12_STAR_MR.LocalBuffer(), sendSize, transposeRank, 0,
A12_STAR_MC.LocalBuffer(), recvSize, transposeRank, 0,
g.VCComm() );
}
}
internal::LocalTrrk
( UPPER, ADJOINT, (F)-1, A12_STAR_MC, A12_STAR_MR, (F)1, A22 );
A12 = A12_STAR_MR;
//--------------------------------------------------------------------//
A12_STAR_MC.FreeAlignments();
A12_STAR_MR.FreeAlignments();
A12_STAR_VR.FreeAlignments();
SlidePartitionDownDiagonal
( ATL, /**/ ATR, A00, A01, /**/ A02,
/**/ A10, A11, /**/ A12,
/*************/ /******************/
//.........這裏部分代碼省略.........
示例13: 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
//.........這裏部分代碼省略.........
示例14: 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)
//.........這裏部分代碼省略.........
示例15: 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
//.........這裏部分代碼省略.........