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

    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:
//.........这里部分代码省略.........

//##ModelId=40A1898902EF
bool LoadBasedScreen1::CalculateModel( PFlowStream1 FeedStream )
  {
  // Local Variables
  int i = 0;
  int j = 0;
  double ScaleFactor = 0.00;
  double Denom = 0.00;

  Vector CombCPPSizing ( nSize );
  CubicSpline FeedSpline;    

  MatrixView FeedSolids      = FeedStream->AccessSolidsMatrix( );

  MatrixView OversizeSolids  = Oversize->AccessSolidsMatrix( );

  MatrixView UndersizeSolids = Undersize->AccessSolidsMatrix( );

  // Ensure correct shape feed matrix

  if( FeedSolids.rowCount() != nSize )
    goto calcFailed;

  if( FeedSolids.columnCount() != nType ) 
    goto calcFailed;

  // Calculate total feed flow rate (solids)

  QF = FeedSolids.sum( );

  // Convert feed matrix into single component size distribution

  for( i = 0; i < nSize; i++ )
    CombRetSizing[ i ] = (FeedSolids.row(i)).sum();

  if (fabs(QF)<1e-8)
    ScaleFactor = 0.00;
  else 
    ScaleFactor = 1/QF;

  // Scale single component distribution to cumulative sizing
  //  scaled to fraction basis: ie. 1 passes top-size

  CombCPPSizing[ nSize - 1 ] = 0;
  for( i = nSize - 2; i >=0; i-- )
    CombCPPSizing[ i ] = CombCPPSizing[i + 1] + ( CombRetSizing[i + 1] * ScaleFactor );

  // Construct cubic spline passing through cumulative curve    

  FeedSpline.SetSpline( Sizes, CombCPPSizing );

  // Use cubic spline to get fraction passing half apperture size

  Feed_HS = FeedSpline.CalculateY( S/2 );

  // Use spline to get fraction passing apperture size    

  U = Feed_US = FeedSpline.CalculateY( S );
  Feed_OS = 1 - Feed_US;

  // Calculate area factors for this screen

  A = CalculateA( S );
  B = CalculateB( Feed_OS );
  C = CalculateC( Feed_HS );
  D = CalculateD( DL );
  E = CalculateE( WS, S );
  H = CalculateH( OT );
  J = CalculateJ( ST, S, OA );
  K = CalculateK( FT );
  L = CalculateL();
  X = CalculateX( CF );

  // Calculate denominator of load equation - ensure not zero

  Denom = A * B * C * D * E * F * H * J * K * L * AF * X;
  if (fabs(Denom)<1e-8)
    Denom = 1e-8;

  // Calculate load

  V = ( 90 * QF * U ) / ( Denom );

  // Calculate efficiency that matches this load

  T = CalculateT( V );

  // Calculate factor G

  G = ( V * T ) / 90;

  // calculate flow to undersize at this efficiency

  QU = QF * U * T;

  // Solve for D50 that matches the undersize flow

  D50 = zbrent( S * 0.1, S * 0.99, TOLERANCE );

  // Calculate splitting of feed water to products
//.........这里部分代码省略.........