C++ MatrixView类代码示例

 void SymBandSVDiv<T>::doMakeInverse(MatrixView<T1> minv) const
 { 
     CallSymSV_Inverse(
         T(),pimpl->U,pimpl->S,pimpl->Vt,pimpl->kmax,
         SymMatrixViewOf(minv,Upper));
     if (pimpl->S.size() > 1)
         minv.lowerTri().offDiag() = minv.upperTri().offDiag().transpose();
 }

 void SymBandSVDiv<T>::doLDiv(
     const GenMatrix<T1>& m, MatrixView<T2> x) const
 {
     TMVAssert(m.rowsize() == x.rowsize());
     TMVAssert(m.colsize() == colsize());
     TMVAssert(x.colsize() == rowsize());
     CallSV_LDiv(T(),pimpl->U,pimpl->S,pimpl->Vt,pimpl->kmax,m,x);
 }

 void LUDiv<T>::doRDiv(
     const GenMatrix<T1>& m1, MatrixView<T2> m0) const
 {
     TMVAssert(m1.colsize() == m0.colsize());
     TMVAssert(m1.rowsize() == m0.rowsize());
     TMVAssert(pimpl->LUx.colsize() == m1.rowsize());
     doRDivEq(m0=m1);
 }

 void HermBandSVDiv<T>::doMakeInverse(MatrixView<T1> minv) const
 { 
     // A^-1 = U S^-1 Ut
     if (isComplex(T1())) minv.diag().imagPart().setZero();
     CallHermSV_Inverse(
         T(),pimpl->U,pimpl->S,pimpl->kmax,HermMatrixViewOf(minv,Upper));
     if (pimpl->S.size() > 1)
         minv.lowerTri().offDiag() = minv.upperTri().offDiag().adjoint();
 }

    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
    }

    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
    }

    static inline void NonLapHessenberg(
        MatrixView<T> A, VectorView<T> Ubeta)
    {
        TMVAssert(A.rowsize() == A.colsize());
        TMVAssert(A.rowsize() > 0);
        TMVAssert(Ubeta.size() == A.rowsize()-1);

#if 0
        if (A.rowsize() > HESS_BLOCKSIZE)
            BlockHessenberg(A,Ubeta,Vbeta,D,E,det);
        else
#endif
            NonBlockHessenberg(A,Ubeta);
    }

    void LU_Inverse(
        const GenBandMatrix<T1>& LUx, const ptrdiff_t* p, MatrixView<T> minv)
    {
        TMVAssert(LUx.isSquare());
        TMVAssert(minv.isSquare());
        TMVAssert(minv.colsize() == LUx.colsize());
#ifdef XDEBUG
        LowerTriMatrix<T,UnitDiag> L0(LUx.colsize());
        LU_PackedPL_Unpack(LUx,p,L0.view());
        UpperTriMatrix<T> U0 = BandMatrixViewOf(LUx,0,LUx.nhi());
        Matrix<T> PLU = L0 * U0;
        if (LUx.nlo() > 0) PLU.reversePermuteRows(p);
        Matrix<T> minv2 = PLU.inverse();
#endif

        if (minv.colsize() > 0) {
            if ( !(minv.iscm() || minv.isrm())) {
                Matrix<T,ColMajor> temp(minv.colsize(),minv.colsize());
                LU_Inverse(LUx,p,temp.view());
                minv = temp;
            } else {
                minv.setZero();
                UpperTriMatrixView<T> U = minv.upperTri();
                U = BandMatrixViewOf(LUx,0,LUx.nhi());
                TriInverse(U,LUx.nhi());
                LU_PackedPL_RDivEq(LUx,p,minv);
            }
        }

#ifdef XDEBUG
        TMV_RealType(T) normdiff = Norm(PLU*minv - T(1));
        TMV_RealType(T) kappa = Norm(PLU)*Norm(minv);
        if (normdiff > 0.001*kappa*minv.colsize()) {
            cerr<<"LUInverse:\n";
            cerr<<"LUx = "<<LUx<<endl;
            cerr<<"p = ";
            for(ptrdiff_t i=0;i<LUx.colsize();i++) cerr<<p[i]<<" ";
            cerr<<endl;
            cerr<<"PLU = "<<PLU<<endl;
            cerr<<"minv = "<<minv<<endl;
            cerr<<"minv2 = "<<minv2<<endl;
            cerr<<"m*minv = "<<PLU*minv<<endl;
            cerr<<"minv*m = "<<minv*PLU<<endl;
            cerr<<"Norm(m*minv - 1) = "<<normdiff<<endl;
            cerr<<"kappa = "<<kappa<<endl;
            abort();
        }
#endif
    }

void seissol::model::applyBoundaryConditionToElasticFluxSolver( enum ::faceType type,
                                                                MatrixView<NUMBER_OF_QUANTITIES, 9> Fneighbor )
{
  if (type == freeSurface) {
    // Gamma is a diagonal matrix
    real Gamma[] = { -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0 };
    // Gamma^T * Fneighbor
    for (unsigned j = 0; j < 9; ++j) {
      for (unsigned i = 0; i < Fneighbor.rows(); ++i) {
        Fneighbor(i,j) *= Gamma[i];
      }
    }
  }
}

void VectorView::set_to_product(const MatrixView& m, const VectorView& v,
                                const bool transpose)
{
  CBLAS_TRANSPOSE tr;
  if (transpose){
    tr = CblasTrans;
    assert(m.cols() == length());
    assert(m.rows() == v.length());
  } else {
    tr = CblasNoTrans;
    assert(m.cols() == v.length());
    assert(m.rows() == length());
  }

  cblas_dgemv(CblasColMajor, tr, m.rows(), m.cols(), 1.0, m.data(),
              m.stride(), v.data(), 1, 0.0, data_, 1);
}

void seissol::model::getTransposedElasticCoefficientMatrix( seissol::model::ElasticMaterial const& i_material,
                                                            unsigned i_dim,
                                                            MatrixView<9, 9> o_M )
{
  o_M.setZero();

  real lambda2mu = i_material.lambda + 2.0 * i_material.mu;
  real rhoInv = 1.0 / i_material.rho;

  switch (i_dim)
  {
    case 0:
      o_M(6,0) = -lambda2mu;
      o_M(6,1) = -i_material.lambda;
      o_M(6,2) = -i_material.lambda;
      o_M(7,3) = -i_material.mu;
      o_M(8,5) = -i_material.mu;
      o_M(0,6) = -rhoInv;
      o_M(3,7) = -rhoInv;
      o_M(5,8) = -rhoInv;
      break;

    case 1:
      o_M(7,0) = -i_material.lambda;
      o_M(7,1) = -lambda2mu;
      o_M(7,2) = -i_material.lambda;
      o_M(6,3) = -i_material.mu;
      o_M(8,4) = -i_material.mu;
      o_M(3,6) = -rhoInv;
      o_M(1,7) = -rhoInv;
      o_M(4,8) = -rhoInv;
      break;

    case 2:
      o_M(8,0) = -i_material.lambda;
      o_M(8,1) = -i_material.lambda;
      o_M(8,2) = -lambda2mu;
      o_M(7,4) = -i_material.mu;
      o_M(6,5) = -i_material.mu;
      o_M(5,6) = -rhoInv;
      o_M(4,7) = -rhoInv;
      o_M(2,8) = -rhoInv;
      break;
      
    default:
      break;
  }
}

    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);
        }
    }

void seissol::model::getTransposedElasticGodunovState( seissol::model::ElasticMaterial const& local,
                                                       seissol::model::ElasticMaterial const& neighbor,
                                                       MatrixView<9, 9> QgodLocal,
                                                       MatrixView<9, 9> QgodNeighbor )
{
  QgodNeighbor.setZero();
  
  real cpL = sqrt((local.lambda + 2.0 * local.mu)       / local.rho);
  real cpN = sqrt((neighbor.lambda + 2.0 * neighbor.mu) / neighbor.rho);
  real csL = sqrt(local.mu / local.rho);
  real csN = sqrt(neighbor.mu / neighbor.rho);
  
  real constP = cpN * (local.lambda + 2.0 * local.mu) + cpL * (neighbor.lambda + 2.0 * neighbor.mu);
  real constS = csN * local.mu + csL * neighbor.mu;
  
  QgodNeighbor(0,0) = cpN * (local.lambda + 2.0 * local.mu) / constP;
  QgodNeighbor(6,0) = (local.lambda + 2.0 * local.mu) * (neighbor.lambda + 2.0 * neighbor.mu) / constP;
  QgodNeighbor(0,1) = cpN * local.lambda / constP;
  QgodNeighbor(6,1) = local.lambda * (neighbor.lambda + 2.0 * neighbor.mu) / constP;
  QgodNeighbor(0,2) = QgodNeighbor(0,1);
  QgodNeighbor(6,2) = QgodNeighbor(6,1);
  QgodNeighbor(3,3) = csN * local.mu / constS;
  QgodNeighbor(7,3) = local.mu * neighbor.mu / constS;
  QgodNeighbor(5,5) = QgodNeighbor(3,3);
  QgodNeighbor(8,5) = QgodNeighbor(7,3);
  QgodNeighbor(0,6) = cpL * cpN / constP;
  QgodNeighbor(6,6) = cpL * (neighbor.lambda + 2.0 * neighbor.mu) / constP;
  QgodNeighbor(3,7) = csL * csN / constS;
  QgodNeighbor(7,7) = csL * neighbor.mu / constS;
  QgodNeighbor(5,8) = QgodNeighbor(3,7);
  QgodNeighbor(8,8) = QgodNeighbor(7,7);

  // QgodLocal = I - QgodNeighbor
  for (unsigned idx = 0; idx < QgodLocal.rows() * QgodLocal.cols(); ++idx) {
    QgodLocal.data[idx] = -QgodNeighbor.data[idx];
  }  
  for (unsigned idx = 0; idx < QgodLocal.rows() && idx < QgodLocal.cols(); ++idx) {
    QgodLocal(idx, idx) += 1.0;
  }
}

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