C++ MatrixView::subMatrix方法代码示例

    static void NonBlockHessenberg(
        MatrixView<T> A, VectorView<T> Ubeta)
    {
#ifdef XDEBUG
        cout<<"Start NonBlock Hessenberg Reduction: A = "<<A<<endl;
        Matrix<T> A0(A);
#endif
        // Decompose A into U H Ut
        // H is a Hessenberg Matrix
        // U is a Unitary Matrix
        // On output, H is stored in the upper-Hessenberg part of A
        // U is stored in compact form in the rest of A along with 
        // the vector Ubeta.
        const ptrdiff_t N = A.rowsize();

        TMVAssert(A.colsize() == A.rowsize());
        TMVAssert(N > 0);
        TMVAssert(Ubeta.size() == N-1);
        TMVAssert(A.iscm() || A.isrm());
        TMVAssert(!Ubeta.isconj());
        TMVAssert(Ubeta.step()==1);

        // We use Householder reflections to reduce A to the Hessenberg form:
        T* Uj = Ubeta.ptr();
        T det = 0; // Ignore Householder det calculations
        for(ptrdiff_t j=0;j<N-1;++j,++Uj) {
#ifdef TMVFLDEBUG
            TMVAssert(Uj >= Ubeta._first);
            TMVAssert(Uj < Ubeta._last);
#endif
            *Uj = Householder_Reflect(A.subMatrix(j+1,N,j,N),det);
            if (*Uj != T(0))
                Householder_LMult(A.col(j+2,N),*Uj,A.subMatrix(0,N,j+1,N).adjoint());
        }

#ifdef XDEBUG
        Matrix<T> U(N,N,T(0));
        U.subMatrix(1,N,1,N) = A.subMatrix(1,N,0,N-1);
        U.upperTri().setZero();
        Vector<T> Ubeta2(N);
        Ubeta2.subVector(1,N) = Ubeta;
        Ubeta2(0) = T(0);
        GetQFromQR(U.view(),Ubeta2);
        Matrix<T> H = A;
        if (N>2) LowerTriMatrixViewOf(H).offDiag(2).setZero();
        Matrix<T> AA = U*H*U.adjoint();
        if (Norm(A0-AA) > 0.001*Norm(A0)) {
            cerr<<"NonBlock Hessenberg: A = "<<Type(A)<<"  "<<A0<<endl;
            cerr<<"A = "<<A<<endl;
            cerr<<"Ubeta = "<<Ubeta<<endl;
            cerr<<"U = "<<U<<endl;
            cerr<<"H = "<<H<<endl;
            cerr<<"UHUt = "<<AA<<endl;
            abort();
        }
#endif
    }

    void SymSquare(MatrixView<T> A)
    {
        const ptrdiff_t N = A.colsize();
        if (N == 1) {
            const T A00 = *A.ptr();
#ifdef TMVFLDEBUG
            TMVAssert(A.ptr() >= A._first);
            TMVAssert(A.ptr() < A._last);
#endif
            if (herm)
                *A.ptr() = TMV_NORM(TMV_REAL(A00));
            else 
                *A.ptr() = TMV_SQR(A00);
        } else {
            const ptrdiff_t K = N/2;
            MatrixView<T> A00 = A.subMatrix(0,K,0,K);
            MatrixView<T> A10 = A.subMatrix(K,N,0,K);
            MatrixView<T> A01 = A.subMatrix(0,K,K,N);
            MatrixView<T> A11 = A.subMatrix(K,N,K,N);
            MatrixView<T> A10t = herm ? A10.adjoint() : A10.transpose();

            // [ A00 A10t ] [ A00 A10t ] 
            // [ A10 A11  ] [ A10 A11  ]
            // = [ A00^2 + A10t A10    A00 A10t + A10t A11 ]
            //   [ A10 A00 + A11 A10   A10 A10t + A11^2    ]

            // A10 stores the actual data for A10
            // We can therefore write to A01 as a temp matrix.

            A01 = A00 * A10t;
            A01 += A10t * A11;

            SymSquare<herm>(A00);
            A00 += A10t*A10;
            SymSquare<herm>(A11);
            A11 += A10*A10t;

            A10t = A01;
        }
    }

    static void BlockHessenberg(
        MatrixView<T> A, VectorView<T> Ubeta)
    {
        // Much like the block version of Bidiagonalize, we try to maintain
        // the operation of several successive Householder matrices in
        // a block form, where the net Block Householder is I - YZYt.
        //
        // But as with the bidiagonlization algorithm (and unlike a simple
        // block QR decomposition), we update the matrix from both the left 
        // and the right, so we also need to keep track of the product
        // ZYtm in addition.
        //
        // The block update at the end of the block loop is
        // m' = (I-YZYt) m (I-YZtYt)
        //
        // The Y matrix is stored in the first K columns of m,
        // and the Hessenberg portion of these columns is updated as we go.
        // For the right-hand-side update, m -= mYZtYt, the m on the right
        // needs to be the full original matrix m, including the original
        // versions of these K columns.  Therefore, we can't wait until 
        // the end for this calculation.  
        //
        // Instead, we keep track of mYZt as we progress, so the final update
        // is:
        //
        // m' = (I-YZYt) (m - mYZt Y)
        //
        // We also need to do this same calculation for each column as we
        // progress through the block.
        //
        const ptrdiff_t N = A.rowsize();

#ifdef XDEBUG
        Matrix<T> A0(A);
#endif

        TMVAssert(A.rowsize() == A.colsize());
        TMVAssert(N > 0);
        TMVAssert(Ubeta.size() == N-1);
        TMVAssert(!Ubeta.isconj());
        TMVAssert(Ubeta.step()==1);

        ptrdiff_t ncolmax = MIN(HESS_BLOCKSIZE,N-1);
        Matrix<T,RowMajor> mYZt_full(N,ncolmax);
        UpperTriMatrix<T,NonUnitDiag|ColMajor> Z_full(ncolmax);

        T det(0); // Ignore Householder Determinant calculations
        T* Uj = Ubeta.ptr();
        for(ptrdiff_t j1=0;j1<N-1;) {
            ptrdiff_t j2 = MIN(N-1,j1+HESS_BLOCKSIZE);
            ptrdiff_t ncols = j2-j1;
            MatrixView<T> mYZt = mYZt_full.subMatrix(0,N-j1,0,ncols);
            UpperTriMatrixView<T> Z = Z_full.subTriMatrix(0,ncols);

            for(ptrdiff_t j=j1,jj=0;j<j2;++j,++jj,++Uj) { // jj = j-j1

                // Update current column of A
                //
                // m' = (I - YZYt) (m - mYZt Yt)
                // A(0:N,j)' = A(0:N,j) - mYZt(0:N,0:j) Y(j,0:j)t
                A.col(j,j1+1,N) -= mYZt.Cols(0,j) * A.row(j,0,j).Conjugate();
                //
                // A(0:N,j)'' = A(0:N,j) - Y Z Yt A(0:N,j)'
                // 
                // Let Y = (L)     where L is unit-diagonal, lower-triangular,
                //         (M)     and M is rectangular
                //
                LowerTriMatrixView<T> L = 
                    LowerTriMatrixViewOf(A.subMatrix(j1+1,j+1,j1,j),UnitDiag);
                MatrixView<T> M = A.subMatrix(j+1,N,j1,j);
                // Use the last column of Z as temporary storage for Yt A(0:N,j)'
                VectorView<T> YtAj = Z.col(jj,0,jj);
                YtAj = L.adjoint() * A.col(j,j1+1,j+1);
                YtAj += M.adjoint() * A.col(j,j+1,N);
                YtAj = Z.subTriMatrix(0,jj) * YtAj;
                A.col(j,j1+1,j+1) -= L * YtAj;
                A.col(j,j+1,N) -= M * YtAj;

                // Do the Householder reflection 
                VectorView<T> u = A.col(j,j+1,N);
                T bu = Householder_Reflect(u,det);
#ifdef TMVFLDEBUG
                TMVAssert(Uj >= Ubeta._first);
                TMVAssert(Uj < Ubeta._last);
#endif
                *Uj = bu;

                // Save the top of the u vector, which isn't actually part of u
                T& Atemp = *u.cptr();
                TMVAssert(IMAG(Atemp) == RealType(T)(0));
                RealType(T) Aorig = REAL(Atemp);
                Atemp = RealType(T)(1);

                // Update Z
                VectorView<T> Zj = Z.col(jj,0,jj);
                Zj = -bu * M.adjoint() * u;
                Zj = Z * Zj;
                Z(jj,jj) = -bu;

                // Update mYtZt:
//.........这里部分代码省略.........