C++ B函数代码示例

 // Don't crash here.
 void f() {
   X x = X();
   (void)noexcept(B(declval<B>()));
 }

bool	Field::initializeBuffers(ID3D11Device* device)
{
	Dictionnary*		dico;
	VertexType*			vertices;
	Buffers*			buffer;
	unsigned long*		indices;
	int					index;	
	int					vertexCount;
	int					indexCount;

	float				fMax = fSpacing * (unSize/2);
	float				fMin = -fMax;
	float				fDiff;
	bool				filled;
	filled = true;
	fDiff = fMax - fMin;
	Generate();
	
	vertexCount = (unSize-1) * (unSize-1) * 20;
	indexCount = vertexCount;

	if ((dico = Dictionnary::getInstance()) == nullptr)
		return false;
	vertices = new VertexType[vertexCount];
	if (vertices == false)
		return (false);
	indices = new unsigned long[indexCount];
	if (indices == false)
	{
		delete []vertices;
		return (false);
	}
	buffer = new Buffers();
	if (buffer == nullptr)
	{
		delete []vertices;
		delete []indices;
		return (false);
	}

	index = 0;
	DirectX::XMVECTOR	vec1;
	DirectX::XMVECTOR	vec2;
	DirectX::XMVECTOR	normal;

	for(unsigned int i = 0; i < (unSize - 1); i++)
	{
		float	fZ = fMin + i * fSpacing; 
		for(unsigned int j = 0; j < (unSize - 1); j++)
		{
			float fX = fMin + j * fSpacing;
			if(filled)
			{
				DirectX::XMFLOAT3	A(fX, vectPoints[i][j], fZ);
				DirectX::XMFLOAT3	B(fX, vectPoints[i + 1][j], fZ + fSpacing);
				DirectX::XMFLOAT3	C(fX + fSpacing, vectPoints[i + 1][j + 1], fZ + fSpacing);
				DirectX::XMFLOAT3	D(fX + fSpacing, vectPoints[i][j + 1], fZ);

				vec1 = DirectX::XMLoadFloat3(&DirectX::XMFLOAT3(B.x - A.x, B.y - A.y, B.z - A.z));
				vec2 = DirectX::XMLoadFloat3(&DirectX::XMFLOAT3(C.x - A.x, C.y - A.y, C.z - A.z));
				normal = DirectX::XMVector3Cross(vec1, vec2);
			setVertices(fX, vectPoints[i][j], fZ, i, j, vertices, indices, &index, normal, (fX + fMax) / fMax * 2, (fZ + fMax) / fMax * 2);
			setVertices(fX, vectPoints[i + 1][j], fZ + fSpacing, i, j, vertices, indices, &index, normal, (fX + fMax) / fMax * 2, (fZ + fSpacing + fMax) / fMax * 2);
			setVertices(fX + fSpacing, vectPoints[i + 1][j + 1], fZ + fSpacing, i, j, vertices, indices, &index, normal, (fX + fSpacing + fMax) / fMax * 2, (fZ + fSpacing + fMax) / fMax * 2);

				/*vec1 = DirectX::XMLoadFloat3(&DirectX::XMFLOAT3(C.x - B.x, C.y - B.y, C.z - B.z));
				vec2 = DirectX::XMLoadFloat3(&DirectX::XMFLOAT3(D.x - B.x, D.y - B.y, D.z - B.z));
				normal = DirectX::XMVector3Cross(vec1, vec2);
			setVertices(fX, vectPoints[i + 1][j], fZ + fSpacing, i, j, vertices, indices, &index, normal, (fX + fMax) / fMax * 2, (fZ + fSpacing + fMax) / fMax * 2);
			setVertices(fX + fSpacing, vectPoints[i + 1][j + 1], fZ + fSpacing, i, j, vertices, indices, &index, normal, (fX + fSpacing + fMax) / fMax * 2, (fZ + fSpacing + fMax) / fMax * 2);
			setVertices(fX + fSpacing, vectPoints[i][j + 1], fZ, i, j, vertices, indices, &index, normal, (fX + fSpacing + fMax) / fMax * 2, (fZ + fMax) / fMax * 2);
*/
				vec1 = DirectX::XMLoadFloat3(&DirectX::XMFLOAT3(D.x - C.x, D.y - C.y, D.z - C.z));
				vec2 = DirectX::XMLoadFloat3(&DirectX::XMFLOAT3(A.x - C.x, A.y - C.y, A.z - C.z));
				normal = DirectX::XMVector3Cross(vec1, vec2);
			setVertices(fX + fSpacing, vectPoints[i + 1][j + 1], fZ + fSpacing, i, j, vertices, indices, &index, normal, (fX + fSpacing + fMax) / fMax * 2, (fZ +fSpacing + fMax) / fMax * 2);
			setVertices(fX + fSpacing, vectPoints[i][j + 1], fZ, i, j, vertices, indices, &index, normal, (fX + fSpacing + fMax) / fMax * 2, (fZ + fMax) / fMax * 2);
			setVertices(fX, vectPoints[i][j], fZ, i, j, vertices, indices, &index, normal, (fX + fMax) / fMax * 2, (fZ + fMax) / fMax * 2);

				/*vec1 = DirectX::XMLoadFloat3(&DirectX::XMFLOAT3(A.x - D.x, A.y - D.y, A.z - D.z));
				vec2 = DirectX::XMLoadFloat3(&DirectX::XMFLOAT3(B.x - D.x, B.y - D.y, B.z - D.z));
				normal = DirectX::XMVector3Cross(vec1, vec2);
			setVertices(fX + fSpacing, vectPoints[i][j + 1], fZ, i, j, vertices, indices, &index, normal, (fX + fSpacing + fMax) / fMax * 2, (fZ + fMax) / fMax * 2);
			setVertices(fX, vectPoints[i][j], fZ, i, j, vertices, indices, &index, normal, (fX + fMax) / fMax * 2, (fZ + fMax) / fMax * 2);
			setVertices(fX, vectPoints[i + 1][j], fZ + fSpacing, i, j, vertices, indices, &index, normal, (fX  + fMax) / fMax * 2, (fZ + fSpacing + fMax) / fMax * 2);
			*/


			
			
			/*setVertices(fX, vectPoints[i + 1][j], fZ + fSpacing, i, j, vertices, indices, &index);
			setVertices(fX, vectPoints[i + 1][j], fZ + fSpacing, i, j, vertices, indices, &index);
			setVertices(fX, vectPoints[i][j + 1], fZ, i, j, vertices, indices, &index);
			setVertices(fX, vectPoints[i + 1][j], fZ + fSpacing, i, j, vertices, indices, &index);*/
			
		
			


			}
//.........这里部分代码省略.........

int ParpackSolver::Solve(int nev) {
    /* Get MPI info */
    int nprocs, me;
    MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
    MPI_Comm_rank(MPI_COMM_WORLD, &me);
    MPI_Fint fcomm = MPI_Comm_c2f(MPI_COMM_WORLD);

    /* Select number of working Ritz vectors */
    if(ncv == -1)
        ncv = 2*nev;
    ncv = std::min(ncv, n-1);

    /* Initialize matrix descriptors */
    xdesc = pcontext->new_descriptor(n, 1, divup(n,nprocs), 1);
    Bdesc = pcontext->new_descriptor(n, ncv, divup(n,nprocs), ncv);
    assert(nloc == Bdesc->num_local_rows() && nloc == xdesc->num_local_rows());
    assert(ncv == Bdesc->num_local_cols() && 1 == xdesc->num_local_cols());

    /* Allocate local memory for eigenvector matrix $B$ */
    Bvalues = (real*) opsec_malloc(Bdesc->local_size() * sizeof(real));

    real sigma;
    int iparam[11], ipntr[11];

    /* Set PARPACK parameters */
    char bmat[] = "I";
    char which[] = "LA";
    char howmny[] = "All";
    iparam[0] = 1;      // ishfts
    iparam[2] = maxitr; // maxitr
    iparam[6] = 1;      // mode

    /* Allocate working memory */
    int lworkl = ncv*(ncv + 8);
    real* workl = (real*) opsec_calloc(lworkl, sizeof(real));
    real* workd = (real*) opsec_calloc(3*nloc, sizeof(real));
    real* resid = (real*) opsec_calloc(nloc, sizeof(real));
    int* select = (int*) opsec_calloc(ncv, sizeof(int));

    /* Begin reverse communication loop */
    int itr = 0;
    int info = 0;
    int ido = 0;
    while(ido != 99) {
        parpack_psaupd(&fcomm, &ido, bmat, &nloc, which, &nev,
                       &tol, resid, &ncv, Bvalues, &nloc, iparam, ipntr,
                       workd, workl, &lworkl, &info);

        if(ido == 1 || ido == -1) {
            /* Compute y = A*x (don't forget Fortran indexing conventions!) */
            slp::Matrix<real> A(Adesc, Avalues);
            slp::Matrix<real> x(xdesc, &workd[ipntr[0] - 1]);
            slp::Matrix<real> y(xdesc, &workd[ipntr[1] - 1]);
            slp::multiply(A, x, y);
        }
    }

    if(me == 0) {
        opsec_info("Number of Implicit Arnoldi update iterations taken is %d\n", iparam[2]);
        opsec_info("  info = %d\n", info);
        opsec_info("  nconv = %d, nev = %d\n", iparam[4], nev);

        time_t t = time(NULL);
        opsec_info("Time: %s\n", ctime(&t));
        opsec_info("Post-processing Ritz values and vectors\n");
    }

    /* Check return code */
    if(info < 0) {
        /* Error encountered.  Abort. */
        if(me == 0)
            opsec_error("parpack_psaupd returned error: info = %d\n", info);
        return info;
    }
    else {
        /* Save number of successfully computed eigenvalues */
        nconv = iparam[4];
        evals.resize(nconv);

        /* Retrieve eigenvalues and eigenvectors */
        int rvec = 1;
        int ierr;
        parpack_pseupd(&fcomm, &rvec, howmny, select, &evals[0], Bvalues, &nloc, &sigma,
                       bmat, &nloc, which, &nev, &tol, resid, &ncv, Bvalues, &nloc,
                       iparam, ipntr, workd, workl, &lworkl, &ierr);

        if(ierr != 0) {
            if(me == 0)
                opsec_error("parpack_pseupd returned error: ierr = %d\n", ierr);
        }
    }

    if(me == 0) {
        time_t t = time(NULL);
        opsec_info("Time: %s\n", ctime(&t));
    }

#if 0
    {
        int i;
//.........这里部分代码省略.........

    void TotalLagrangian::CalculateAll( MatrixType& rLeftHandSideMatrix,
                                        VectorType& rRightHandSideVector,
                                        ProcessInfo& rCurrentProcessInfo,
                                        bool CalculateStiffnessMatrixFlag,
                                        bool CalculateResidualVectorFlag )
    {
        KRATOS_TRY
        const unsigned int number_of_nodes = GetGeometry().size();
        const unsigned int dim = GetGeometry().WorkingSpaceDimension();
        unsigned int StrainSize;

        if ( dim == 2 )
            StrainSize = 3;
        else
            StrainSize = 6;

        Matrix B( StrainSize, number_of_nodes * dim );

        Matrix F( dim, dim );

        Matrix D( StrainSize, StrainSize );

        Matrix C( dim, dim );

        Vector StrainVector( StrainSize );

        Vector StressVector( StrainSize );

        Matrix DN_DX( number_of_nodes, dim );



        //resizing as needed the LHS
        unsigned int MatSize = number_of_nodes * dim;

        if ( CalculateStiffnessMatrixFlag == true ) //calculation of the matrix is required
        {
            if ( rLeftHandSideMatrix.size1() != MatSize )
                rLeftHandSideMatrix.resize( MatSize, MatSize, false );

            noalias( rLeftHandSideMatrix ) = ZeroMatrix( MatSize, MatSize ); //resetting LHS
        }


        //resizing as needed the RHS
        if ( CalculateResidualVectorFlag == true ) //calculation of the matrix is required
        {
            if ( rRightHandSideVector.size() != MatSize )
                rRightHandSideVector.resize( MatSize, false );

            rRightHandSideVector = ZeroVector( MatSize ); //resetting RHS
        }

        //reading integration points and local gradients
        const GeometryType::IntegrationPointsArrayType& integration_points = GetGeometry().IntegrationPoints( mThisIntegrationMethod );

        const GeometryType::ShapeFunctionsGradientsType& DN_De = GetGeometry().ShapeFunctionsLocalGradients( mThisIntegrationMethod );


        const Matrix& Ncontainer = GetGeometry().ShapeFunctionsValues( mThisIntegrationMethod );


        //calculating actual jacobian
        GeometryType::JacobiansType J;

        GetGeometry().Jacobian( J );


        //KRATOS_WATCH(J)

        //auxiliary terms
        Vector BodyForce;

        for ( unsigned int PointNumber = 0; PointNumber < integration_points.size(); PointNumber++ )
        {
            //Calculating the cartesian derivatives (it is avoided storing them to minimize storage)
            noalias( DN_DX ) = prod( DN_De[PointNumber], mInvJ0[PointNumber] );


            //deformation gradient
            noalias( F ) = prod( J[PointNumber], mInvJ0[PointNumber] );


            //strain calculation
            noalias( C ) = prod( trans( F ), F );


            CalculateStrain( C, StrainVector );


            Comprobate_State_Vector( StrainVector );
            mConstitutiveLawVector[PointNumber]->CalculateMaterialResponse(
                StrainVector,
                F,
                StressVector,
                D,
                rCurrentProcessInfo,
                GetProperties(),
                GetGeometry(),
                row( Ncontainer, PointNumber ),
//.........这里部分代码省略.........

static void panel_eid_fixup(uint16_t *mfr_name, uint16_t *product_code)
{
	B("%s: enter.\n", __func__);
	*mfr_name = 0x0101;
	*product_code = 0x0;
}

struct A;

struct B
{
    operator A();
};

struct A
{
    A(A const&){}
};

A a{B()};

/**
    Purpose
    -------
    CHEGST reduces a complex Hermitian-definite generalized
    eigenproblem to standard form.
    
    If ITYPE = 1, the problem is A*x = lambda*B*x,
    and A is overwritten by inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H)
    
    If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or
    B*A*x = lambda*x, and A is overwritten by U*A*U**H or L**H*A*L.
    
    B must have been previously factorized as U**H*U or L*L**H by CPOTRF.
    
    Arguments
    ---------
    @param[in]
    itype   INTEGER
            = 1: compute inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H);
            = 2 or 3: compute U*A*U**H or L**H*A*L.
    
    @param[in]
    uplo    magma_uplo_t
      -     = MagmaUpper:  Upper triangle of A is stored and B is factored as U**H*U;
      -     = MagmaLower:  Lower triangle of A is stored and B is factored as L*L**H.
    
    @param[in]
    n       INTEGER
            The order of the matrices A and B.  N >= 0.
    
    @param[in,out]
    A       COMPLEX array, dimension (LDA,N)
            On entry, the Hermitian matrix A.  If UPLO = MagmaUpper, the leading
            N-by-N upper triangular part of A contains the upper
            triangular part of the matrix A, and the strictly lower
            triangular part of A is not referenced.  If UPLO = MagmaLower, the
            leading N-by-N lower triangular part of A contains the lower
            triangular part of the matrix A, and the strictly upper
            triangular part of A is not referenced.
    \n
            On exit, if INFO = 0, the transformed matrix, stored in the
            same format as A.
    
    @param[in]
    lda     INTEGER
            The leading dimension of the array A.  LDA >= max(1,N).
    
    @param[in]
    B       COMPLEX array, dimension (LDB,N)
            The triangular factor from the Cholesky factorization of B,
            as returned by CPOTRF.
    
    @param[in]
    ldb     INTEGER
            The leading dimension of the array B.  LDB >= max(1,N).
    
    @param[out]
    info    INTEGER
      -     = 0:  successful exit
      -     < 0:  if INFO = -i, the i-th argument had an illegal value

    @ingroup magma_cheev_comp
    ********************************************************************/
extern "C" magma_int_t
magma_chegst(
    magma_int_t itype, magma_uplo_t uplo, magma_int_t n,
    magmaFloatComplex *A, magma_int_t lda,
    magmaFloatComplex *B, magma_int_t ldb,
    magma_int_t *info)
{
#define A(i, j) (A + (j)*lda + (i))
#define B(i, j) (B + (j)*ldb + (i))

#define dA(i, j) (dw + (j)*ldda + (i))
#define dB(i, j) (dw + n*ldda + (j)*lddb + (i))

    const char* uplo_ = lapack_uplo_const( uplo );
    magma_int_t        nb;
    magma_int_t        k, kb, kb2;
    magmaFloatComplex    c_one      = MAGMA_C_ONE;
    magmaFloatComplex    c_neg_one  = MAGMA_C_NEG_ONE;
    magmaFloatComplex    c_half     = MAGMA_C_HALF;
    magmaFloatComplex    c_neg_half = MAGMA_C_NEG_HALF;
    magmaFloatComplex   *dw;
    magma_int_t        ldda = n;
    magma_int_t        lddb = n;
    float             d_one = 1.0;
    int upper = (uplo == MagmaUpper);
    
    /* Test the input parameters. */
    *info = 0;
    if (itype < 1 || itype > 3) {
        *info = -1;
    } else if (! upper && uplo != MagmaLower) {
        *info = -2;
    } else if (n < 0) {
        *info = -3;
    } else if (lda < max(1,n)) {
        *info = -5;
    } else if (ldb < max(1,n)) {
//.........这里部分代码省略.........

#include "../konoha_mpi.h"

/* ------------------------------------------------------------------------ */
/* TYPEMAP */

TYPEMAP Int_MPIData(CTX ctx, ksfp_t *sfp _RIX)
{
	MPID(data, new_O(MPIData, knh_getcid(ctx, B("konoha.mpi.MPIData"))));
	MPID_INIT(data, new_Int(ctx, sfp[1].ivalue), MPI_LONG, CLASS_Int, O_cid(sfp[1].o));
	RETURN_(data);
}

TYPEMAP int___MPIData(CTX ctx, ksfp_t *sfp _RIX)
{
	MPID(data, new_O(MPIData, knh_getcid(ctx, B("konoha.mpi.MPIData"))));
	MPID_INIT(data, sfp[1].a, MPI_LONG, CLASS_Array, O_cid(sfp[1].o));
	RETURN_(data);
}

TYPEMAP Float_MPIData(CTX ctx, ksfp_t *sfp _RIX)
{
	MPID(data, new_O(MPIData, knh_getcid(ctx, B("konoha.mpi.MPIData"))));
	MPID_INIT(data, new_Float_(ctx, CLASS_Float, sfp[1].fvalue), MPI_DOUBLE, CLASS_Float, O_cid(sfp[1].o));
	RETURN_(data);
}

TYPEMAP float___MPIData(CTX ctx, ksfp_t *sfp _RIX)
{
	MPID(data, new_O(MPIData, knh_getcid(ctx, B("konoha.mpi.MPIData"))));
	MPID_INIT(data, sfp[1].a, MPI_DOUBLE, CLASS_Array, O_cid(sfp[1].o));
	RETURN_(data);

int A(int x)
{
	int arr[5];
    return B(arr[4]);
}

int main(int argc, char *argv[])
{
  int my_ID, myrow, mycol, /* my index and row and column index    */
      root=0,           /* ID of root rank                         */
      Num_procs,        /* number of ranks                         */
      nprow, npcol,     /* row, column dimensions of rank grid     */
      order,            /* matrix order                            */
      mynrows, myfrow,  /* my number of rows and index of first row*/
      mylrow,           /* and last row                            */
    /*myncols,*/ myfcol,/* my number of cols and index of first row*/
      mylcol,           /* and last row                            */
      *mm,              /* arrays that hold m_i's and n_j's        */
      *nn,
    /*nb,*/             /* block factor for SUMMA                  */
      inner_block_flag, /* flag to select local DGEMM blocking     */
      error=0,          /* error flag                              */
      *ranks,           /* work array for row and column ranks     */
    /*lda, ldb, ldc,*/  /* leading array dimensions of a, b, and c */
      iter, iterations;
  long lda, ldb, ldc,
       nb, myncols;     /* make long to avoid integer overflow     */
  double RESTRICT *a, *b, *c,    /* arrays that hold local a, b, c */
      *work1, *work2,   /* work arrays to pass to dpmmmult         */
      local_dgemm_time, /* timing parameters                       */
      dgemm_time,
      avgtime; 
  double
      forder, nflops,   /* float matrix order + total flops        */
      checksum,         /* array checksum for verification test    */
      checksum_local=0.0,
      ref_checksum;     /* reference checkcum for verification     */
  MPI_Group world_group, 
      temp_group;
  MPI_Comm comm_row,    /* communicators for row and column ranks  */
      comm_col;         /* of rank grid                            */
  int shortcut;         /* true if only doing initialization       */

  /* initialize                                                    */
  MPI_Init(&argc,&argv);
  MPI_Comm_rank( MPI_COMM_WORLD, &my_ID );
  MPI_Comm_size( MPI_COMM_WORLD, &Num_procs );

/*********************************************************************
** process, test and broadcast input parameters
*********************************************************************/

  if (my_ID == root) {
    printf("Parallel Research Kernels version %s\n", PRKVERSION);
    printf("MPI Dense matrix-matrix multiplication: C = A x B\n");

    if (argc != 5) {
      printf("Usage: %s <# iterations> <matrix order> <outer block size> ",
                                                               *argv);
      printf("<local block flag (non-zero=yes, zero=no)>\n");
      error = 1;
      goto ENDOFTESTS;
    }

    iterations  = atoi(*++argv);
    if(iterations < 1){
      printf("ERROR: iterations must be positive: %d \n",iterations);
      error = 1;
      goto ENDOFTESTS;
    }

    order = atoi(*++argv);
    if (order < 0) {
      shortcut = 1;
      order    = -order;
    } else shortcut = 0;
    if (order < Num_procs) {
      printf("ERROR: matrix order too small: %d\n", order);
      error = 1;
      goto ENDOFTESTS;
    }

    nb = atol(*++argv);
    /* a non-positive tile size means no outer level tiling        */

    inner_block_flag = atoi(*++argv);
    
    ENDOFTESTS:;
  }
  bail_out(error);

  MPI_Bcast(&order,            1, MPI_INT,  root, MPI_COMM_WORLD);
  MPI_Bcast(&iterations,       1, MPI_INT,  root, MPI_COMM_WORLD);
  MPI_Bcast(&nb,               1, MPI_LONG, root, MPI_COMM_WORLD);
  MPI_Bcast(&shortcut,         1, MPI_INT,  root, MPI_COMM_WORLD);
  MPI_Bcast(&inner_block_flag, 1, MPI_INT,  root, MPI_COMM_WORLD);

  /* compute rank grid to most closely match a square; to do so,
     compute largest divisor of Num_procs, using hare-brained method. 
     The small term epsilon is used to guard against roundoff errors 
     in case Num_procs is a perfect square                         */
  nprow = (int) (sqrt((double) Num_procs + epsilon));
  while (Num_procs%nprow) nprow--;
  npcol = Num_procs/nprow;

  if (my_ID == root) {
//.........这里部分代码省略.........

unsigned int fast_inclusion(const tetra *t, const point *p) {

    auto &a = *t->p[0],
         &b = *t->p[1],
         &c = *t->p[2],
         &d = *t->p[3];

    Eigen::Matrix<double, 3, 3> A;

    A.col(0) << b.x - a.x, b.y - a.y, b.z - a.z;
    A.col(1) << c.x - a.x, c.y - a.y, c.z - a.z;
    A.col(2) << d.x - a.x, d.y - a.y, d.z - a.z;

    if (verbose == 3)
        std::cout << A << std::endl;

    Eigen::Matrix<double, 3, 1> x, B(p->x - a.x, p->y - a.y, p->z - a.z);

    x = A.inverse() * B;

    double sum = 0;

    for (unsigned int i = 0; i < 3; ++i) {

        if (std::abs(x[i]) < 1e-10)
            x[i] = 0;

        if (x[i] < 0)
            return 0; // outside
        else
            sum += x[i];
    }


    if (std::abs(sum - 1) < 1e-10)
        sum = 1;//return exact_inclusion(t, p);

    if (std::abs(sum) < 1e-10)
        sum = 0;

    if (sum > 1)
        return 0; // outside

    if (sum == 0)
        return 1; // vertex 0

    double u(x[0]), v(x[1]), w(x[2]);

    if (u == 1) {

        return 2; // vertex 1
    }

    else if (u > 0) {

        if (v > 0) {
            
            if (w > 0) {

                if (sum == 1)
                    return 14; // surface 321
                else
                    return 15; // inside
            }

            else {

                if (sum == 1)
                    return 6; // edge 21
                else
                    return 7; // surface 012
            }
        }

        else {
            
            if (w > 0) {

                if (sum == 1)
                    return 10; // edge 31
                else
                    return 11; // surface 031
            }

            else {

                return 3; // edge 10
            }
        } 
    } else {

        if (v == 1)
            return 4; // vertex 2

        else if (v > 0) {

            if (w > 0) {

                if (sum == 1)
                    return 12; // edge 32
//.........这里部分代码省略.........

extern "C" magma_int_t
magma_zhegst_gpu(magma_int_t itype, char uplo, magma_int_t n,
                 cuDoubleComplex *da, magma_int_t ldda,
                 cuDoubleComplex *db, magma_int_t lddb, magma_int_t *info)
{
/*
  -- MAGMA (version 1.3.0) --
     Univ. of Tennessee, Knoxville
     Univ. of California, Berkeley
     Univ. of Colorado, Denver
     November 2012
 
   Purpose
   =======
   ZHEGST_GPU reduces a complex Hermitian-definite generalized
   eigenproblem to standard form.
   
   If ITYPE = 1, the problem is A*x = lambda*B*x,
   and A is overwritten by inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H)
   
   If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or
   B*A*x = lambda*x, and A is overwritten by U*A*U**H or L**H*A*L.
   
   B must have been previously factorized as U**H*U or L*L**H by ZPOTRF.
   
   Arguments
   =========
   ITYPE   (input) INTEGER
           = 1: compute inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H);
           = 2 or 3: compute U*A*U**H or L**H*A*L.
   
   UPLO    (input) CHARACTER*1
           = 'U':  Upper triangle of A is stored and B is factored as
                   U**H*U;
           = 'L':  Lower triangle of A is stored and B is factored as
                   L*L**H.
   
   N       (input) INTEGER
           The order of the matrices A and B.  N >= 0.
   
   DA      (device input/output) COMPLEX*16 array, dimension (LDA,N)
           On entry, the Hermitian matrix A.  If UPLO = 'U', the leading
           N-by-N upper triangular part of A contains the upper
           triangular part of the matrix A, and the strictly lower
           triangular part of A is not referenced.  If UPLO = 'L', the
           leading N-by-N lower triangular part of A contains the lower
           triangular part of the matrix A, and the strictly upper
           triangular part of A is not referenced.
   
           On exit, if INFO = 0, the transformed matrix, stored in the
           same format as A.
   
   LDDA    (input) INTEGER
           The leading dimension of the array A.  LDA >= max(1,N).
   
   DB      (device input) COMPLEX*16 array, dimension (LDB,N)
           The triangular factor from the Cholesky factorization of B,
           as returned by ZPOTRF.
   
   LDDB    (input) INTEGER
           The leading dimension of the array B.  LDB >= max(1,N).
   
   INFO    (output) INTEGER
           = 0:  successful exit
           < 0:  if INFO = -i, the i-th argument had an illegal value
   =====================================================================*/
  
  char uplo_[2] = {uplo, 0};
  magma_int_t        nb;
  magma_int_t        k, kb, kb2;
  cuDoubleComplex    c_one      = MAGMA_Z_ONE;
  cuDoubleComplex    c_neg_one  = MAGMA_Z_NEG_ONE;
  cuDoubleComplex    c_half     = MAGMA_Z_HALF;
  cuDoubleComplex    c_neg_half = MAGMA_Z_NEG_HALF;
  cuDoubleComplex   *w;
  magma_int_t        lda;
  magma_int_t        ldb;
  double             d_one = 1.0;
  int upper = lapackf77_lsame(uplo_, "U");
  
  /* Test the input parameters. */
  *info = 0;
  if (itype<1 || itype>3){
    *info = -1;
  }else if ((! upper) && (! lapackf77_lsame(uplo_, "L"))) {
    *info = -2;
  } else if (n < 0) {
    *info = -3;
  } else if (ldda < max(1,n)) {
    *info = -5;
  }else if (lddb < max(1,n)) {
    *info = -7;
  }
    if (*info != 0) {
        magma_xerbla( __func__, -(*info) );
        return *info;
    }
  
  /* Quick return */
  if ( n == 0 )
//.........这里部分代码省略.........

 void foo (int) {
   struct B { virtual ~B() {} };
   B();
 }

LR::LRSplineSurface* ASMu2D::scRecovery (const IntegrandBase& integrand) const
{
  PROFILE2("ASMu2D::scRecovery");

  const int m = integrand.derivativeOrder();
  const int p1 = lrspline->order(0);
  const int p2 = lrspline->order(1);

  // Get Gaussian quadrature point coordinates
  const int ng1 = p1 - m;
  const int ng2 = p2 - m;
  const double* xg = GaussQuadrature::getCoord(ng1);
  const double* yg = GaussQuadrature::getCoord(ng2);
  if (!xg || !yg) return nullptr;

  // Compute parameter values of the Greville points
  std::array<RealArray,2> gpar;
  if (!this->getGrevilleParameters(gpar[0],0)) return nullptr;
  if (!this->getGrevilleParameters(gpar[1],1)) return nullptr;

  const int n1 = p1 - m + 1; // Patch size in first parameter direction
  const int n2 = p2 - m + 1; // Patch size in second parameter direction

  const size_t nCmp = integrand.getNoFields(); // Number of result components
  const size_t nPol = n1*n2; // Number of terms in polynomial expansion

  Matrix sValues(nCmp,gpar[0].size());
  Vector P(nPol);
  Go::Point X, G;

  // Loop over all Greville points (one for each basis function)
  size_t k, l, ip = 0;
  std::vector<LR::Element*>::const_iterator elStart, elEnd, el;
  std::vector<LR::Element*> supportElements;
  for (LR::Basisfunction *b : lrspline->getAllBasisfunctions())
  {
#if SP_DEBUG > 2
    std::cout <<"Basis: "<< *b <<"\n  ng1 ="<< ng1 <<"\n  ng2 ="<< ng2
              <<"\n  nPol="<< nPol << std::endl;
#endif

    // Special case for basis functions with too many zero knot spans by using
    // the extended support
    // if(nel*ng1*ng2 < nPol)
    if(true)
    {
      // KMO: Here I'm not sure how this will change when m > 1.
      // In that case I think we would need smaller patches (as in the tensor
      // splines case). But how to do that???
      supportElements = b->getExtendedSupport();
      elStart = supportElements.begin();
      elEnd   = supportElements.end();
#if SP_DEBUG > 2
      std::cout <<"Extended basis:";
      for (el = elStart; el != elEnd; el++)
        std::cout <<"\n  " << **el;
      std::cout << std::endl;
#endif
    }
    else
    {
      elStart = b->supportedElementBegin();
      elEnd   = b->supportedElementEnd();
    }

    // Physical coordinates of current Greville point
    lrspline->point(G,gpar[0][ip],gpar[1][ip]);

    // Set up the local projection matrices
    DenseMatrix A(nPol,nPol);
    Matrix B(nPol,nCmp);

    // Loop over all non-zero knot-spans in the support of
    // the basis function associated with current Greville point
    for (el = elStart; el != elEnd; el++)
    {
      int iel = (**el).getId()+1;

      // evaluate all gauss points for this element
      std::array<RealArray,2> gaussPt, unstrGauss;
      this->getGaussPointParameters(gaussPt[0],0,ng1,iel,xg);
      this->getGaussPointParameters(gaussPt[1],1,ng2,iel,yg);

#if SP_DEBUG > 2
      std::cout << "Element " << **el << std::endl;
#endif

      // convert to unstructred mesh representation
      expandTensorGrid(gaussPt.data(),unstrGauss.data());

      // Evaluate the secondary solution at all Gauss points
      Matrix sField;
      if (!this->evalSolution(sField,integrand,unstrGauss.data()))
        return nullptr;

      // Loop over the Gauss points in current knot-span
      int i, j, ig = 1;
      for (j = 0; j < ng2; j++)
	for (i = 0; i < ng1; i++, ig++)
	{
//.........这里部分代码省略.........

void AddDot4x4( int k, double *a, int lda,  double *b, int ldb, double *c, int ldc )
{
  /* So, this routine computes a 4x4 block of matrix A

           C( 0, 0 ), C( 0, 1 ), C( 0, 2 ), C( 0, 3 ).  
           C( 1, 0 ), C( 1, 1 ), C( 1, 2 ), C( 1, 3 ).  
           C( 2, 0 ), C( 2, 1 ), C( 2, 2 ), C( 2, 3 ).  
           C( 3, 0 ), C( 3, 1 ), C( 3, 2 ), C( 3, 3 ).  

     Notice that this routine is called with c = C( i, j ) in the
     previous routine, so these are actually the elements 

           C( i  , j ), C( i  , j+1 ), C( i  , j+2 ), C( i  , j+3 ) 
           C( i+1, j ), C( i+1, j+1 ), C( i+1, j+2 ), C( i+1, j+3 ) 
           C( i+2, j ), C( i+2, j+1 ), C( i+2, j+2 ), C( i+2, j+3 ) 
           C( i+3, j ), C( i+3, j+1 ), C( i+3, j+2 ), C( i+3, j+3 ) 
	  
     in the original matrix C 

     In this version, we merge each set of four loops, computing four
     inner products simultaneously. */

  int p;

  for ( p=0; p<k; p++ ){
    /* First row */
    C( 0, 0 ) += A( 0, p ) * B( p, 0 );     
    C( 0, 1 ) += A( 0, p ) * B( p, 1 );     
    C( 0, 2 ) += A( 0, p ) * B( p, 2 );     
    C( 0, 3 ) += A( 0, p ) * B( p, 3 );     

    /* Second row */
    C( 1, 0 ) += A( 1, p ) * B( p, 0 );     
    C( 1, 1 ) += A( 1, p ) * B( p, 1 );     
    C( 1, 2 ) += A( 1, p ) * B( p, 2 );     
    C( 1, 3 ) += A( 1, p ) * B( p, 3 );     

    /* Third row */
    C( 2, 0 ) += A( 2, p ) * B( p, 0 );     
    C( 2, 1 ) += A( 2, p ) * B( p, 1 );     
    C( 2, 2 ) += A( 2, p ) * B( p, 2 );     
    C( 2, 3 ) += A( 2, p ) * B( p, 3 );     

    /* Fourth row */
    C( 3, 0 ) += A( 3, p ) * B( p, 0 );     
    C( 3, 1 ) += A( 3, p ) * B( p, 1 );     
    C( 3, 2 ) += A( 3, p ) * B( p, 2 );     
    C( 3, 3 ) += A( 3, p ) * B( p, 3 );     
  }
}