C++ Transpose函数代码示例

void FoxLi( Matrix<Complex<Real>>& A, Int n, Real omega )
{
    DEBUG_CSE
    typedef Complex<Real> C;
    const Real pi = 4*Atan( Real(1) );
    const C phi = Sqrt( C(0,omega/pi) ); 
    
    // Compute Gauss quadrature points and weights
    Matrix<Real> d, e; 
    Zeros( d, n, 1 );
    e.Resize( n-1, 1 );
    for( Int j=0; j<n-1; ++j )
    {
        const Real betaInv = 2*Sqrt(1-Pow(j+Real(1),-2)/4);
        e(j) = 1/betaInv;
    }
    Matrix<Real> x, Z;
    HermitianTridiagEig( d, e, x, Z, UNSORTED );
    auto z = Z( IR(0), ALL );
    Matrix<Real> sqrtWeights( z ), sqrtWeightsTrans;
    for( Int j=0; j<n; ++j )
        sqrtWeights(0,j) = Sqrt(Real(2))*Abs(sqrtWeights(0,j));
    herm_eig::Sort( x, sqrtWeights, ASCENDING );
    Transpose( sqrtWeights, sqrtWeightsTrans );

    // Form the integral operator
    A.Resize( n, n );
    for( Int j=0; j<n; ++j )
    {
        for( Int i=0; i<n; ++i )
        {
            const Real theta = -omega*Pow(x(i)-x(j),2);
            const Real realPart = Cos(theta);
            const Real imagPart = Sin(theta);
            A(i,j) = phi*C(realPart,imagPart);
        }
    }

    // Apply the weighting
    DiagonalScale( LEFT, NORMAL, sqrtWeightsTrans, A );
    DiagonalScale( RIGHT, NORMAL, sqrtWeightsTrans, A );
}

// Rotation kring godtycklig axel (enbart rotationen)
void ArbRotate(Point3D *axis, GLfloat fi, GLfloat *m)
{
	Point3D x, y, z, a;
	GLfloat R[16], Rt[16], Raxel[16], RtRx[16];
	
// Kolla ocksΠom parallell med Z-axel!
	if (axis->x < 0.0000001) // Under nŒgon tillrŠckligt liten grŠns
	if (axis->x > -0.0000001)
	if (axis->y < 0.0000001)
	if (axis->y > -0.0000001)
		if (axis->z > 0)
		{
			Rz(fi, m);
			return;
		}
		else
		{
			Rz(-fi, m);
			return;
		}

	x = *axis;
	Normalize(&x); // |x|
	SetVector(0,0,1, &z); // Temp z
	CrossProduct(&z, &x, &y);
	Normalize(&y); // y' = z^ x x'
	CrossProduct(&x, &y, &z); // z' = x x y

	R[0] = x.x; R[4] = x.y; R[8] = x.z;  R[12] = 0.0;
	R[1] = y.x; R[5] = y.y; R[9] = y.z;  R[13] = 0.0;
	R[2] = z.x; R[6] = z.y; R[10] = z.z;  R[14] = 0.0;

	R[3] = 0.0; R[7] = 0.0; R[11] = 0.0;  R[15] = 1.0;

	Transpose(&R, &Rt); // Transpose = Invert -> felet ej i Transpose, och det Šr en ortonormal matris
	
	Rx(fi, &Raxel); // Rotate around x axis
	
	// m := Rt * Rx * R
	Mult(&Rt, &Raxel, &RtRx);
	Mult(&RtRx, &R, m);
}

void SUMMA_NTA
( Orientation orientB,
  T alpha,
  const AbstractDistMatrix<T>& APre,
  const AbstractDistMatrix<T>& BPre,
        AbstractDistMatrix<T>& CPre )
{
    EL_DEBUG_CSE
    const Int n = CPre.Width();
    const Int bsize = Blocksize();
    const Grid& g = APre.Grid();
    const bool conjugate = ( orientB == ADJOINT );

    DistMatrixReadProxy<T,T,MC,MR> AProx( APre );
    DistMatrixReadProxy<T,T,MC,MR> BProx( BPre );
    DistMatrixReadWriteProxy<T,T,MC,MR> CProx( CPre );
    auto& A = AProx.GetLocked();
    auto& B = BProx.GetLocked();
    auto& C = CProx.Get();

    // Temporary distributions
    DistMatrix<T,MR,STAR> B1Trans_MR_STAR(g);
    DistMatrix<T,MC,STAR> D1_MC_STAR(g);

    B1Trans_MR_STAR.AlignWith( A );
    D1_MC_STAR.AlignWith( A );

    for( Int k=0; k<n; k+=bsize )
    {
        const Int nb = Min(bsize,n-k);
        auto B1 = B( IR(k,k+nb), ALL        );
        auto C1 = C( ALL,        IR(k,k+nb) );

        // C1[MC,*] := alpha A[MC,MR] (B1^[T/H])[MR,*]
        Transpose( B1, B1Trans_MR_STAR, conjugate );
        LocalGemm( NORMAL, NORMAL, alpha, A, B1Trans_MR_STAR, D1_MC_STAR );

        // C1[MC,MR] += scattered result of D1[MC,*] summed over grid rows
        AxpyContract( T(1), D1_MC_STAR, C1 );
    }
}

// AnimatedTransform Method Definitions
void AnimatedTransform::Decompose(const Matrix4x4 &m, Vector *T,
                                  Quaternion *Rquat, Matrix4x4 *S) {
    // Extract translation _T_ from transformation matrix
    T->x = m.m[0][3];
    T->y = m.m[1][3];
    T->z = m.m[2][3];

    // Compute new transformation matrix _M_ without translation
    Matrix4x4 M = m;
    for (int i = 0; i < 3; ++i)
        M.m[i][3] = M.m[3][i] = 0.f;
    M.m[3][3] = 1.f;

    // Extract rotation _R_ from transformation matrix
    float norm;
    int count = 0;
    Matrix4x4 R = M;
    do {
        // Compute next matrix _Rnext_ in series
        Matrix4x4 Rnext;
        Matrix4x4 Rit = Inverse(Transpose(R));
        for (int i = 0; i < 4; ++i)
            for (int j = 0; j < 4; ++j)
                Rnext.m[i][j] = 0.5f * (R.m[i][j] + Rit.m[i][j]);

        // Compute norm of difference between _R_ and _Rnext_
        norm = 0.f;
        for (int i = 0; i < 3; ++i) {
            float n = fabsf(R.m[i][0] - Rnext.m[i][0]) +
                      fabsf(R.m[i][1] - Rnext.m[i][1]) +
                      fabsf(R.m[i][2] - Rnext.m[i][2]);
            norm = max(norm, n);
        }
        R = Rnext;
    } while (++count < 100 && norm > .0001f);
    // XXX TODO FIXME deal with flip...
    *Rquat = Quaternion(R);

    // Compute scale _S_ using rotation and original matrix
    *S = Matrix4x4::Mul(Inverse(R), M);
}

int CalcSphereCenter (const Point<3> ** pts, Point<3> & c)
{
  Vec3d row1 (*pts[0], *pts[1]);
  Vec3d row2 (*pts[0], *pts[2]);
  Vec3d row3 (*pts[0], *pts[3]);

  Vec3d rhs(0.5 * (row1*row1),
	    0.5 * (row2*row2),
	    0.5 * (row3*row3));
  Transpose (row1, row2, row3);
  
  Vec3d sol;
  if (SolveLinearSystem (row1, row2, row3, rhs, sol))
    {
      (*testout) << "CalcSphereCenter: degenerated" << endl;
      return 1;
    }

  c = *pts[0] + sol;
  return 0;
}

Matrix Matrix::Inverse() const
{
	if (Determinant() == 0)
	{
		return Matrix(
				Vector3D(0, 0, 0),
				Vector3D(0, 0, 0),
				Vector3D(0, 0, 0));
	}
	else
	{
		Matrix matrix = Transpose();
		Float blah = Determinant();
		Float invDet = 1.0/Determinant();

		return Matrix(
				matrix.Rows[1].CrossProduct(matrix.Rows[2]) * invDet,
				matrix.Rows[2].CrossProduct(matrix.Rows[0]) * invDet,
				matrix.Rows[0].CrossProduct(matrix.Rows[1]) * invDet);
	}
}

void LapackInvAndDet(cDMatrix& theMatrix, cDMatrix& theInvMatrix, double& theDet)
{
uint myNCol = theMatrix.GetNCols() ;

double  *myAP = new double[myNCol*(myNCol + 1)/2],
                *myW = new double[myNCol],
                *myZ = new double[myNCol*myNCol],
                *myWork = new double[myNCol * 3] ;
int myInfo,
        myN = (int)(myNCol),
        myldz = (int)(myNCol) ;

        for (register int i = 0 ; i < myN ; i++)
                for (register int j = i ; j < myldz ; j++)
                        myAP[i+(j+1)*j/2]  = theMatrix[i][j] ;

        F77_NAME(dspev)("V", "U", &myN, myAP, myW, myZ, &myldz, myWork, &myInfo) ;

        if (myInfo != 0)
                throw cOTError("Non inversible matrix") ;
        theDet = 1.0L ;
cDVector myInvEigenValue = cDVector(myNCol) ;

cDMatrix myEigenVector(myNCol, myNCol) ;
        for (register uint i = 0 ; i < myNCol ; i++)
        {       theDet *= myW[i] ;
                myInvEigenValue[i] = 1.0 /myW[i] ;
                for (register int j = 0 ; j < myN ; j++)
                        myEigenVector[i][j] = myZ[i + j*myN] ;
        }
        theInvMatrix =  myEigenVector ;
cDMatrix myAuxMat1 = Diag(myInvEigenValue), myAuxMat2 = Transpose(myEigenVector) ;
cDMatrix myAuxMat = myAuxMat1 * myAuxMat2 ;
        theInvMatrix = theInvMatrix * myAuxMat ;
        
        delete myAP ;
        delete myW ;
        delete myZ ;
        delete myWork ;
}

void Camera::Precompute()
{
    //*********************************************************
    //Compute m_mKRt
    //*********************************************************
    KRt_ = K_ * Transpose(R_);

    //*********************************************************
    //Compute KRtT
    //*********************************************************
    KRtT_ = KRt_ * t_;

    //*********************************************************
    //Compute VPN
    //*********************************************************

    VPN_ = R_.getColumn(2);

    //*********************************************************
    //Compute VPd
    //*********************************************************

    VPd_ = 0.0;

    for(int i = 0; i < t_.getSize(); i++)
    {
        VPd_ += t_(i)*(VPN_(i)*(-1.0));
    }

    //*********************************************************
    //Compute GPN and GPD
    //*********************************************************
    GPD_ = GP_(3);
    GPN_.setSize(3);

    for(int i = 0; i < GPN_.getSize(); i++)
    {
        GPN_(i) = GP_(i);
    }
}

void SUMMA_TNC
( Orientation orientA,
  T alpha,
  const AbstractDistMatrix<T>& APre,
  const AbstractDistMatrix<T>& BPre,
        AbstractDistMatrix<T>& CPre )
{
    DEBUG_CSE
    const Int sumDim = BPre.Height();
    const Int bsize = Blocksize();
    const Grid& g = APre.Grid();

    DistMatrixReadProxy<T,T,MC,MR> AProx( APre );
    DistMatrixReadProxy<T,T,MC,MR> BProx( BPre );
    DistMatrixReadWriteProxy<T,T,MC,MR> CProx( CPre );
    auto& A = AProx.GetLocked();
    auto& B = BProx.GetLocked();
    auto& C = CProx.Get();

    // Temporary distributions
    DistMatrix<T,STAR,MC> A1_STAR_MC(g);
    DistMatrix<T,MR,STAR> B1Trans_MR_STAR(g);

    A1_STAR_MC.AlignWith( C );
    B1Trans_MR_STAR.AlignWith( C );

    for( Int k=0; k<sumDim; k+=bsize )
    {
        const Int nb = Min(bsize,sumDim-k);
        auto A1 = A( IR(k,k+nb), ALL );
        auto B1 = B( IR(k,k+nb), ALL );

        // C[MC,MR] += alpha (A1[*,MC])^T B1[*,MR]
        //           = alpha (A1^T)[MC,*] B1[*,MR]
        A1_STAR_MC = A1; 
        Transpose( B1, B1Trans_MR_STAR );
        LocalGemm
        ( orientA, TRANSPOSE, alpha, A1_STAR_MC, B1Trans_MR_STAR, T(1), C );
    }
}

Vector<double> AncillaryMethods::PlaneToCam(const Camera& camera)
{
    Vector<double> plane = camera.get_GP();

    Vector<double> pv(plane(0), plane(1), plane(2));

    Matrix<double> cam_rot_trans = Transpose(camera.get_R());

    pv = cam_rot_trans * pv;

    Vector<double> t = cam_rot_trans * camera.get_t();

    double d = plane(3) + pv(0)*t(0) + pv(1)*t(1) + pv(2)*t(2);

    Vector<double> gp_in_camera(4);
    gp_in_camera(0) = pv(0);
    gp_in_camera(1) = pv(1);
    gp_in_camera(2) = pv(2);
    gp_in_camera(3) = d;

    return gp_in_camera;
}

int main()
{
  double Av[9] = {2.0, 1.0, 1.0,
                  1.0, 2.0, 1.0,
                  1.0, 1.0, 2.0, };

  Matrix<double> A(3, 3, Av);
  Matrix<double> B = 0.5 * Transpose(A) * A;
  printf("Matrix:\n");
  B.Print();
  char outputFilename[96] = "BMatrix";
  B.Save(outputFilename);

  Matrix<double> Q(3,3); // will hold eigenvectors
  Matrix<double> Lambda(3,1); // will hold eigenvalues
  B.SymmetricEigenDecomposition(Q, Lambda);

  printf("Eigenvalues:\n");
  Lambda.Print();

  return 0;
}

 PhaseEnumerationCache
 ( const Matrix<Field>& B,
   const Matrix<Base<Field>>& d,
   const Matrix<Field>& N,
   const Matrix<Base<Field>>& normUpperBounds,
         Int batchSize=256,
         Int blocksize=32,
         bool useTranspose=true )
 : B_(B),
   d_(d),
   N_(N),
   normUpperBounds_(normUpperBounds),
   foundVector_(false),
   numQueued_(0),
   insertionBound_(normUpperBounds.Height()),
   blocksize_(blocksize),
   useTranspose_(useTranspose)
 { 
     Zeros( Y_, N.Height(), batchSize );   
     if( useTranspose )
         Transpose( N, NTrans_ );
 }

void SBVAR_symmetric::SetupSBVAR_symmetric(void)
{
  prior_YY=TransposeMultiply(prior_Y,prior_Y);
  prior_XX=TransposeMultiply(prior_X,prior_X);
  prior_XY=TransposeMultiply(prior_X,prior_Y);

  if (flat_prior)
    log_prior_constant=0.0;
  else
    {
      TDenseMatrix S(n_vars+n_predetermined,n_vars+n_predetermined);
      S.Insert(0,0,prior_YY);
      S.Insert(n_vars,0,-prior_XY);
      S.Insert(0,n_vars,-Transpose(prior_XY));
      S.Insert(n_vars,n_vars,prior_XX);
      log_prior_constant=n_vars*(-0.918938533204673*(n_vars+n_predetermined) + 0.5*LogAbsDeterminant(S));  // 0.918938533204673 = 0.5*ln(2*pi)
    }

  prior_YY*=lambda_bar;
  prior_XX*=lambda_bar;
  prior_XY*=lambda_bar;
  log_prior_constant*=lambda_bar;
}

void CallFunction(FunctionCall fc)
{
	int function = GetFunction(fc->function);
	
	switch (function)
	{
		case VAR : NewVariable(fc); break;
		case NMX : NewMatrix(fc); break;
		case ADD : Addition(fc); break;
		case SUB : Substraction(fc); break;
		case MUL : Multiplication(fc); break;
		case MSC : Scalar_Mult(fc); break;
		case EXP : Exponentiation(fc); break;
		case TRA : Transpose(fc); break;
		case DET : Determinant(fc); break;
		case DLU : Decomposition(fc); break;
		case SOL : Solve(fc); break;
		case INV : Inversion(fc); break;
		case RNK : Rank(fc); break;
		case DSP : Display(fc); break;
		case NOF : // default
		default :
		{
			if (GetFunction(fc->name)==SPT) SpeedTest(fc);
			else if (IndexVariable(fc->function)!=-1) NewVariable(fc);
			else if (IndexMatrix(fc->function)!=-1) NewMatrix(fc);
			else
			{
				printf("\t%s : Function Not Implemented\n", fc->function);
				fni++;
			}
			break;
		}
	}
	
	if (function!=NOF && function !=VAR) fni = 0;
}

    //-----------------------------------------------------------------------------
    //  Update
    //  Updates the object
    //  TODO: Pre- and Post- updates?
    //-----------------------------------------------------------------------------
    void CView::Update( void )
    {
        RVector4 x, y, z;

        z = m_vLook = Normalize( m_vLook );
        x = m_vRight = Normalize( CrossProduct( m_vUp, z ) );
        y = CrossProduct( z, x );
        
        m_mView.r0 = x;
        m_mView.r1 = y;
        m_mView.r2 = z;
        m_mView.r3 = RVector4Zero();
        m_mView = Transpose( m_mView );
        
        m_mView.r3 = RVector4( -DotProduct( x, m_vPosition), -DotProduct( y, m_vPosition), -DotProduct( z, m_vPosition), 1.0f );

        
        //z = Normalize( RQuatGetZAxis(m_Transform.orientation) );
        //x = Normalize( CrossProduct( RVector3(0.0f,1.0f,0.0f), z ) );
        //y = CrossProduct( z, x );
        //
        //m_mView.r0 = Homogonize( x );
        //m_mView.r1 = Homogonize( y );
        //m_mView.r2 = Homogonize( z );
        //m_mView.r3 = RVector4Zero();
        //m_mView = Transpose( m_mView );
        //
        //m_mView.r3 = RVector4( -DotProduct( x, m_Transform.position), -DotProduct( y, m_Transform.position), -DotProduct( z, m_Transform.position), 1.0f );
    
        char szCameraData[256] = { 0 };

        sprintf( szCameraData, "Pos:  %f, %f, %f", m_vPosition.x, m_vPosition.y, m_vPosition.z );
        Engine::GetRenderer()->DrawString( 200, 16, szCameraData );
        sprintf( szCameraData, "Look: %f, %f, %f", m_vLook.x, m_vLook.y, m_vLook.z );
        Engine::GetRenderer()->DrawString( 200, 32, szCameraData );
    }