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

    template <> void LapHessenberg(
        MatrixView<std::complex<float> > A, 
        VectorView<std::complex<float> > Ubeta) 
    {
        TMVAssert(A.iscm());
        TMVAssert(A.colsize() >= A.rowsize());
        TMVAssert(Ubeta.size() == A.rowsize());
        TMVAssert(A.ct()==NonConj);

        int n = A.rowsize();
        int ilo = 1;
        int ihi = n;
        int lda = A.stepj();
        int Lap_info=0;
#ifndef LAPNOWORK
        int lwork = n*LAP_BLOCKSIZE;
        std::complex<float>* work = LAP_CWork(lwork);
#endif
        LAPNAME(cgehrd) (
            LAPCM LAPV(n),LAPV(ilo),LAPV(ihi),
            LAPP(A.ptr()),LAPV(lda),LAPP(Ubeta.ptr())
            LAPWK(work) LAPVWK(lwork) LAPINFO);
        Ubeta.ConjugateSelf();
#ifdef LAPNOWORK
        LAP_Results(Lap_info,"cgehrd");
#else
        LAP_Results(Lap_info,int(REAL(work[0])),m,n,lwork,"cgehrd");
#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:
//.........这里部分代码省略.........