double quadform(double *x, double *A, int N, int incx, int LDA)
{
  
  //quadform2 seems to be faster
  return( quadform2(x, A, N, incx, LDA) );
  
  int Nsqr = N*N,info,i=0,j=0;
  double *B = Calloc(Nsqr,double);
  //double one=1;
  //double zero=0;
  double sumSq=0;
  double y[N];
  int iOne=1;

  for(i=0;i<N;i++){
    y[i] = x[i*incx];
  }
  for(i=0;i<N;i++){
	Memcpy(&B[i*N],&A[i*LDA],N);
  }

  F77_NAME(dpotrf)("U", &N, B, &N, &info);
  F77_NAME(dtrmv)("U","N","N", &N, B, &N, y, &iOne);
  
  for(i=0;i<N;i++){
    sumSq += y[i]*y[i];
  }
  
  Free(B);
  
  return(sumSq);
}

VertexBufferPtr D3D9VideoBufferManager::CreateVertexBuffer(
	int size, int stride, USAGE usage, bool cpuAccess, const void * initData)
{
    d_assert (size > 0);

    HRESULT hr;
    DWORD D3DUsage = D3D9Mapping::GetD3DUsage(usage);
    D3DPOOL D3DPool = D3D9Mapping::GetD3DPool(usage);
    IDirect3DVertexBuffer9 * pD3DVB;

    hr = mD3D9Device->CreateVertexBuffer(size, D3DUsage, 0, D3DPool, &pD3DVB, NULL);

    if (FAILED(hr))
    {
        EXCEPTION("D3D Error: CreateVertexBuffer failed, desc: " + D3D9Mapping::GetD3DErrorDescription(hr));
    }

    D3D9VertexBuffer * pVB = new D3D9VertexBuffer(mD3D9Device);

    pVB->mD3D9VertexBuffer = pD3DVB;
    pVB->mSize = size;
	pVB->mStride = stride;
    pVB->mUsage = usage;

    mVertexBuffers.PushBack(pVB);

	if (initData != NULL)
	{
		void * data = pVB->Lock(0, 0, LOCK_DISCARD);
		Memcpy(data, initData, size);
		pVB->Unlock();
	}
	
    return VertexBufferPtr(pVB);
}

static
double* check_gv(SEXP gr, SEXP hs, SEXP rho, int n, double *gv, double *hv)
{
    SEXP gval = PROTECT(coerceVector(eval(gr, rho), REALSXP));
    if (LENGTH(gval) != n)
	error(_("gradient function must return a numeric vector of length %d"), n);
    Memcpy(gv, REAL(gval), n);
    for (int i = 0; i < n; i++)
	if(ISNAN(gv[i])) error("NA/NaN gradient evaluation");
    if (hv) {
	SEXP hval = PROTECT(eval(hs, rho));
	SEXP dim = getAttrib(hval, R_DimSymbol);
	int i, j, pos;
	double *rhval = REAL(hval);

	if (!isReal(hval) || LENGTH(dim) != 2 ||
	    INTEGER(dim)[0] != n || INTEGER(dim)[1] != n)
	    error(_("Hessian function must return a square numeric matrix of order %d"),
		  n);
	for (i = 0, pos = 0; i < n; i++) /* copy lower triangle row-wise */
	    for (j = 0; j <= i; j++) {
		hv[pos] = rhval[i + j * n];
		if(ISNAN(hv[pos])) error("NA/NaN Hessian evaluation");
		pos++;
	    }
	UNPROTECT(1);
    }
    UNPROTECT(1);
    return gv;
}

void CopyVertex(PolygonGroup *srcGroup, uint32 srcPos, PolygonGroup *dstGroup, uint32 dstPos)
{
    int32 srcFormat = srcGroup->GetFormat();
    int32 dstFormat = dstGroup->GetFormat();
    int32 copyFormat = srcFormat&dstFormat;    //most common format;   

    uint8 *srcData = srcGroup->meshData+srcPos*GetVertexSize(srcFormat);
    uint8 *dstData = dstGroup->meshData+dstPos*GetVertexSize(dstFormat);
        
    for (uint32 mask = EVF_LOWER_BIT; mask <= EVF_HIGHER_BIT; mask = mask << 1)
    {
        int32 vertexAttribSize = GetVertexSize(mask);
        if (mask&copyFormat)
            Memcpy(dstData, srcData, vertexAttribSize);

        if (mask&srcFormat)
            srcData+=vertexAttribSize;
        
        if (mask&dstFormat)
            dstData+=vertexAttribSize;

         copyFormat&=~mask;
    }    
    
    /*unsupported stream*/
    DVASSERT((copyFormat == 0)&&"Unsupported attribute stream in copy");
    
}

UINT CMemFile::Read(void* lpBuf, UINT nCount)
{
	ASSERT_VALID(this);

	if (nCount == 0)
		return 0;

	ASSERT(lpBuf != NULL);
	ASSERT(AfxIsValidAddress(lpBuf, nCount));

	if (m_nPosition > m_nFileSize)
		return 0;

	UINT nRead;
	if (m_nPosition + nCount > m_nFileSize)
		nRead = (UINT)(m_nFileSize - m_nPosition);
	else
		nRead = nCount;

	Memcpy((BYTE*)lpBuf, (BYTE*)m_lpBuffer + m_nPosition, nRead);
	m_nPosition += nRead;

	ASSERT_VALID(this);

	return nRead;
}

 void CD3DDevice::UpdateBuffer( IGfxBuffer* pBuffer, void* pData, uint nSize )
 {
     D3D11_MAPPED_SUBRESOURCE pMappedBuffer;
     m_pContext->Map( ((CD3DBuffer*)pBuffer)->m_pBuffer, 0, D3D11_MAP_WRITE_DISCARD, 0, &pMappedBuffer );
     Memcpy( pMappedBuffer.pData, pData, nSize );
     m_pContext->Unmap( ((CD3DBuffer*)pBuffer)->m_pBuffer, 0 );
 }

int32
cache_writesector(
    IN  tcache_t * pcache,
    IN  int32 addr,
    IN  ubyte * ptr)
{
	int16 secindex;
	int32 ret;

	ret = CACHE_OK;

	if (!pcache->use_cache) {
		return HAI_readsector(pcache->hdev, addr, ptr);
	}

	if (!_find_in_list(pcache, addr, &secindex)) {
		if ((secindex = _do_cache_miss(pcache, secindex, addr)) < 0) {
			ret = secindex;
		}
	}
	else {
		_do_cache_hint(pcache, secindex);
	}

	if (ret >= 0) {
		Memcpy(pcache->secbufs[secindex].secbuf, ptr, pcache->sector_size);
		pcache->secbufs[secindex].is_dirty = 1;
	}

	return ret;
}

// Compute determinant of an N by N positive definite matrix A
double matrixDet(double *A, int N, int LDA, int doLog, int *info)
{
//SUBROUTINE DPOTRF( UPLO, N, A, LDA, INFO )
	int i=0;
	double B[N*N], logDet=0;
	
	//Memcpy(B,A,N*N);
	for(i=0;i<N;i++){
		Memcpy(&B[i*N],&A[i*LDA],N);
	}
	
	F77_CALL(dpotrf)("U", &N, B, &N, info);
	if(*info){
		Rprintf("Cholesky decomposition in matrixDet() returned nonzero info %d.\n",info);
	}
	
	for(i=0;i<N;i++)
	{
		logDet += 2 * log(B[i*N+i]);
	}
	
	if(doLog){
		return(logDet);
	}else{
		return(exp(logDet));
	}
}

// Modified version of Tim Davis's cs_qr_mex.c file for MATLAB (in CSparse)
//  Usage: [V,beta,p,R,q] = cs_qr(A) ;
SEXP dgCMatrix_QR(SEXP Ap, SEXP order)
{
    CSP A = AS_CSP__(Ap), D;
    int io = INTEGER(order)[0];
    Rboolean verbose = (io < 0);
    int m = A->m, n = A->n, ord = asLogical(order) ? 3 : 0, *p;
    R_CheckStack();

    if (m < n) error(_("A must have #{rows} >= #{columns}")) ;
    SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("sparseQR")));
    int *dims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));
    dims[0] = m; dims[1] = n;
    css *S = cs_sqr(ord, A, 1);	/* symbolic QR ordering & analysis*/
    if (!S) error(_("cs_sqr failed"));
    if(verbose && S->m2 > m) // in ./cs.h , m2 := # of rows for QR, after adding fictitious rows
	Rprintf("Symbolic QR(): Matrix structurally rank deficient (m2-m = %d)\n",
		S->m2 - m);
    csn *N = cs_qr(A, S);		/* numeric QR factorization */
    if (!N) error(_("cs_qr failed")) ;
    cs_dropzeros(N->L);		/* drop zeros from V and sort */
    D = cs_transpose(N->L, 1); cs_spfree(N->L);
    N->L = cs_transpose(D, 1); cs_spfree(D);
    cs_dropzeros(N->U);		/* drop zeros from R and sort */
    D = cs_transpose(N->U, 1); cs_spfree(N->U) ;
    N->U = cs_transpose(D, 1); cs_spfree(D);
    m = N->L->m;		/* m may be larger now */
    // MM: m := S->m2  also counting the ficticious rows (Tim Davis, p.72, 74f)
    p = cs_pinv(S->pinv, m);	/* p = pinv' */
    SET_SLOT(ans, install("V"),
	     Matrix_cs_to_SEXP(N->L, "dgCMatrix", 0));
    Memcpy(REAL(ALLOC_SLOT(ans, install("beta"),
			   REALSXP, n)), N->B, n);
    Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_pSym,
			      INTSXP, m)), p, m);
    SET_SLOT(ans, install("R"),
	     Matrix_cs_to_SEXP(N->U, "dgCMatrix", 0));
    if (ord)
	Memcpy(INTEGER(ALLOC_SLOT(ans, install("q"),
				  INTSXP, n)), S->q, n);
    else
	ALLOC_SLOT(ans, install("q"), INTSXP, 0);
    cs_nfree(N);
    cs_sfree(S);
    cs_free(p);
    UNPROTECT(1);
    return ans;
}

/**
   Copy cholmod_triplet to an R_alloc()ed version of it
 */
static void chTr2Ralloc(CHM_TR dest, CHM_TR src)
{
    int nnz;

    /* copy all the (non-pointer) characteristics of src to dest */
    memcpy(dest, src, sizeof(cholmod_triplet));

    /* R_alloc the vector storage for dest and copy the contents from src */
    nnz = src->nnz;
    dest->i = (void*) Memcpy((int*)R_alloc(nnz, sizeof(int)),
			     (int*)(src->i), nnz);
    dest->j = (void*) Memcpy((int*)R_alloc(nnz, sizeof(int)),
			     (int*)(src->j), nnz);
    if(src->xtype)
	dest->x = (void*) Memcpy((double*)R_alloc(nnz, sizeof(double)),
				 (double*)(src->x), nnz);
}

SEXP dtCMatrix_sparse_solve(SEXP a, SEXP b)
{
    SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgCMatrix")));
    CSP A = AS_CSP(a), B = AS_CSP(b);
    int *xp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, (B->n) + 1)),
	xnz = 10 * B->p[B->n];	/* initial estimate of nnz in x */
    int *ti = Calloc(xnz, int), k, lo = uplo_P(a)[0] == 'L', pos = 0;
    double *tx = Calloc(xnz, double);
    double  *wrk = Alloca(A->n, double);
    int *xi = Alloca(2*A->n, int);	/* for cs_reach */
    R_CheckStack();

    if (A->m != A->n || B->n < 1 || A->n < 1 || A->n != B->m)
	error(_("Dimensions of system to be solved are inconsistent"));
    slot_dup(ans, b, Matrix_DimSym);
    SET_DimNames(ans, b);
    xp[0] = 0;
    for (k = 0; k < B->n; k++) {
	int top = cs_spsolve (A, B, k, xi, wrk, (int *)NULL, lo);
	int nz = A->n - top, p;

	xp[k + 1] = nz + xp[k];
	if (xp[k + 1] > xnz) {
	    while (xp[k + 1] > xnz) xnz *= 2;
	    ti = Realloc(ti, xnz, int);
	    tx = Realloc(tx, xnz, double);
	}
	if (lo)			/* increasing row order */
	    for(p = top; p < A->n; p++, pos++) {
		ti[pos] = xi[p];
		tx[pos] = wrk[xi[p]];
	    }
	else			/* decreasing order, reverse copy */
	    for(p = A->n - 1; p >= top; p--, pos++) {
		ti[pos] = xi[p];
		tx[pos] = wrk[xi[p]];
	    }
    }
    xnz = xp[B->n];
    Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP,  xnz)), ti, xnz);
    Memcpy(   REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, xnz)), tx, xnz);

    Free(ti); Free(tx);
    UNPROTECT(1);
    return ans;
}

/*parse the response message and manipulate the send message*/
static int handlemessage(struct dhcp_msg*rmsg,struct dhcp_msg*smsg)
{
	int i,j,p;
	/*find DNS info*/
	if(getdns){
		p=messagefindoption(rmsg,OPT_DNS_SERVER);
		if(p>=0){
			for(i=0;i<rmsg->msg_opt[p].opt_dns_server.num_dns;i++)
				addaddr(dnsservers,rmsg->msg_opt[p].opt_dns_server.addr[i]);
		}
		p=messagefindoption(rmsg,OPT_DNS_NAME);
		if(p>=0){
			for(i=0;i<rmsg->msg_opt[p].opt_dns_name.num_dns;i++)
				adddomain(rmsg->msg_opt[p].opt_dns_name.namelist[i]);
		}
	}
	/*find PREFIX info*/
	if(getprefix){
		p=messagefindoption(rmsg,OPT_IAPD);
		if(p>=0)
		for(i=0;i<rmsg->msg_opt[p].opt_numopts;i++)
			if(rmsg->msg_opt[p].subopt[i].opt_type==OPT_IAPREFIX){
				j=addaddr(prefixes,rmsg->msg_opt[p].subopt[i].opt_iaprefix.prefix);
				if(j>=0)
					prefixlens[j]=rmsg->msg_opt[p].subopt[i].opt_iaprefix.prefixlen;
			}
	}
	/*find IANA info*/
	if(getaddress){
		p=messagefindoption(rmsg,OPT_IANA);
		if(p>=0)
		for(i=0;i<rmsg->msg_opt[p].opt_numopts;i++)
			if(rmsg->msg_opt[p].subopt[i].opt_type==OPT_IAADDR)
				addaddr(addresses,rmsg->msg_opt[p].subopt[i].opt_iaaddress.addr);
	}
	/*copy server address*/
	Memcpy(&dhcpserver,&rmsg->msg_peer.sin6_addr,16);
	/*check for rapid commit or type=REPLY; if so: tell caller it can stop now*/
	if(rmsg->msg_type==MSG_REPLY)return 0;
	if(messagefindoption(rmsg,OPT_RAPIDCOMMIT)>=0)return 0;
	/*otherwise we need to continue*/
	/*correct message type & id*/
	clearrecvfilter();
	if(getprefix||getaddress){
		addrecvfilter(MSG_REPLY);
		smsg->msg_type=MSG_REQUEST;
	}else{
		addrecvfilter(MSG_REPLY);
		smsg->msg_type=MSG_IREQUEST;
	}
	smsg->msg_id++; /*elapsed time continues to count*/
	/*rapid commit is no longer applicable*/
	messageremoveoption(smsg,OPT_RAPIDCOMMIT);
	/*append server ID*/
	p=messagefindoption(rmsg,OPT_SERVERID);
	if(p>=0)messageappendopt(smsg,&rmsg->msg_opt[p]);
	return 1;
}

/* TODO(gauravsh): This could easily be integrated into KeyBlockCreate()
 * since the code is almost a mirror - I have kept it as such to avoid changing
 * the existing interface. */
VbKeyBlockHeader* KeyBlockCreate_external(const VbPublicKey* data_key,
                                          const char* signing_key_pem_file,
                                          uint64_t algorithm,
                                          uint64_t flags,
                                          const char* external_signer) {
  VbKeyBlockHeader* h;
  uint64_t signed_size = sizeof(VbKeyBlockHeader) + data_key->key_size;
  uint64_t block_size = (signed_size + SHA512_DIGEST_SIZE +
                         siglen_map[algorithm]);
  uint8_t* data_key_dest;
  uint8_t* block_sig_dest;
  uint8_t* block_chk_dest;
  VbSignature *sigtmp;

  /* Allocate key block */
  h = (VbKeyBlockHeader*)malloc(block_size);
  if (!h)
    return NULL;
  if (!signing_key_pem_file || !data_key || !external_signer)
    return NULL;

  data_key_dest = (uint8_t*)(h + 1);
  block_chk_dest = data_key_dest + data_key->key_size;
  block_sig_dest = block_chk_dest + SHA512_DIGEST_SIZE;

  Memcpy(h->magic, KEY_BLOCK_MAGIC, KEY_BLOCK_MAGIC_SIZE);
  h->header_version_major = KEY_BLOCK_HEADER_VERSION_MAJOR;
  h->header_version_minor = KEY_BLOCK_HEADER_VERSION_MINOR;
  h->key_block_size = block_size;
  h->key_block_flags = flags;

  /* Copy data key */
  PublicKeyInit(&h->data_key, data_key_dest, data_key->key_size);
  PublicKeyCopy(&h->data_key, data_key);

  /* Set up signature structs so we can calculate the signatures */
  SignatureInit(&h->key_block_checksum, block_chk_dest,
                SHA512_DIGEST_SIZE, signed_size);
  SignatureInit(&h->key_block_signature, block_sig_dest,
                siglen_map[algorithm], signed_size);

  /* Calculate checksum */
  sigtmp = CalculateChecksum((uint8_t*)h, signed_size);
  SignatureCopy(&h->key_block_checksum, sigtmp);
  free(sigtmp);

  /* Calculate signature */
  sigtmp = CalculateSignature_external((uint8_t*)h, signed_size,
                                       signing_key_pem_file, algorithm,
                                       external_signer);
  SignatureCopy(&h->key_block_signature, sigtmp);
  free(sigtmp);

  /* Return the header */
  return h;
}

void CMemFile::Write( const void *lpBuf, UINT nCount )
/****************************************************/
{
    if( m_nPosition + nCount >= m_nFileSize ) {
        m_nFileSize = m_nPosition + nCount;
        GrowFile( m_nFileSize );
    }
    Memcpy( m_lpBuffer + m_nPosition, (const BYTE *)lpBuf, nCount );
    m_nPosition += nCount;
}

void* HcpNI::Memcpy(void* pDest, const void* pSource, const hcp_Size_t Length, void* pContext) {
	void* result = NULL;

	if (pContext != NULL) {
		auto hni = static_cast<HcpNI*>(pContext);
		result = hni->Memcpy(pDest, pSource, Length);
	}

	return result;
}

/* dv7prm... applies reverse permutation to vector.  */
void F77_NAME(dv7prm)(int *n, const int ip[], double x[])
{
    /* permute x so that x[ip[i]] := x[i]. */
    int i, nn = *n;
    double *xcp = Calloc(nn, double);

    for (i = 0; i < nn; i++) xcp[ip[i] - 1] = x[i]; /* ip contains 1-based indices */
    Memcpy(x, xcp, nn);
    Free(xcp);
}

int PublicKeyCopy(VbPublicKey *dest, const VbPublicKey *src)
{
	if (dest->key_size < src->key_size)
		return 1;

	dest->key_size = src->key_size;
	dest->algorithm = src->algorithm;
	dest->key_version = src->key_version;
	Memcpy(GetPublicKeyData(dest), GetPublicKeyDataC(src), src->key_size);
	return 0;
}

/* coerce a vector to REAL and copy the result to freshly R_alloc'd memory */
static void *RallocedREAL(SEXP x)
{
    SEXP rx = PROTECT(coerceVector(x, REALSXP));
    int lx = LENGTH(rx);
    /* We over-allocate the memory chunk so that it is never NULL. */
    /* The CHOLMOD code checks for a NULL pointer even in the length-0 case. */
    double *ans = Memcpy((double*)R_alloc(lx + 1, sizeof(double)),
			 REAL(rx), lx);
    UNPROTECT(1);
    return (void*)ans;
}

SEXP dgCMatrix_matrix_solve(SEXP Ap, SEXP b, SEXP give_sparse)
// FIXME:  add  'keep_dimnames' as argument
{
    Rboolean sparse = asLogical(give_sparse);
    if(sparse) {
	// FIXME: implement this
	error(_("dgCMatrix_matrix_solve(.., sparse=TRUE) not yet implemented"));

	/* Idea: in the  for(j = 0; j < nrhs ..) loop below, build the *sparse* result matrix
	 * ----- *column* wise -- which is perfect for dgCMatrix
	 * --> build (i,p,x) slots "increasingly" [well, allocate in batches ..]
	 *
	 * --> maybe first a protoype in R
	 */

    }
    SEXP ans = PROTECT(dup_mMatrix_as_dgeMatrix(b)),
	lu, qslot;
    CSP L, U;
    int *bdims = INTEGER(GET_SLOT(ans, Matrix_DimSym)), *p, *q;
    int j, n = bdims[0], nrhs = bdims[1];
    double *x, *ax = REAL(GET_SLOT(ans, Matrix_xSym));
    C_or_Alloca_TO(x, n, double);

    if (isNull(lu = get_factors(Ap, "LU"))) {
	install_lu(Ap, /* order = */ 1, /* tol = */ 1.0, /* err_sing = */ TRUE,
		   /* keep_dimnames = */ TRUE);
	lu = get_factors(Ap, "LU");
    }
    qslot = GET_SLOT(lu, install("q"));
    L = AS_CSP__(GET_SLOT(lu, install("L")));
    U = AS_CSP__(GET_SLOT(lu, install("U")));
    R_CheckStack();
    if (U->n != n)
	error(_("Dimensions of system to be solved are inconsistent"));
    if(nrhs >= 1 && n >= 1) {
	p = INTEGER(GET_SLOT(lu, Matrix_pSym));
	q = LENGTH(qslot) ? INTEGER(qslot) : (int *) NULL;

	for (j = 0; j < nrhs; j++) {
	    cs_pvec(p, ax + j * n, x, n);  /* x = b(p) */
	    cs_lsolve(L, x);	       /* x = L\x */
	    cs_usolve(U, x);	       /* x = U\x */
	    if (q)		       /* r(q) = x , hence
					  r = Q' U{^-1} L{^-1} P b = A^{-1} b */
		cs_ipvec(q, x, ax + j * n, n);
	    else
		Memcpy(ax + j * n, x, n);
	}
    }
    if(n >= SMALL_4_Alloca) Free(x);
    UNPROTECT(1);
    return ans;
}

/**
 * Establish a fill-reducing permutation for the sparse symmetric
 * matrix of order n represented by the column pointers Tp and row
 * indices Ti.
 *
 * @param n  order of the sparse symmetric matrix
 * @param Tp  column pointers (total length n + 1)
 * @param Ti  row indices (total length Tp[n])
 * @param Perm array of length n to hold the permutation
 * @param iPerm array of length n to hold the inverse permutation
 *
 */
void ssc_metis_order(int n, const int Tp [], const int Ti [],
                     int Perm[], int iPerm[])
{
    int  j, num_flag = 0, options_flag = 0;
    idxtype
    *perm = Calloc(n, idxtype), /* in case idxtype != int */
     *iperm = Calloc(n, idxtype),
      *xadj = Calloc(n+1, idxtype),
       *adj = Calloc(2 * (Tp[n] - n), idxtype);

    /* check row indices for correct range */
    for (j = 0; j < Tp[n]; j++)
        if (Ti[j] < 0 || Ti[j] >= n)
            error(_("row index Ti[%d] = %d is out of range [0,%d]"),
                  j, Ti[j], n - 1);
    /* temporarily use perm to store lengths */
    AZERO(perm, n);
    for (j = 0; j < n; j++) {
        int ip, p2 = Tp[j+1];
        for (ip = Tp[j]; ip < p2; ip++) {
            int i = Ti[ip];
            if (i != j) {
                perm[i]++;
                perm[j]++;
            }
        }
    }
    xadj[0] = 0;
    for (j = 0; j < n; j++) xadj[j+1] = xadj[j] + perm[j];
    /* temporarily use perm to store pointers */
    Memcpy(perm, xadj, n);
    for (j = 0; j < n; j++) {
        int ip, p2 = Tp[j+1];
        for (ip = Tp[j]; ip < p2; ip++) {
            int i = Ti[ip];
            if (i != j) {
                adj[perm[i]] = j;
                adj[perm[j]] = i;
                perm[i]++;
                perm[j]++;
            }
        }
    }
    METIS_NodeND(&n, xadj, adj, &num_flag, &options_flag, perm, iperm);
    for (j = 0; j < n; j++) {
        Perm[j] = (int) perm[j];
        iPerm[j] = (int) iperm[j];
    }
    Free(iperm);
    Free(perm);
    Free(xadj);
    Free(adj);
}