C++ PopCallStack函数代码示例

inline void
Zero( Matrix<T>& A )
{
#ifndef RELEASE
    PushCallStack("Zero");
#endif
    const int height = A.Height();
    const int width = A.Width();
#ifdef HAVE_OPENMP
    #pragma omp parallel for
#endif
    for( int j=0; j<width; ++j )
        MemZero( A.Buffer(0,j), height );
#ifndef RELEASE
    PopCallStack();
#endif
}

inline void
Pseudoinverse( DistMatrix<F>& A )
{
#ifndef RELEASE
    PushCallStack("Pseudoinverse");
#endif
    typedef typename Base<F>::type R;

    const Grid& g = A.Grid();
    const int m = A.Height();
    const int n = A.Width();
    const int k = std::max(m,n);

    // Get the SVD of A
    DistMatrix<R,VR,STAR> s(g);
    DistMatrix<F> U(g), V(g);
    U = A;
    SVD( U, s, V );

    // Compute the two-norm of A as the maximum singular value
    const R twoNorm = Norm( s, INFINITY_NORM );

    // Set the tolerance equal to k ||A||_2 eps and invert above tolerance
    const R eps = lapack::MachineEpsilon<R>();
    const R tolerance = k*twoNorm*eps;
    const int numLocalVals = s.LocalHeight();
    for( int iLocal=0; iLocal<numLocalVals; ++iLocal )
    {
        const R sigma = s.GetLocal(iLocal,0);
        if( sigma < tolerance )
            s.SetLocal(iLocal,0,0);
        else
            s.SetLocal(iLocal,0,1/sigma);
    }

    // Scale U with the singular values, U := U Sigma
    DiagonalScale( RIGHT, NORMAL, s, U );

    // Form pinvA = (U Sigma V^H)^H = V (U Sigma)^H
    Zeros( n, m, A );
    Gemm( NORMAL, ADJOINT, F(1), V, U, F(0), A );
#ifndef RELEASE
    PopCallStack();
#endif
}

inline void
Walsh( int k, Matrix<T>& A, bool binary )
{
#ifndef RELEASE
    PushCallStack("Walsh");
#endif
    if( k < 1 )
        throw std::logic_error("Walsh matrices are only defined for k>=1");

    const unsigned n = 1u<<k;
    A.ResizeTo( n, n );

    // Run a simple O(n^2 log n) algorithm for computing the entries
    // based upon successive sign flips
    const T onValue = 1;
    const T offValue = ( binary ? 0 : -1 );
    for( unsigned j=0; j<n; ++j )
    {
        for( unsigned i=0; i<n; ++i )
        {
            // Recurse on the quadtree, flipping the sign of the entry each
            // time we are in the bottom-right quadrant
            unsigned r = i;     
            unsigned s = j;
            unsigned t = n;
            bool on = true;
            while( t != 1u )
            {
                t >>= 1;
                if( r >= t && s >= t )
                    on = !on;
                r %= t;
                s %= t;
            }

            if( on )
                A.Set( i, j, onValue );
            else
                A.Set( i, j, offValue );
        }
    }
#ifndef RELEASE
    PopCallStack();
#endif
}

inline void
LQ( Matrix<Real>& A )
{
#ifndef RELEASE
    PushCallStack("LQ");
#endif
    if( IsComplex<Real>::val )
        throw std::logic_error("Called real routine with complex datatype");

    // Matrix views
    Matrix<Real>
        ATL, ATR,  A00, A01, A02,  ATopPan, ABottomPan,
        ABL, ABR,  A10, A11, A12,
                   A20, A21, A22;

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

        View1x2( ATopPan, A11, A12 );
        View1x2( ABottomPan, 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
}

inline void
HPSDSquareRoot( UpperOrLower uplo, DistMatrix<R,MC,MR>& A )
{
#ifndef RELEASE
    PushCallStack("HPSDSquareRoot");
#endif
    // Get the EVD of A
    const Grid& g = A.Grid();
    DistMatrix<R,VR,STAR> w(g);
    DistMatrix<R,MC,MR> Z(g);
    HermitianEig( uplo, A, w, Z );

    // Compute the two-norm of A as the maximum absolute value 
    // of its eigenvalues
    R maxLocalAbsEig = 0;
    const int localHeight = w.LocalHeight();
    for( int iLocal=0; iLocal<localHeight; ++iLocal )
        maxLocalAbsEig = 
            std::max(maxLocalAbsEig,Abs(w.GetLocalEntry(iLocal,0)));
    R twoNorm;
    mpi::AllReduce( &maxLocalAbsEig, &twoNorm, 1, mpi::MAX, g.VCComm() );

    // Compute the smallest eigenvalue of A
    R minLocalEig = twoNorm;
    for( int iLocal=0; iLocal<localHeight; ++iLocal )
        minLocalEig = std::min(minLocalEig,w.GetLocalEntry(iLocal,0));
    R minEig;
    mpi::AllReduce( &minLocalEig, &minEig, 1, mpi::MIN, g.VCComm() );

    // Set the tolerance equal to n ||A||_2 eps
    const int n = A.Height();
    const R eps = lapack::MachineEpsilon<R>();
    const R tolerance = n*twoNorm*eps;

    // Ensure that the minimum eigenvalue is not less than - n ||A||_2 eps
    if( minEig < -tolerance )
        throw NonHPSDMatrixException();

    // Form the pseudoinverse
    square_root::Functor<R> f( tolerance );
    hermitian_function::ReformHermitianMatrix( uplo, A, w, Z, f );
#ifndef RELEASE
    PopCallStack();
#endif
}

inline typename Base<F>::type 
LogDetDivergence( UpperOrLower uplo, const Matrix<F>& A, const Matrix<F>& B )
{
#ifndef RELEASE
    PushCallStack("LogDetDivergence");
#endif
    if( A.Height() != A.Width() || B.Height() != B.Width() ||
        A.Height() != B.Height() )
        throw std::logic_error
        ("A and B must be square matrices of the same size");

    typedef typename Base<F>::type R;
    const int n = A.Height();

    Matrix<F> ACopy( A );
    Matrix<F> BCopy( B );

    Cholesky( uplo, ACopy );
    Cholesky( uplo, BCopy );

    if( uplo == LOWER )
    {
        Trtrsm( LEFT, uplo, NORMAL, NON_UNIT, F(1), BCopy, ACopy );
    }
    else
    {
        MakeTrapezoidal( LEFT, uplo, 0, ACopy );
        Trsm( LEFT, uplo, NORMAL, NON_UNIT, F(1), BCopy, ACopy );
    }

    MakeTrapezoidal( LEFT, uplo, 0, ACopy );
    const R frobNorm = Norm( ACopy, FROBENIUS_NORM );

    Matrix<F> d;
    ACopy.GetDiagonal( d );
    R logDet(0);
    for( int i=0; i<n; ++i )
        logDet += 2*Log( RealPart(d.Get(i,0)) );

    const R logDetDiv = frobNorm*frobNorm - logDet - R(n);
#ifndef RELEASE
    PopCallStack();
#endif
    return logDetDiv;
}

inline void
TrrkNNKernel
( UpperOrLower uplo, 
  T alpha, const Matrix<T>& A, const Matrix<T>& B,
  T beta,        Matrix<T>& C )
{
#ifndef RELEASE
    PushCallStack("TrrkNNKernel");
    CheckInputNN( A, B, C );
#endif
    Matrix<T> AT,
              AB;
    Matrix<T> BL, BR;
    Matrix<T> CTL, CTR,
              CBL, CBR;
    Matrix<T> DTL, DBR;

    const int half = C.Height()/2;
    Scale( beta, C );
    LockedPartitionDown
    ( A, AT,
         AB, half );
    LockedPartitionRight( B, BL, BR, half );
    PartitionDownDiagonal
    ( C, CTL, CTR,
         CBL, CBR, half );

    DTL.ResizeTo( CTL.Height(), CTL.Width() );
    DBR.ResizeTo( CBR.Height(), CBR.Width() );
    //------------------------------------------------------------------------//
    if( uplo == LOWER )
        Gemm( NORMAL, NORMAL, alpha, AB, BL, T(1), CBL );
    else
        Gemm( NORMAL, NORMAL, alpha, AT, BR, T(1), CTR );

    Gemm( NORMAL, NORMAL, alpha, AT, BL, T(0), DTL );
    AxpyTriangle( uplo, T(1), DTL, CTL );

    Gemm( NORMAL, NORMAL, alpha, AB, BR, T(0), DBR );
    AxpyTriangle( uplo, T(1), DBR, CBR );
    //------------------------------------------------------------------------//
#ifndef RELEASE
    PopCallStack();
#endif
}

inline void
Hankel( int m, int n, const std::vector<T>& a, Matrix<T>& A )
{
#ifndef RELEASE
    PushCallStack("Hankel");
#endif
    const int length = m+n-1;
    if( a.size() != (unsigned)length )
        throw std::logic_error("a was the wrong size");
    A.ResizeTo( m, n );

    for( int j=0; j<n; ++j )
        for( int i=0; i<n; ++i )
            A.Set( i, j, a[i+j] );
#ifndef RELEASE
    PopCallStack();
#endif
}

inline void
internal::TrmmLUN
( UnitOrNonUnit diag,
  T alpha, const DistMatrix<T,MC,MR>& U,
                 DistMatrix<T,MC,MR>& X )
{
#ifndef RELEASE
    PushCallStack("internal::TrmmLUN");
#endif
    // TODO: Come up with a better routing mechanism
    if( U.Height() > 5*X.Width() )
        internal::TrmmLUNA( diag, alpha, U, X );
    else
        internal::TrmmLUNC( diag, alpha, U, X );
#ifndef RELEASE
    PopCallStack();
#endif
}

inline typename Base<F>::type
SymmetricNorm( UpperOrLower uplo, const DistMatrix<F>& A, NormType type )
{
#ifndef RELEASE
    PushCallStack("SymmetricNorm");
#endif
    typename Base<F>::type norm = 0;
    if( type == NUCLEAR_NORM )
        norm = internal::SymmetricNuclearNorm( uplo, A );
    else if( type == TWO_NORM )
        norm = internal::SymmetricTwoNorm( uplo, A );
    else
        norm = HermitianNorm( uplo, A );
#ifndef RELEASE
    PopCallStack();
#endif
    return norm;
}

inline typename Base<F>::type
EntrywiseOneNorm( const Matrix<F>& A )
{
#ifndef RELEASE
    PushCallStack("EntrywiseOneNorm");
#endif
    typedef typename Base<F>::type R;
    R norm = 0;
    const int width = A.Width();
    const int height = A.Height();
    for( int j=0; j<width; ++j )
        for( int i=0; i<height; ++i )
            norm += Abs(A.Get(i,j));
#ifndef RELEASE
    PopCallStack();
#endif
    return norm;
}

inline void
DistMatrix<T,MD,STAR,Int>::ResizeTo( Int height, Int width )
{
#ifndef RELEASE
    PushCallStack("[MD,* ]::ResizeTo");
    this->AssertNotLockedView();
    if( height < 0 || width < 0 )
        throw std::logic_error("Height and width must be non-negative");
#endif
    this->height_ = height;
    this->width_ = width;
    if( this->Participating() )
        this->localMatrix_.ResizeTo
        ( LocalLength(height,this->ColShift(),this->Grid().LCM()), width );
#ifndef RELEASE
    PopCallStack();
#endif
}

inline void
SymmLL
( T alpha, const DistMatrix<T,MC,MR>& A,
           const DistMatrix<T,MC,MR>& B,
  T beta,        DistMatrix<T,MC,MR>& C )
{
#ifndef RELEASE
    PushCallStack("internal::SymmLL");
#endif
    // TODO: Come up with a better routing mechanism
    if( A.Height() > 5*B.Width() )
        SymmLLA( alpha, A, B, beta, C );
    else
        SymmLLC( alpha, A, B, beta, C );
#ifndef RELEASE
    PopCallStack();
#endif
}

inline void
TrmmRLN
( UnitOrNonUnit diag,
  T alpha, const DistMatrix<T>& L,
                 DistMatrix<T>& X )
{
#ifndef RELEASE
    PushCallStack("internal::TrmmRLN");
#endif
    // TODO: Come up with a better routing mechanism
    if( L.Height() > 5*X.Height() )
        TrmmRLNA( diag, alpha, L, X );
    else
        TrmmRLNC( diag, alpha, L, X );
#ifndef RELEASE
    PopCallStack();
#endif
}

inline void
SymmRU
( T alpha, const DistMatrix<T>& A, const DistMatrix<T>& B,
  T beta,        DistMatrix<T>& C,
  bool conjugate=false )
{
#ifndef RELEASE
    PushCallStack("internal::SymmRU");
#endif
    // TODO: Come up with a better routing mechanism
    if( A.Height() > 5*B.Height() )
        SymmRUA( alpha, A, B, beta, C, conjugate );
    else
        SymmRUC( alpha, A, B, beta, C, conjugate );
#ifndef RELEASE
    PopCallStack();
#endif
}