本文整理汇总了C++中VectorView::step方法的典型用法代码示例。如果您正苦于以下问题:C++ VectorView::step方法的具体用法?C++ VectorView::step怎么用?C++ VectorView::step使用的例子?那么恭喜您, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在类VectorView
的用法示例。
在下文中一共展示了VectorView::step方法的5个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: NonBlockHessenberg
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
}
示例2: Hessenberg
static inline void Hessenberg(
MatrixView<T> A, VectorView<T> Ubeta)
{
TMVAssert(A.colsize() == A.rowsize());
TMVAssert(Ubeta.size() == A.rowsize()-1);
TMVAssert(A.isrm() || A.iscm());
TMVAssert(A.ct()==NonConj);
TMVAssert(Ubeta.step() == 1);
if (A.rowsize() > 0) {
#ifdef LAP
if (A.iscm())
LapHessenberg(A,Ubeta);
else
#endif
NonLapHessenberg(A,Ubeta);
}
}
示例3: MultMV
void MultMV(
const T alpha, const GenDiagMatrix<Ta>& A, const GenVector<Tx>& x,
VectorView<T> y)
// y (+)= alpha * A * x
// yi (+)= alpha * Ai * xi
{
TMVAssert(A.size() == x.size());
TMVAssert(A.size() == y.size());
#ifdef XDEBUG
//cout<<"MultMV: \n";
//cout<<"alpha = "<<alpha<<endl;
//cout<<"A = "<<TMV_Text(A)<<" "<<A<<endl;
//cout<<"x = "<<TMV_Text(x)<<" "<<x<<endl;
//cout<<"y = "<<TMV_Text(y)<<" "<<y<<endl;
Vector<T> y0 = y;
Vector<Tx> x0 = x;
Matrix<Ta> A0 = A;
Vector<T> y2 = alpha*A0*x0;
if (add) y2 += y0;
#endif
ElemMultVV<add>(alpha,A.diag(),x,y);
#ifdef XDEBUG
if (!(Norm(y-y2) <=
0.001*(TMV_ABS(alpha)*Norm(A0)*Norm(x0)+
(add?Norm(y0):TMV_RealType(T)(0))))) {
cerr<<"MultMV: alpha = "<<alpha<<endl;
cerr<<"add = "<<add<<endl;
cerr<<"A = "<<TMV_Text(A)<<" step "<<A.diag().step()<<" "<<A0<<endl;
cerr<<"x = "<<TMV_Text(x)<<" step "<<x.step()<<" "<<x0<<endl;
cerr<<"y = "<<TMV_Text(y)<<" step "<<y.step()<<" "<<y0<<endl;
cerr<<"-> y = "<<y<<endl;
cerr<<"y2 = "<<y2<<endl;
abort();
}
#endif
}
示例4: BlockHessenberg
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:
//.........这里部分代码省略.........
示例5: DoNonLapCH_Decompose
static void DoNonLapCH_Decompose(SymBandMatrixView<T> A)
{
// Cholesky decompostion for a banded Hermitian matrix follows the
// same structure as a regular Cholesky decomposition, but we
// take advantage of the fact that the column or row lengths are
// shorter.
//
TMVAssert(A.uplo() == Lower);
TMVAssert(A.ct() == NonConj);
TMVAssert(isReal(T()) || A.isherm());
TMVAssert(A.iscm() || A.isrm());
TMVAssert(cm == A.iscm());
const ptrdiff_t N = A.size();
#ifdef XDEBUG
Matrix<T> A0(A);
//cout<<"CHDecompose:\n";
//cout<<"A = "<<TMV_Text(A)<<" "<<A<<endl;
#endif
VectorView<TMV_RealType(T)> Adiag = A.diag().realPart();
const ptrdiff_t nlo = A.nlo();
if (nlo == 0) {
TMV_RealType(T)* Ajj= Adiag.ptr();
const ptrdiff_t ds = Adiag.step();
for(ptrdiff_t j=0;j<N;++j,Ajj+=ds) {
#ifdef TMVFLDEBUG
TMVAssert(Ajj >= A.realPart()._first);
TMVAssert(Ajj < A.realPart()._last);
#endif
if (*Ajj <= TMV_RealType(T)(0)) {
#ifdef NOTHROW
std::cerr<<"Non Posdef HermBandMatrix found\n";
exit(1);
#else
throw NonPosDefHermBandMatrix<T>(A);
#endif
}
*Ajj = TMV_SQRT(*Ajj);
}
} else if (cm) {
TMV_RealType(T)* Ajj= Adiag.ptr();
const ptrdiff_t ds = Adiag.step();
ptrdiff_t endcol = nlo+1;
for(ptrdiff_t j=0;j<N-1;++j,Ajj+=ds) {
#ifdef TMVFLDEBUG
TMVAssert(Ajj >= A.realPart()._first);
TMVAssert(Ajj < A.realPart()._last);
#endif
if (*Ajj <= TMV_RealType(T)(0)) {
#ifdef NOTHROW
std::cerr<<"Non Posdef HermBandMatrix found\n";
exit(1);
#else
throw NonPosDefHermBandMatrix<T>(A);
#endif
}
*Ajj = TMV_SQRT(*Ajj);
A.col(j,j+1,endcol) /= *Ajj;
A.subSymMatrix(j+1,endcol) -=
A.col(j,j+1,endcol) ^ A.col(j,j+1,endcol).conjugate();
if (endcol < N) ++endcol;
}
#ifdef TMVFLDEBUG
TMVAssert(Ajj >= A.realPart()._first);
TMVAssert(Ajj < A.realPart()._last);
#endif
if (*Ajj <= TMV_RealType(T)(0)) {
#ifdef NOTHROW
std::cerr<<"Non Posdef HermBandMatrix found\n";
exit(1);
#else
throw NonPosDefHermBandMatrix<T>(A);
#endif
}
*Ajj = TMV_SQRT(*Ajj);
} else {
TMV_RealType(T)* Aii = Adiag.ptr();
const ptrdiff_t ds = Adiag.step();
ptrdiff_t startrow = 0;
#ifdef TMVFLDEBUG
TMVAssert(Aii >= A.realPart()._first);
TMVAssert(Aii < A.realPart()._last);
#endif
if (*Aii <= TMV_RealType(T)(0)) {
#ifdef NOTHROW
std::cerr<<"Non Posdef HermBandMatrix found\n";
exit(1);
#else
throw NonPosDefHermBandMatrix<T>(A);
#endif
}
*Aii = TMV_SQRT(*Aii);
for(ptrdiff_t i=1;i<N;++i) {
if (i > nlo) ++startrow;
Aii+=ds;
A.row(i,startrow,i) %=
A.subSymBandMatrix(startrow,i).lowerBand().adjoint();
#ifdef TMVFLDEBUG
TMVAssert(Aii >= A.realPart()._first);
//.........这里部分代码省略.........