C++ PartitionDownDiagonal函数代码示例

void Sylvester
( Int m,
  Matrix<F>& W,
  Matrix<F>& X,
  SignCtrl<Base<F>> ctrl )
{
    EL_DEBUG_CSE
    Sign( W, ctrl );
    Matrix<F> WTL, WTR,
              WBL, WBR;
    PartitionDownDiagonal
    ( W, WTL, WTR,
         WBL, WBR, m );
    // WTL and WBR should be the positive and negative identity, WBL should be 
    // zero, and WTR should be -2 X
    X = WTR;
    X *= -F(1)/F(2);

    // TODO: Think of how to probe for checks on other quadrants.
    /*
    typedef Base<F> Real; 
    UpdateDiagonal( WTL, F(-1) );
    const Real errorWTL = FrobeniusNorm( WTL );
    const Int n = W.Height() - m;
    UpdateDiagonal( WBR, F(1) );
    const Real errorWBR = FrobeniusNorm( WBR );
    const Real errorWBL = FrobeniusNorm( WBL );
    */
}

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 );
    }
}

void Sylvester
( Int m,
  ElementalMatrix<F>& WPre,
  ElementalMatrix<F>& X, 
  SignCtrl<Base<F>> ctrl )
{
    EL_DEBUG_CSE

    DistMatrixReadProxy<F,F,MC,MR> WProx( WPre );
    auto& W = WProx.Get();

    const Grid& g = W.Grid();
    Sign( W, ctrl );
    DistMatrix<F> WTL(g), WTR(g),
                  WBL(g), WBR(g);
    PartitionDownDiagonal
    ( W, WTL, WTR,
         WBL, WBR, m );
    // WTL and WBR should be the positive and negative identity, WBL should be 
    // zero, and WTR should be -2 X
    Copy( WTR, X );
    X *= -F(1)/F(2);

    // TODO: Think of how to probe for checks on other quadrants.
    //       Add UpdateDiagonal routine to avoid explicit identity Axpy?
    /*
    typedef Base<F> Real; 
    UpdateDiagonal( WTL, F(-1) );
    const Real errorWTL = FrobeniusNorm( WTL );
    const Int n = W.Height() - m;
    UpdateDiagonal( WBR, F(1) );
    const Real errorWBR = FrobeniusNorm( WBR );
    const Real errorWBL = FrobeniusNorm( WBL );
    */
}

inline int
Lyapunov( const DistMatrix<F>& A, const DistMatrix<F>& C, DistMatrix<F>& X )
{
#ifndef RELEASE
    CallStackEntry cse("Sylvester");
    if( A.Height() != A.Width() )
        LogicError("A must be square");
    if( C.Height() != A.Height() || C.Width() != A.Height() )
        LogicError("C must conform with A");
    if( A.Grid() != C.Grid() )
        LogicError("A and C must have the same grid");
#endif
    const Grid& g = A.Grid();
    const Int m = A.Height();
    DistMatrix<F> W(g), WTL(g), WTR(g),
                        WBL(g), WBR(g);
    Zeros( W, 2*m, 2*m );
    PartitionDownDiagonal
    ( W, WTL, WTR,
         WBL, WBR, m );
    WTL = A;
    Adjoint( A, WBR ); Scale( F(-1), WBR );
    WTR = C; Scale( F(-1), WTR );
    return Sylvester( m, W, X );
}

void Riccati
( ElementalMatrix<F>& WPre,
  ElementalMatrix<F>& X, 
  SignCtrl<Base<F>> ctrl )
{
    DEBUG_CSE

    DistMatrixReadProxy<F,F,MC,MR> WProx( WPre );
    auto& W = WProx.Get();

    const Grid& g = W.Grid();
    Sign( W, ctrl );
    const Int n = W.Height()/2;
    DistMatrix<F> WTL(g), WTR(g),
                  WBL(g), WBR(g);
    PartitionDownDiagonal
    ( W, WTL, WTR,
         WBL, WBR, n );

    // (ML, MR) = sgn(W) - I
    ShiftDiagonal( W, F(-1) );

    // Solve for X in ML X = -MR
    DistMatrix<F> ML(g), MR(g);
    PartitionRight( W, ML, MR, n );
    MR *= -1;
    ls::Overwrite( NORMAL, ML, MR, X );
}

inline void
Householder( Matrix<F>& A, Matrix<F>& t )
{
#ifndef RELEASE
    CallStackEntry entry("lq::Householder");
#endif
    t.ResizeTo( Min(A.Height(),A.Width()), 1 );

    // Matrix views
    Matrix<F>
    ATL, ATR,  A00, A01, A02,  ATopPan, ABottomPan,
         ABL, ABR,  A10, A11, A12,
         A20, A21, A22;
    Matrix<F>
    tT,  t0,
    tB,  t1,
    t2;

    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 );

        View1x2( ATopPan, A11, A12 );
        View1x2( ABottomPan, A21, A22 );
        //--------------------------------------------------------------------//
        PanelHouseholder( ATopPan, t1 );
        ApplyQ( RIGHT, ADJOINT, ATopPan, t1, ABottomPan );
        //--------------------------------------------------------------------//

        SlidePartitionDown
        ( tT,  t0,
          t1,
          /**/ /**/
          tB,  t2 );

        SlidePartitionDownDiagonal
        ( ATL, /**/ ATR,  A00, A01, /**/ A02,
          /**/       A10, A11, /**/ A12,
          /*************/ /******************/
          ABL, /**/ ABR,  A20, A21, /**/ A22 );
    }
}

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 );
    }
}

inline void
LocalTrrkKernel
( UpperOrLower uplo,
  Orientation orientationOfA,
  Orientation orientationOfB,
  T alpha, const DistMatrix<T,STAR,MC  >& A,
           const DistMatrix<T,MR,  STAR>& B,
  T beta,        DistMatrix<T,MC,  MR  >& C )
{
#ifndef RELEASE
    PushCallStack("LocalTrrkKernel");
    CheckInput( orientationOfA, orientationOfB, A, B, C );
#endif
    const Grid& g = C.Grid();

    DistMatrix<T,STAR,MC> AL(g), AR(g);
    DistMatrix<T,MR,STAR> BT(g), 
                          BB(g);
    DistMatrix<T,MC,MR> CTL(g), CTR(g),
                        CBL(g), CBR(g);
    DistMatrix<T,MC,MR> DTL(g), DBR(g);

    const int half = C.Height()/2;
    ScaleTrapezoid( beta, LEFT, uplo, 0, C );
    LockedPartitionRight( A, AL, AR, half );
    LockedPartitionDown
    ( B, BT, 
         BB, half );
    PartitionDownDiagonal
    ( C, CTL, CTR,
         CBL, CBR, half );

    DTL.AlignWith( CTL );
    DBR.AlignWith( CBR );
    DTL.ResizeTo( CTL.Height(), CTL.Width() );
    DBR.ResizeTo( CBR.Height(), CBR.Width() );
    //------------------------------------------------------------------------//
    if( uplo == LOWER )
        internal::LocalGemm
        ( orientationOfA, orientationOfB, alpha, AR, BT, T(1), CBL );
    else
        internal::LocalGemm
        ( orientationOfA, orientationOfB, alpha, AL, BB, T(1), CTR );

    internal::LocalGemm
    ( orientationOfA, orientationOfB, alpha, AL, BT, T(0), DTL );
    AxpyTriangle( uplo, T(1), DTL, CTL );

    internal::LocalGemm
    ( orientationOfA, orientationOfB, alpha, AR, BB, T(0), DBR );
    AxpyTriangle( uplo, T(1), DBR, CBR );
    //------------------------------------------------------------------------//
#ifndef RELEASE
    PopCallStack();
#endif
}

inline void HermitianTridiagL( Matrix<R>& A )
{
#ifndef RELEASE
    PushCallStack("HermitianTridiagL");
    if( A.Height() != A.Width() )
        throw std::logic_error("A must be square");
#endif
    // Matrix views 
    Matrix<R>
        ATL, ATR,  A00, a01,     A02,  alpha21T,
        ABL, ABR,  a10, alpha11, a12,  a21B,
                   A20, a21,     A22;

    // Temporary matrices
    Matrix<R> w21;

    PushBlocksizeStack( 1 );
    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 );

        PartitionDown
        ( a21, alpha21T,
               a21B,     1 );

        w21.ResizeTo( a21.Height(), 1 );
        //--------------------------------------------------------------------//
        const R tau = Reflector( alpha21T, a21B );
        const R epsilon1 = alpha21T.Get(0,0);
        alpha21T.Set(0,0,R(1));

        Symv( LOWER, tau, A22, a21, R(0), w21 );
        const R alpha = -tau*Dot( w21, a21 )/R(2);
        Axpy( alpha, a21, w21 );
        Syr2( LOWER, R(-1), a21, w21, A22 );
        alpha21T.Set(0,0,epsilon1);
        //--------------------------------------------------------------------//

        SlidePartitionDownDiagonal
        ( ATL, /**/ ATR,  A00, a01,     /**/ A02,
               /**/       a10, alpha11, /**/ a12,
         /*************/ /**********************/
          ABL, /**/ ABR,  A20, a21,     /**/ A22 );
    }
    PopBlocksizeStack();
#ifndef RELEASE
    PopCallStack();
#endif
}

inline void
RepartitionDownDiagonal
( DM& ATL, DM& ATR, DM& A00, DM& A01, DM& A02,
                    DM& A10, DM& A11, DM& A12,
  DM& ABL, DM& ABR, DM& A20, DM& A21, DM& A22, Int bsize=Blocksize() )
{
    DEBUG_ONLY(CallStackEntry cse("RepartitionDownDiagonal"))
    View( A00, ATL );
    PartitionDownDiagonal( ABR, A11, A12,
                                A21, A22, bsize );
    PartitionDown( ABL, A10, A20, A11.Height() );
    PartitionRight( ATR, A01, A02, A11.Width() );
}

inline void
CholeskyUVar2( Matrix<F>& A )
{
#ifndef RELEASE
    PushCallStack("hpd_inverse::CholeskyUVar2");
    if( A.Height() != A.Width() )
        throw std::logic_error("Nonsquare matrices cannot be triangular");
#endif
    // Matrix views
    Matrix<F> 
        ATL, ATR,  A00, A01, A02,
        ABL, ABR,  A10, A11, A12,
                   A20, A21, A22;

    // Start the algorithm
    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 );

        //--------------------------------------------------------------------//
        Cholesky( UPPER, A11 );
        Trsm( RIGHT, UPPER, NORMAL, NON_UNIT, F(1), A11, A01 );
        Trsm( LEFT, UPPER, ADJOINT, NON_UNIT, F(1), A11, A12 );
        Herk( UPPER, NORMAL, F(1), A01, F(1), A00 );
        Gemm( NORMAL, NORMAL, F(-1), A01, A12, F(1), A02 );
        Herk( UPPER, ADJOINT, F(-1), A12, F(1), A22 );
        Trsm( RIGHT, UPPER, ADJOINT, NON_UNIT, F(1), A11, A01 );
        Trsm( LEFT, UPPER, NORMAL, NON_UNIT, F(-1), A11, A12 );
        TriangularInverse( UPPER, NON_UNIT, A11 );
        Trtrmm( ADJOINT, UPPER, A11 );
        //--------------------------------------------------------------------//

        SlidePartitionDownDiagonal
        ( ATL, /**/ ATR,  A00, A01, /**/ A02,
               /**/       A10, A11, /**/ A12,
         /*************/ /******************/
          ABL, /**/ ABR,  A20, A21, /**/ A22 );
    }
#ifndef RELEASE
    PopCallStack();
#endif
}

inline void
RepartitionDownDiagonal
( DM& ATL, DM& ATR, DM& A00, DM& A01, DM& A02,
                    DM& A10, DM& A11, DM& A12,
  DM& ABL, DM& ABR, DM& A20, DM& A21, DM& A22, Int bsize )
{
#ifndef RELEASE
    CallStackEntry cse("RepartitionDownDiagonal [DistMatrix]");
#endif
    View( A00, ATL );
    PartitionDownDiagonal( ABR, A11, A12,
                                A21, A22, bsize );
    PartitionDown( ABL, A10, A20, A11.Height() );
    PartitionRight( ATR, A01, A02, A11.Width() );
}

inline void
TrdtrmmUVar1( Orientation orientation, Matrix<F>& U )
{
#ifndef RELEASE
    PushCallStack("internal::TrtdrmmUVar1");
    if( U.Height() != U.Width() )
        throw std::logic_error("U must be square");
    if( orientation == NORMAL )
        throw std::logic_error("Orientation must be (conjugate-)transpose");
#endif
    Matrix<F>
        UTL, UTR,  U00, U01, U02,
        UBL, UBR,  U10, U11, U12,
                   U20, U21, U22;
    Matrix<F> d1, S01;

    PartitionDownDiagonal
    ( U, UTL, UTR,
         UBL, UBR, 0 );
    while( UTL.Height() < U.Height() && UTL.Width() < U.Height() )
    {
        RepartitionDownDiagonal
        ( UTL, /**/ UTR,  U00, /**/ U01, U02,
         /*************/ /******************/
               /**/       U10, /**/ U11, U12,
          UBL, /**/ UBR,  U20, /**/ U21, U22 );

        //--------------------------------------------------------------------/
        U11.GetDiagonal( d1 );
        S01 = U01;
        DiagonalSolve( LEFT, NORMAL, d1, U01, true );
        Trrk( UPPER, NORMAL, orientation, F(1), U01, S01, F(1), U00 );
        Trmm( RIGHT, UPPER, ADJOINT, UNIT, F(1), U11, U01 );
        TrdtrmmUUnblocked( orientation, U11 );
        //--------------------------------------------------------------------/

        SlidePartitionDownDiagonal
        ( UTL, /**/ UTR,  U00, U01, /**/ U02,
               /**/       U10, U11, /**/ U12,
         /*************/ /******************/
          UBL, /**/ UBR,  U20, U21, /**/ U22 );
    }
#ifndef RELEASE
    PopCallStack();
#endif
}

inline void
CholeskyLVar2( Matrix<F>& A )
{
#ifndef RELEASE
    PushCallStack("internal::CholeskyLVar2");
    if( A.Height() != A.Width() )
        throw std::logic_error
        ("Can only compute Cholesky factor of square matrices");
#endif
    // Matrix views
    Matrix<F> 
        ATL, ATR,   A00, A01, A02,
        ABL, ABR,   A10, A11, A12,
                    A20, A21, A22;

    // Start the algorithm
    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 );

        //--------------------------------------------------------------------//
        Herk( LOWER, NORMAL, F(-1), A10, F(1), A11 );
        CholeskyLVar3Unb( A11 );
        Gemm( NORMAL, ADJOINT, F(-1), A20, A10, F(1), A21 );
        Trsm( RIGHT, LOWER, ADJOINT, NON_UNIT, F(1), A11, A21 );
        //--------------------------------------------------------------------//

        SlidePartitionDownDiagonal
        ( ATL, /**/ ATR,  A00, A01, /**/ A02,
               /**/       A10, A11, /**/ A12,
         /*************/ /******************/
          ABL, /**/ ABR,  A20, A21, /**/ A22 );
    }
#ifndef RELEASE
    PopCallStack();
#endif
}

inline void
TrdtrmmLVar1( Orientation orientation, Matrix<F>& L )
{
#ifndef RELEASE
    CallStackEntry entry("internal::TrdtrmmLVar1");
    if( L.Height() != L.Width() )
        LogicError("L must be square");
    if( orientation == NORMAL )
        LogicError("Orientation must be (conjugate-)transpose");
#endif
    Matrix<F>
        LTL, LTR,  L00, L01, L02,
        LBL, LBR,  L10, L11, L12,
                   L20, L21, L22;
    Matrix<F> d1, S10;

    PartitionDownDiagonal
    ( L, LTL, LTR,
         LBL, LBR, 0 );
    while( LTL.Height() < L.Height() && LTL.Width() < L.Height() )
    {
        RepartitionDownDiagonal
        ( LTL, /**/ LTR,  L00, /**/ L01, L02,
         /*************/ /******************/
               /**/       L10, /**/ L11, L12,
          LBL, /**/ LBR,  L20, /**/ L21, L22 );

        //--------------------------------------------------------------------/
        L11.GetDiagonal( d1 );
        S10 = L10;
        DiagonalSolve( LEFT, NORMAL, d1, L10, true );
        Trrk( LOWER, orientation, NORMAL, F(1), S10, L10, F(1), L00 );
        Trmm( LEFT, LOWER, orientation, UNIT, F(1), L11, L10 );
        TrdtrmmLUnblocked( orientation, L11 );
        //--------------------------------------------------------------------/

        SlidePartitionDownDiagonal
        ( LTL, /**/ LTR,  L00, L01, /**/ L02,
               /**/       L10, L11, /**/ L12,
         /*************/ /******************/
          LBL, /**/ LBR,  L20, L21, /**/ L22 );
    }
}