本文整理汇总了C++中MatrixView::stepi方法的典型用法代码示例。如果您正苦于以下问题:C++ MatrixView::stepi方法的具体用法?C++ MatrixView::stepi怎么用?C++ MatrixView::stepi使用的例子?那么, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在类MatrixView
的用法示例。
在下文中一共展示了MatrixView::stepi方法的1个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: 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
}