当前位置: 首页>>代码示例>>C++>>正文


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

本文整理汇总了C++中MatrixView::stepj方法的典型用法代码示例。如果您正苦于以下问题:C++ MatrixView::stepj方法的具体用法?C++ MatrixView::stepj怎么用?C++ MatrixView::stepj使用的例子?那么, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在MatrixView的用法示例。


在下文中一共展示了MatrixView::stepj方法的2个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的C++代码示例。

示例1: LAPV

    template <> void LapHessenberg(
        MatrixView<double> A, VectorView<double> Ubeta)
    {
        TMVAssert(A.iscm());
        TMVAssert(A.colsize() == A.rowsize());
        TMVAssert(Ubeta.size() == A.rowsize()-1);
        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;
        double* work = LAP_DWork(lwork);
#endif
        LAPNAME(dgehrd) (
            LAPCM LAPV(n),LAPV(ilo),LAPV(ihi),
            LAPP(A.ptr()),LAPV(lda),LAPP(Ubeta.ptr())
            LAPWK(work) LAPVWK(lwork) LAPINFO);
#ifdef LAPNOWORK
        LAP_Results(Lap_info,"dgehrd");
#else
        LAP_Results(Lap_info,int(work[0]),m,n,lwork,"dgehrd");
#endif
    }
开发者ID:rmjarvis,项目名称:tmv,代码行数:27,代码来源:TMV_Eigen_Hessenberg.cpp

示例2: BlockHessenberg


//.........这里部分代码省略.........
                //
                // 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:
                //
                // mYZt(0:N,j) = m(0:N,0:N) Y(0:N,0:j) Zt(0:j,j)
                //             = m(0:N,j+1:N) Y(j+1:N,j) Zt(j,j)
                //             = bu* m(0:N,j+1:N) u 
                //
                mYZt.col(jj) = CONJ(bu) * A.subMatrix(j1,N,j+1,N) * u;

                // Restore Aorig, which is actually part of the Hessenberg matrix.
                Atemp = Aorig;
            }

            // Update the rest of the matrix:
            // A(j2,j2-1) needs to be temporarily changed to 1 for use in Y
            T& Atemp = *(A.ptr() + j2*A.stepi() + (j2-1)*A.stepj());
            TMVAssert(IMAG(Atemp) == RealType(T)(0));
            RealType(T) Aorig = Atemp;
            Atemp = RealType(T)(1);

            // m' = (I-YZYt) (m - mYZt Y)
            MatrixView<T> m = A.subMatrix(j1,N,j2,N);
            ConstMatrixView<T> Y = A.subMatrix(j2+1,N,j1,j2);
            m -= mYZt * Y.adjoint();
            BlockHouseholder_LMult(Y,Z,m);

            // Restore A(j2,j2-1)
            Atemp = Aorig;
            j1 = j2;
        }

#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();
        U(0,0) = T(1);
        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;
            Matrix<T,ColMajor> A2 = A0;
            Vector<T> Ub2(Ubeta.size());
            NonBlockHessenberg(A2.view(),Ub2.view());
            cerr<<"cf NonBlock: A -> "<<A2<<endl;
            cerr<<"Ubeta = "<<Ub2<<endl;
            abort();
        }
#endif
    }
开发者ID:rmjarvis,项目名称:tmv,代码行数:101,代码来源:TMV_Eigen_Hessenberg.cpp


注:本文中的MatrixView::stepj方法示例由纯净天空整理自Github/MSDocs等开源代码及文档管理平台,相关代码片段筛选自各路编程大神贡献的开源项目,源码版权归原作者所有,传播和使用请参考对应项目的License;未经允许,请勿转载。