C++ Matrix::Get方法代码示例

void PushPairInto
(       Matrix<Real>& s, 
        Matrix<Real>& z,
  const Matrix<Real>& w,
  const Matrix<Int>& orders, 
  const Matrix<Int>& firstInds,
  Real wMaxNormLimit )
{
    DEBUG_ONLY(CSE cse("soc::PushPairInto"))

    Matrix<Real> sLower, zLower;
    soc::LowerNorms( s, sLower, orders, firstInds );
    soc::LowerNorms( z, zLower, orders, firstInds );

    const Int height = s.Height();
    for( Int i=0; i<height; ++i )
    {
        const Real w0 = w.Get(i,0);
        const Int firstInd = firstInds.Get(i,0);
        if( i == firstInd && w0 > wMaxNormLimit )
        {
            // TODO: Switch to a non-adhoc modification     
            s.Update( i, 0, Real(1)/wMaxNormLimit );
            z.Update( i, 0, Real(1)/wMaxNormLimit );
        }
    }
}

double gen_adjsrc
( const Matrix<Complex<double>>& S, 
  const Matrix<Complex<double>>& O, 
        Matrix<Complex<double>>& A )
{
    // TODO: Compute local misfits and then sum the result with an
    //       MPI_Allreduce summation
    // NOTE: This routine should already work but may not be as fast
    //       as possible.
    double l2MisfitSquared=0.0; 
    const int nsrc = S.Height();
    const int nrec = S.Width();

    A.Resize( nsrc, nrec );
    for (int isrc=0;isrc<nsrc;isrc++ ) { 
        for (int irec=0;irec<nrec;irec++ ) { 

	        Complex<double> obs = O.Get(irec,isrc);
	        Complex<double> syn = S.Get(irec,isrc);
	        A.Set( irec, isrc, obs-syn );

            l2MisfitSquared += Abs(A.Get(irec,isrc));
        }
    }
    return l2MisfitSquared;
}

void SOCSquareRoot
( const Matrix<Real>& x, 
        Matrix<Real>& xRoot,
  const Matrix<Int>& orders, 
  const Matrix<Int>& firstInds )
{
    DEBUG_ONLY(CSE cse("SOCSquareRoot"))

    Matrix<Real> d;
    SOCDets( x, d, orders, firstInds );
    ConeBroadcast( d, orders, firstInds );

    const Int height = x.Height();
    Zeros( xRoot, height, 1 );
    for( Int i=0; i<height; )
    {
        const Int order = orders.Get(i,0);
        const Int firstInd = firstInds.Get(i,0);
        if( i != firstInd )       
            LogicError("Inconsistency in orders and firstInds");

        const Real eta0 = Sqrt(x.Get(i,0)+Sqrt(d.Get(i,0)))/Sqrt(Real(2));
        xRoot.Set( i, 0, eta0 );
        for( Int k=1; k<order; ++k )
            xRoot.Set( i+k, 0, x.Get(i+k,0)/(2*eta0) );

        i += order;
    }
}

Matrix operator-(Matrix a,Matrix b){
	Matrix result(a.Rows(),a.Cols());
	for(int i = 0 ; i < a.Rows() ; i++)
		for( int j = 0 ; j < a.Cols() ; j++)
			result.Set(j,i,a.Get(j,i) - b.Get(j,i));
	return result;
}

void
QuasiDiagonalSolve
( LeftOrRight side, UpperOrLower uplo, 
  const Matrix<FMain>& d,
  const Matrix<F>& dSub, 
        Matrix<F>& X,
  bool conjugated )
{
    DEBUG_CSE
    const Int m = X.Height();
    const Int n = X.Width();

    Matrix<F> D( 2, 2 );
    if( side == LEFT && uplo == LOWER )
    {
        if( m == 0 )
            return;
        DEBUG_ONLY(
          if( d.Height() != m )
              LogicError
              ("d was the wrong size: m=",m,",n=",n,", d ~ ",
               d.Height()," x ",d.Width());
          if( dSub.Height() != m-1 )
              LogicError
              ("dSub was the wrong size: m=",m,",n=",n,", dSub ~ ",
               dSub.Height()," x ",dSub.Width());
        )
        Int i=0;
        while( i < m )
        {
            Int nb;
            if( i < m-1 && Abs(dSub.Get(i,0)) > 0 )
                nb = 2;
            else
                nb = 1;

            auto XRow = X( IR(i,i+nb), IR(0,n) );
            if( nb == 1 )
            {
                Scale( F(1)/d.Get(i,0), XRow );
            }
            else
            {
                D.Set(0,0,d.Get(i,0));
                D.Set(1,1,d.Get(i+1,0));
                D.Set(1,0,dSub.Get(i,0));
                Symmetric2x2Inv( LOWER, D, conjugated );
                MakeSymmetric( LOWER, D, conjugated );

                Transform2x2Rows( D, X, i, i+1 );
            }

            i += nb;
        }
    }

inline Complex<R>
Reflector( Matrix<Complex<R> >& chi, Matrix<Complex<R> >& x )
{
#ifndef RELEASE
    CallStackEntry entry("Reflector");
#endif
    typedef Complex<R> C;

    R norm = Nrm2( x );
    C alpha = chi.Get(0,0);

    if( norm == 0 && alpha.imag == R(0) )
    {
        chi.Set(0,0,-chi.Get(0,0));
        return C(2);
    }

    R beta;
    if( alpha.real <= 0 )
        beta = lapack::SafeNorm( alpha.real, alpha.imag, norm );
    else
        beta = -lapack::SafeNorm( alpha.real, alpha.imag, norm );

    const R one = 1;
    const R safeMin = lapack::MachineSafeMin<R>();
    const R epsilon = lapack::MachineEpsilon<R>();
    const R safeInv = safeMin/epsilon;
    Int count = 0;
    if( Abs(beta) < safeInv )
    {
        R invOfSafeInv = one/safeInv;
        do
        {
            ++count;
            Scale( invOfSafeInv, x );
            alpha *= invOfSafeInv;
            beta *= invOfSafeInv;
        } while( Abs(beta) < safeInv );

        norm = Nrm2( x );
        if( alpha.real <= 0 )
            beta = lapack::SafeNorm( alpha.real, alpha.imag, norm );
        else
            beta = -lapack::SafeNorm( alpha.real, alpha.imag, norm );
    }

    C tau = C( (beta-alpha.real)/beta, -alpha.imag/beta );
    Scale( one/(alpha-beta), x );

    for( Int j=0; j<count; ++j )
        beta *= safeInv;
    chi.Set(0,0,beta);
    return tau;
}

void SOCMaxEig
( const Matrix<Real>& x, 
        Matrix<Real>& maxEigs,
  const Matrix<Int>& orders, 
  const Matrix<Int>& firstInds )
{
    DEBUG_ONLY(CSE cse("SOCMaxEig"))
    SOCLowerNorms( x, maxEigs, orders, firstInds );
    const Int height = x.Height();
    for( Int i=0; i<height; ++i )
        if( i == firstInds.Get(i,0) ) 
            maxEigs.Set( i, 0, x.Get(i,0)+maxEigs.Get(i,0) );
}

inline void
SortEig( Matrix<R>& w, Matrix<R>& Z )
{
#ifndef RELEASE
    PushCallStack("SortEig");
#endif
    const int n = Z.Height();
    const int k = Z.Width();

    // Initialize the pairs of indices and eigenvalues
    std::vector<internal::IndexValuePair<R> > pairs( k );
    for( int i=0; i<k; ++i )
    {
        pairs[i].index = i;
        pairs[i].value = w.Get(i,0);
    }

    // Sort the eigenvalues and simultaneously form the permutation
    std::sort
    ( pairs.begin(), pairs.end(), internal::IndexValuePair<R>::Compare );

    // Reorder the eigenvectors and eigenvalues using the new ordering
    Matrix<R> ZPerm( n, k );
    for( int j=0; j<k; ++j )
    {
        const int source = pairs[j].index;
        MemCopy( ZPerm.Buffer(0,j), Z.LockedBuffer(0,source), n );
        w.Set(j,0,pairs[j].value);
    }
    Z = ZPerm;
#ifndef RELEASE
    PopCallStack();
#endif
}

inline SafeProduct<F> 
AfterLUPartialPiv( const Matrix<F>& A, const Matrix<Int>& p )
{
#ifndef RELEASE
    CallStackEntry entry("determinant::AfterLUPartialPiv");
#endif
    if( A.Height() != A.Width() )
        LogicError("Cannot compute determinant of nonsquare matrix");
    if( A.Height() != p.Height() )
        LogicError("Pivot vector is incorrect length");

    typedef BASE(F) R;
    const Int n = A.Height();

    Matrix<F> d;
    A.GetDiagonal( d );
    const R scale(n);
    SafeProduct<F> det( n );
    for( Int i=0; i<n; ++i )
    {
        const F delta = d.Get(i,0);
        R alpha = Abs(delta);
        det.rho *= delta/alpha;
        det.kappa += Log(alpha)/scale;
    }
    const bool isOdd = PivotParity( p );
    if( isOdd )
        det.rho = -det.rho;

    return det;
}

FrobeniusNorm( const Matrix<F>& A )
{
#ifndef RELEASE
    CallStackEntry entry("FrobeniusNorm");
#endif
    typedef BASE(F) R;
    R scale = 0;
    R scaledSquare = 1;
    const Int width = A.Width();
    const Int height = A.Height();
    for( Int j=0; j<width; ++j )
    {
        for( Int i=0; i<height; ++i )
        {
            const R alphaAbs = Abs(A.Get(i,j));
            if( alphaAbs != 0 )
            {
                if( alphaAbs <= scale )
                {
                    const R relScale = alphaAbs/scale;
                    scaledSquare += relScale*relScale;
                }
                else
                {
                    const R relScale = scale/alphaAbs;
                    scaledSquare = scaledSquare*relScale*relScale + 1;
                    scale = alphaAbs;
                }
            }
        }
    }
    return scale*Sqrt(scaledSquare);
}

double Matrix::dot(Matrix other) {
	double dot = 0;
	for(int i = 0 ; i < 3 ; i ++ ){
		dot = dot + (Get(i,0) * other.Get(i,0));
	}
	return dot;
}

void Transform2x2( const Matrix<T>& G, Matrix<T>& a1, Matrix<T>& a2 )
{
    DEBUG_CSE
    T* a1Buf = a1.Buffer();
    T* a2Buf = a2.Buffer();
    const Int inc1 = ( a1.Height() == 1 ? a1.LDim() : 1 );
    const Int inc2 = ( a2.Height() == 1 ? a2.LDim() : 1 );
    const Int n = ( a1.Height() == 1 ? a1.Width() : a1.Height() );

    const T gamma11 = G.Get(0,0);
    const T gamma12 = G.Get(0,1);
    const T gamma21 = G.Get(1,0);
    const T gamma22 = G.Get(1,1);
    Transform2x2
    ( n, gamma11, gamma12, gamma21, gamma22, a1Buf, inc1, a2Buf, inc2 );
}

Matrix::Matrix(Matrix& x) {
	int i,j;
	AllocMxMem(x.GetZSize(),x.GetSSize());
	for (i=1;i<=zSize;i++)
		for(j=1;j<=sSize;j++)
			Set(i,j,x.Get(i,j));
}

void LLN( const Matrix<F>& L, Matrix<F>& X, bool checkIfSingular )
{
    DEBUG_CSE
    const Int m = X.Height();
    const Int bsize = Blocksize();

    for( Int k=0; k<m; k+=bsize )
    {
        const Int nbProp = Min(bsize,m-k);
        const bool in2x2 = ( k+nbProp<m && L.Get(k+nbProp-1,k+nbProp) != F(0) );
        const Int nb = ( in2x2 ? nbProp+1 : nbProp );

        const Range<Int> ind1( k,    k+nb ),
                         ind2( k+nb, m    );

        auto L11 = L( ind1, ind1 );
        auto L21 = L( ind2, ind1 );

        auto X1 = X( ind1, ALL );
        auto X2 = X( ind2, ALL );

        LLNUnb( L11, X1, checkIfSingular );
        Gemm( NORMAL, NORMAL, F(-1), L21, X1, F(1), X2 );
    }
}

		// ==================================================================================================================
		int mainloop() {
			this->render.BeginScene();
			{		
				Matrix worldMatrix;
				worldMatrix.Identity();
				worldMatrix.RotateRPY(t, .3455*t, .7346*t);

				struct {
					matrix worldMatrix;
					matrix viewMatrix;
					matrix projectionMatrix;
				} viewMatrices;

				m_camera->Calculate(render);

				viewMatrices.worldMatrix = XMMatrixTranspose(worldMatrix.Get());
				viewMatrices.viewMatrix = XMMatrixTranspose(m_camera->GetViewMatrix().Get());
				viewMatrices.projectionMatrix = XMMatrixTranspose(m_camera->GetProjectionMatrix().Get());

				m_vertexShader->GetParam("MatrixBuffer").SetP(&viewMatrices);

				m_fragmentShader->GetBRes("diffuse").Set(m_texture->GetTextureResource());
				m_fragmentShader->GetBRes("SampleType").Set(m_textureSampler->GetSamplerState());

				m_vertexShader->Render(this->render);
				m_fragmentShader->Render(this->render);

				m_model->RenderMesh(this->render);

				this->t += 0.01;
			}

			this->render.EndScene();
			return 0;
		};