static void CG_LCP (int m, int nb, dRealMutablePtr J, int *jb, dxBody * const *body,
	dRealPtr invI, dRealMutablePtr lambda, dRealMutablePtr fc, dRealMutablePtr b,
	dRealMutablePtr lo, dRealMutablePtr hi, dRealPtr cfm, int *findex,
	dxQuickStepParameters *qs)
{
	int i,j;
	const int num_iterations = qs->num_iterations;

	// precompute iMJ = inv(M)*J'
	dRealAllocaArray (iMJ,m*12);
	compute_invM_JT (m,J,iMJ,jb,body,invI);

	dReal last_rho = 0;
	dRealAllocaArray (r,m);
	dRealAllocaArray (z,m);
	dRealAllocaArray (p,m);
	dRealAllocaArray (q,m);

	// precompute 1 / diagonals of A
	dRealAllocaArray (Ad,m);
	dRealPtr iMJ_ptr = iMJ;
	dRealPtr J_ptr = J;
	for (i=0; i<m; i++) {
		dReal sum = 0;
		for (j=0; j<6; j++) sum += iMJ_ptr[j] * J_ptr[j];
		if (jb[i*2+1] >= 0) {
			for (j=6; j<12; j++) sum += iMJ_ptr[j] * J_ptr[j];
		}
		iMJ_ptr += 12;
		J_ptr += 12;
		Ad[i] = REAL(1.0) / (sum + cfm[i]);
	}

#ifdef WARM_STARTING
	// compute residual r = b - A*lambda
	multiply_J_invM_JT (m,nb,J,iMJ,jb,cfm,fc,lambda,r);
	for (i=0; i<m; i++) r[i] = b[i] - r[i];
#else
	dSetZero (lambda,m);
	memcpy (r,b,m*sizeof(dReal));		// residual r = b - A*lambda
#endif

	for (int iteration=0; iteration < num_iterations; iteration++) {
		for (i=0; i<m; i++) z[i] = r[i]*Ad[i];	// z = inv(M)*r
		dReal rho = dot (m,r,z);		// rho = r'*z

		// @@@
		// we must check for convergence, otherwise rho will go to 0 if
		// we get an exact solution, which will introduce NaNs into the equations.
		if (rho < 1e-10) {
			printf ("CG returned at iteration %d\n",iteration);
			break;
		}

		if (iteration==0) {
			memcpy (p,z,m*sizeof(dReal));	// p = z
		}
		else {
			add (m,p,z,p,rho/last_rho);	// p = z + (rho/last_rho)*p
		}

		// compute q = (J*inv(M)*J')*p
		multiply_J_invM_JT (m,nb,J,iMJ,jb,cfm,fc,p,q);

		dReal alpha = rho/dot (m,p,q);		// alpha = rho/(p'*q)
		add (m,lambda,lambda,p,alpha);		// lambda = lambda + alpha*p
		add (m,r,r,q,-alpha);			// r = r - alpha*q
		last_rho = rho;
	}

	// compute fc = inv(M)*J'*lambda
	multiply_invM_JT (m,nb,iMJ,jb,lambda,fc);

#if 0
	// measure solution error
	multiply_J_invM_JT (m,nb,J,iMJ,jb,cfm,fc,lambda,r);
	dReal error = 0;
	for (i=0; i<m; i++) error += dFabs(r[i] - b[i]);
	printf ("lambda error = %10.6e\n",error);
#endif
}

//.........这里部分代码省略.........
		error("Error returned from miget_data_type.\n");
	}
	/* append to return list ... */
	list_index++;
	SET_VECTOR_ELT(rtnList, list_index, ScalarInteger(volume_type));
	SET_STRING_ELT(listNames, list_index, mkChar("volumeDataType"));


	/* retrieve the volume space type (talairach, native, etc) */
	result = miget_space_name(minc_volume, &space_type);
	if ( result == MI_NOERROR ) { error("Error returned from miget_space_name.\n"); }
	/* append to return list ... */
	list_index++;
	SET_VECTOR_ELT(rtnList, list_index, mkString(space_type));
	SET_STRING_ELT(listNames, list_index, mkChar("spaceType"));


	/* retrieve the total number of dimensions in this volume */
	if ( miget_volume_dimension_count(minc_volume, MI_DIMCLASS_ANY, MI_DIMATTR_ALL, &n_dimensions) != MI_NOERROR ){
		error("Error returned from miget_volume_dimension_count.\n");
	}
	/* append to return list ... */
	list_index++;
	SET_VECTOR_ELT(rtnList, list_index, ScalarInteger(n_dimensions));
	SET_STRING_ELT(listNames, list_index, mkChar("nDimensions"));


	/* load up dimension-related information */
	//
	/* first allocate the R variables */
	PROTECT( xDimSizes=allocVector(INTSXP,n_dimensions) );
	PROTECT( xDimNames=allocVector(STRSXP,n_dimensions) );
	PROTECT( xDimUnits=allocVector(STRSXP,n_dimensions) );
	PROTECT( xDimStarts=allocVector(REALSXP,n_dimensions) );
	PROTECT( xDimSteps=allocVector(REALSXP,n_dimensions) );
	n_protects = n_protects +5;

	/* next, load up the midimension struct for all dimensions*/
	dimensions = (midimhandle_t *) malloc( sizeof( midimhandle_t ) * n_dimensions );
	result = miget_volume_dimensions(minc_volume, MI_DIMCLASS_ANY, MI_DIMATTR_ALL, MI_DIMORDER_APPARENT, n_dimensions, dimensions);
	// need to check against MI_ERROR, as "result" will contain nDimensions if OK
	if ( result == MI_ERROR ) { error("Error code(%d) returned from miget_volume_dimensions.\n", result); }


	/* get the dimension sizes for all dimensions */
	result = miget_dimension_sizes(dimensions, n_dimensions, dim_sizes);
	if ( result != MI_NOERROR ) { error("Error returned from miget_dimension_sizes.\n"); }
	/* add to R vector ... */
	for (i=0; i<n_dimensions; ++i){
		INTEGER(xDimSizes)[i] = dim_sizes[i];
	}
	list_index++;
	SET_VECTOR_ELT(rtnList, list_index, xDimSizes);
	SET_STRING_ELT(listNames, list_index, mkChar("dimSizes"));


	/* get the dimension START values for all dimensions */
	result = miget_dimension_starts(dimensions, MI_ORDER_FILE, n_dimensions, dim_starts);
	if ( result == MI_ERROR ) { error("Error returned from miget_dimension_starts.\n"); }
	/* add to R vector ... */
	for (i=0; i<n_dimensions; ++i){
		REAL(xDimStarts)[i] = dim_starts[i];
	}
	list_index++;
	SET_VECTOR_ELT(rtnList, list_index, xDimStarts);
	SET_STRING_ELT(listNames, list_index, mkChar("dimStarts"));

void fft6(t_fft *x, int n1, int N1, int r, t_fft *t, int dir)
{
  int i = r+n1;
  real *x0 = x[i]; i+=N1;
  real *x1 = x[i]; i+=N1;
  real *x2 = x[i]; i+=N1;
  real *x3 = x[i]; i+=N1;
  real *x4 = x[i]; i+=N1;
  real *x5 = x[i];

  real za00 = x2[0] + x4[0];
  real za01 = x2[1] + x4[1];
  
  real za10 = x0[0] - T601*za00;
  real za11 = x0[1] - T601*za01;

  real za20;
  real za21;
  if(dir==1) {
    za20 = T600*(x4[0] - x2[0]);
    za21 = T600*(x4[1] - x2[1]);
  } else {
    za20 = T600*(x2[0] - x4[0]);
    za21 = T600*(x2[1] - x4[1]);
  }

  real a00 = x0[0] + za00;
  real a01 = x0[1] + za01;

  real a10 = za10 - za21;
  real a11 = za11 + za20;

  real a20 = za10 + za21;
  real a21 = za11 - za20;

  real zb00 = x1[0] + x5[0];
  real zb01 = x1[1] + x5[1];
  
  real zb10 = x3[0] - T601*zb00;
  real zb11 = x3[1] - T601*zb01;

  real zb20;
  real zb21;
  if(dir==1) {
    zb20 = T600*(x1[0] - x5[0]);
    zb21 = T600*(x1[1] - x5[1]);
  } else {
    zb20 = T600*(x5[0] - x1[0]);
    zb21 = T600*(x5[1] - x1[1]);
  }

  real b00 = x3[0] + zb00;
  real b01 = x3[1] + zb01;

  real b10 = zb10 - zb21;
  real b11 = zb11 + zb20;

  real b20 = zb10 + zb21;
  real b21 = zb11 - zb20;

  x0[0] = a00 + b00;
  x0[1] = a01 + b01;

  real y10 = a10 - b10;
  real y11 = a11 - b11;
  real *t1 = t[n1];
  real t10 = t1[0];
  real t11 = t1[1];
  x1[0] = REAL(t10,t11,y10,y11);
  x1[1] = IMAG(t10,t11,y10,y11);

  real y20 = a20 + b20;
  real y21 = a21 + b21;
  real *t2 = t[n1<<1];
  real t20 = t2[0];
  real t21 = t2[1];
  x2[0] = REAL(t20,t21,y20,y21);
  x2[1] = IMAG(t20,t21,y20,y21);

  real y30 = a00 - b00;
  real y31 = a01 - b01;
  real *t3 = t[n1*3];
  real t30 = t3[0];
  real t31 = t3[1];
  x3[0] = REAL(t30,t31,y30,y31);
  x3[1] = IMAG(t30,t31,y30,y31);

  real y40 = a10 + b10;
  real y41 = a11 + b11;
  real *t4 = t[n1<<2];
  real t40 = t4[0];
  real t41 = t4[1];
  x4[0] = REAL(t40,t41,y40,y41);
  x4[1] = IMAG(t40,t41,y40,y41);

  real y50 = a20 - b20;
  real y51 = a21 - b21;
  real *t5 = t[n1*5];
  real t50 = t5[0];
  real t51 = t5[1];
//.........这里部分代码省略.........

int dCollideRTL(dxGeom* g1, dxGeom* RayGeom, int Flags, dContactGeom* Contacts, int Stride){
	dIASSERT (Stride >= (int)sizeof(dContactGeom));
	dIASSERT (g1->type == dTriMeshClass);
	dIASSERT (RayGeom->type == dRayClass);
	dIASSERT ((Flags & NUMC_MASK) >= 1);

	dxTriMesh* TriMesh = (dxTriMesh*)g1;

	const dVector3& TLPosition = *(const dVector3*)dGeomGetPosition(TriMesh);
	const dMatrix3& TLRotation = *(const dMatrix3*)dGeomGetRotation(TriMesh);

	const unsigned uiTLSKind = TriMesh->getParentSpaceTLSKind();
	dIASSERT(uiTLSKind == RayGeom->getParentSpaceTLSKind()); // The colliding spaces must use matching cleanup method
	TrimeshCollidersCache *pccColliderCache = GetTrimeshCollidersCache(uiTLSKind);
	RayCollider& Collider = pccColliderCache->_RayCollider;

	dReal Length = dGeomRayGetLength(RayGeom);

	int FirstContact, BackfaceCull;
	dGeomRayGetParams(RayGeom, &FirstContact, &BackfaceCull);
	int ClosestHit = dGeomRayGetClosestHit(RayGeom);

	Collider.SetFirstContact(FirstContact != 0);
	Collider.SetClosestHit(ClosestHit != 0);
	Collider.SetCulling(BackfaceCull != 0);
	Collider.SetMaxDist(Length);

	dVector3 Origin, Direction;
	dGeomRayGet(RayGeom, Origin, Direction);

	/* Make Ray */
	Ray WorldRay;
	WorldRay.mOrig.x = Origin[0];
	WorldRay.mOrig.y = Origin[1];
	WorldRay.mOrig.z = Origin[2];
	WorldRay.mDir.x = Direction[0];
	WorldRay.mDir.y = Direction[1];
	WorldRay.mDir.z = Direction[2];

	/* Intersect */
	Matrix4x4 amatrix;
        int TriCount = 0;
        if (Collider.Collide(WorldRay, TriMesh->Data->BVTree, &MakeMatrix(TLPosition, TLRotation, amatrix))) {
                TriCount = pccColliderCache->Faces.GetNbFaces();
        }

        if (TriCount == 0) {
                return 0;
        }
	
	const CollisionFace* Faces = pccColliderCache->Faces.GetFaces();

	int OutTriCount = 0;
	for (int i = 0; i < TriCount; i++) {
		if (TriMesh->RayCallback == null ||
                    TriMesh->RayCallback(TriMesh, RayGeom, Faces[i].mFaceID,
                                         Faces[i].mU, Faces[i].mV)) {
			const int& TriIndex = Faces[i].mFaceID;
			if (!Callback(TriMesh, RayGeom, TriIndex)) {
                                continue;
                        }

			dContactGeom* Contact = SAFECONTACT(Flags, Contacts, OutTriCount, Stride);

			dVector3 dv[3];
			FetchTriangle(TriMesh, TriIndex, TLPosition, TLRotation, dv);

			dVector3 vu;
			vu[0] = dv[1][0] - dv[0][0];
			vu[1] = dv[1][1] - dv[0][1];
			vu[2] = dv[1][2] - dv[0][2];
			vu[3] = REAL(0.0);
				
			dVector3 vv;
			vv[0] = dv[2][0] - dv[0][0];
			vv[1] = dv[2][1] - dv[0][1];
			vv[2] = dv[2][2] - dv[0][2];
			vv[3] = REAL(0.0);

			dCalcVectorCross3(Contact->normal, vv, vu);	// Reversed

			// Even though all triangles might be initially valid, 
			// a triangle may degenerate into a segment after applying 
			// space transformation.
			if (dSafeNormalize3(Contact->normal))
			{
				// No sense to save on single type conversion in algorithm of this size.
				// If there would be a custom typedef for distance type it could be used 
				// instead of dReal. However using float directly is the loss of abstraction 
				// and possible loss of precision in future.
				/*float*/ dReal T = Faces[i].mDistance;
				Contact->pos[0] = Origin[0] + (Direction[0] * T);
				Contact->pos[1] = Origin[1] + (Direction[1] * T);
				Contact->pos[2] = Origin[2] + (Direction[2] * T);
				Contact->pos[3] = REAL(0.0);

				Contact->depth = T;
				Contact->g1 = TriMesh;
				Contact->g2 = RayGeom;
				Contact->side1 = TriIndex;
//.........这里部分代码省略.........

SEXP calc_het(SEXP C, SEXP LG, SEXP n, SEXP ncol, SEXP nrow)
/*********************************************************************** 
 * Functie om de heterogeniteit van een 3x3 focal area (neighbourhood)
 * in de LG kaart te berekenen.
 *
 * INVOER:
 *
 * C    = celnummers van geselecteerde rastercellen
 * LG   = LG kaart (geselecteerde cellen) als vector (de LG kaart mag 
 *        alleen de codes 1 t/m 5 of een NA bevatten!)
 * n    = lengte van 'C'
 * ncol = aantal kolommen in LG-kaart
 * nrow = aantal rijen in LG-kaart
 *
 * UITVOER:
 *
 * Een vector met lengte n en van type INTEGER.
 *
 ***********************************************************************/
{
    double *tmp;    // tijdelijke vector voor volle LG kaart
    double code;    // code op basis van LG waarden in focal area
    int i, k;       // iteratoren
    int index;
    double *xlg = REAL(LG);
    int *xc = INTEGER(C);
    int *xn = INTEGER(n), *xncol = INTEGER(ncol), *xnrow = INTEGER(nrow);
    int m = *xncol * *xnrow;
    int focal[9] = {-(*xncol+1), -*xncol, -(*xncol-1), -1, 0, 1, *xncol-1, *xncol, *xncol+1};

    tmp = Calloc(m,double);

    for (i = 0; i < *xn; i++) 
    {
        if (!ISNA(xlg[i]))
        {
            tmp[xc[i]] = xlg[i];
        }
    }

    SEXP HET;
    PROTECT(HET=allocVector(INTSXP,*xn));
    int *het = INTEGER(HET);

    for (i = 0; i < *xn; i++) 
    {
        if (!ISNA(xlg[i])) 
        {
            code = 0;
            for (k = 0; k < 9; k++)
            {
                index = xc[i] + focal[k];
                if (!(index < 1 || index > m))
                {
                    code += pow(10.0, tmp[index-1]);
                }
            }
            het[i] = code2het(code);
        } 
        else 
        {
            het[i] = NA_INTEGER;
        }
    }

    Free(tmp);
    UNPROTECT(1);
    return(HET);
}

/**
 * Define variables and write data
 */
SEXP R_write(SEXP R_filename,
             SEXP R_group,
             SEXP R_groupname,
             SEXP R_nvars,          // number of vars
             SEXP R_varname_list,   // var names
             SEXP R_var_list,       // var values
             SEXP R_varlength_list, // length of var values
             SEXP R_ndim,           // number of dims
             SEXP R_type, 
             SEXP R_comm,
             SEXP R_p,
             SEXP R_adios_rank)
{
    const char *filename = CHARPT(R_filename, 0); 
    int64_t m_adios_group = (int64_t)(REAL(R_group)[0]);
    const char *groupname = CHARPT(R_groupname, 0);
    int nvars = asInteger(R_nvars);
    MPI_Comm comm = MPI_Comm_f2c(INTEGER(R_comm)[0]);
    int size = asInteger(R_p);
    int rank = asInteger(R_adios_rank);
    
    int i, j;
    int Global_bounds, Offsets; 
    uint64_t adios_groupsize, adios_totalsize;
    int64_t m_adios_file;

    // variable to store the value converted from integer
    char str[256];

    // Define variables
    for(i = 0; i < nvars; i++) {
        const char *varname = CHAR(asChar(VECTOR_ELT(R_varname_list,i)));
        int *length = INTEGER(VECTOR_ELT(R_varlength_list, i));
        int *vndim = INTEGER(VECTOR_ELT(R_ndim, i));
        int *typetag = INTEGER(VECTOR_ELT(R_type, i));

        if((length[0] == 1) && (vndim[0] == 1)){
            // scalar
            if(typetag[0] == 0) {
                adios_define_var (m_adios_group, varname, "", adios_integer, 0, 0, 0);
            }else {
                adios_define_var (m_adios_group, varname, "", adios_double, 0, 0, 0);
            }
        }else {
            // define dimensions, global_dimensions, local_offsets and the variable
            int temp_var_length = strlen(varname) + 8;
            char* local_var = (char*)malloc(vndim[0]*temp_var_length);
            char* global_var = (char*)malloc(vndim[0]*temp_var_length);
            char* offset_var = (char*)malloc(vndim[0]*temp_var_length);

            // initialize char variables
            strcpy(local_var, "");
            strcpy(global_var, "");
            strcpy(offset_var, "");

            // j = 0
            j = 0;
            sprintf(str, "%d", j);

            char* local = (char*)malloc(temp_var_length);
            strcpy(local, varname);
            strcat(local, "_nx_");
            strcat(local, str);
            strcat(local_var, local);

            char* global = (char*)malloc(temp_var_length);
            strcpy(global, varname);
            strcat(global, "_gx_");
            strcat(global, str);
            strcat(global_var, global);

            char* offset = (char*)malloc(temp_var_length);
            strcpy(offset, varname);
            strcat(offset, "_off_");
            strcat(offset, str);
            strcat(offset_var, offset);

            // define local dim, global dim and offset for each dimension
            adios_define_var (m_adios_group, local,
                          "", adios_integer, 0, 0, 0);
            adios_define_var (m_adios_group, global,
                          "", adios_integer, 0, 0, 0);
            adios_define_var (m_adios_group, offset,
                          "", adios_integer, 0, 0, 0);

            Free(local);
            Free(global);
            Free(offset);

            for(j = 1; j < vndim[0]; j++) {
                sprintf(str, "%d", j);

                strcat(local_var, ",");
                char* local = (char*)malloc(temp_var_length);
                strcpy(local, varname);
                strcat(local, "_nx_");
                strcat(local, str);
//.........这里部分代码省略.........

 static inline double * Cast( SEXP Rvec ) { assert(Rvec); return REAL(Rvec); }

static void klu_l_demo (Long n, Long *Ap, Long *Ai, double *Ax, Long isreal)
{
    double rnorm ;
    klu_l_common Common ;
    double *B, *X, *R ;
    Long i, lunz ;

    printf ("KLU: %s, version: %d.%d.%d\n", KLU_DATE, KLU_MAIN_VERSION,
        KLU_SUB_VERSION, KLU_SUBSUB_VERSION) ;

    /* ---------------------------------------------------------------------- */
    /* set defaults */
    /* ---------------------------------------------------------------------- */

    klu_l_defaults (&Common) ;

    /* ---------------------------------------------------------------------- */
    /* create a right-hand-side */
    /* ---------------------------------------------------------------------- */

    if (isreal)
    {
        /* B = 1 + (1:n)/n */
        B = klu_l_malloc (n, sizeof (double), &Common) ;
        X = klu_l_malloc (n, sizeof (double), &Common) ;
        R = klu_l_malloc (n, sizeof (double), &Common) ;
        if (B)
        {
            for (i = 0 ; i < n ; i++)
            {
                B [i] = 1 + ((double) i+1) / ((double) n) ;
            }
        }
    }
    else
    {
        /* real (B) = 1 + (1:n)/n, imag(B) = (n:-1:1)/n */
        B = klu_l_malloc (n, 2 * sizeof (double), &Common) ;
        X = klu_l_malloc (n, 2 * sizeof (double), &Common) ;
        R = klu_l_malloc (n, 2 * sizeof (double), &Common) ;
        if (B)
        {
            for (i = 0 ; i < n ; i++)
            {
                REAL (B, i) = 1 + ((double) i+1) / ((double) n) ;
                IMAG (B, i) = ((double) n-i) / ((double) n) ;
            }
        }
    }

    /* ---------------------------------------------------------------------- */
    /* X = A\b using KLU and print statistics */
    /* ---------------------------------------------------------------------- */

    if (!klu_l_backslash (n, Ap, Ai, Ax, isreal, B, X, R, &lunz, &rnorm,
        &Common))
    {
        printf ("KLU failed\n") ;
    }
    else
    {
        printf ("n %ld nnz(A) %ld nnz(L+U+F) %ld resid %g\n"
            "recip growth %g condest %g rcond %g flops %g\n",
            n, Ap [n], lunz, rnorm, Common.rgrowth, Common.condest,
            Common.rcond, Common.flops) ;
    }

    /* ---------------------------------------------------------------------- */
    /* free the problem */
    /* ---------------------------------------------------------------------- */

    if (isreal)
    {
        klu_l_free (B, n, sizeof (double), &Common) ;
        klu_l_free (X, n, sizeof (double), &Common) ;
        klu_l_free (R, n, sizeof (double), &Common) ;
    }
    else
    {
        klu_l_free (B, 2*n, sizeof (double), &Common) ;
        klu_l_free (X, 2*n, sizeof (double), &Common) ;
        klu_l_free (R, 2*n, sizeof (double), &Common) ;
    }
    printf ("peak memory usage: %g bytes\n\n", (double) (Common.mempeak)) ;
}

//.........这里部分代码省略.........
        for (j = 0 ; j < n ; j++)
        {
            asum = 0 ;
            for (p = Ap [j] ; p < Ap [j+1] ; p++)
            {
                /* R (i) -= A (i,j) * X (j) */
                R [Ai [p]] -= Ax [p] * X [j] ;
                asum += fabs (Ax [p]) ;
            }
            anorm = MAX (anorm, asum) ;
        }
        *rnorm = 0 ;
        for (i = 0 ; i < n ; i++)
        {
            *rnorm = MAX (*rnorm, fabs (R [i])) ;
        }

        /* ------------------------------------------------------------------ */
        /* free numeric factorization */
        /* ------------------------------------------------------------------ */

        klu_l_free_numeric (&Numeric, Common) ;

    }
    else
    {

        /* ------------------------------------------------------------------ */
        /* statistics (not required to solve Ax=b) */
        /* ------------------------------------------------------------------ */

        Numeric = klu_zl_factor (Ap, Ai, Ax, Symbolic, Common) ;
        if (!Numeric)
        {
            klu_l_free_symbolic (&Symbolic, Common) ;
            return (0) ;
        }

        /* ------------------------------------------------------------------ */
        /* statistics */
        /* ------------------------------------------------------------------ */

        klu_zl_rgrowth (Ap, Ai, Ax, Symbolic, Numeric, Common) ;
        klu_zl_condest (Ap, Ax, Symbolic, Numeric, Common) ;
        klu_zl_rcond (Symbolic, Numeric, Common) ;
        klu_zl_flops (Symbolic, Numeric, Common) ;
        *lunz = Numeric->lnz + Numeric->unz - n + 
            ((Numeric->Offp) ? (Numeric->Offp [n]) : 0) ;

        /* ------------------------------------------------------------------ */
        /* solve Ax=b */
        /* ------------------------------------------------------------------ */

        for (i = 0 ; i < 2*n ; i++)
        {
            X [i] = B [i] ;
        }
        klu_zl_solve (Symbolic, Numeric, n, 1, X, Common) ;

        /* ------------------------------------------------------------------ */
        /* compute residual, rnorm = norm(b-Ax,1) / norm(A,1) */
        /* ------------------------------------------------------------------ */

        for (i = 0 ; i < 2*n ; i++)
        {
            R [i] = B [i] ;
        }
        for (j = 0 ; j < n ; j++)
        {
            asum = 0 ;
            for (p = Ap [j] ; p < Ap [j+1] ; p++)
            {
                /* R (i) -= A (i,j) * X (j) */
                i = Ai [p] ;
                REAL (R,i) -= REAL(Ax,p) * REAL(X,j) - IMAG(Ax,p) * IMAG(X,j) ;
                IMAG (R,i) -= IMAG(Ax,p) * REAL(X,j) + REAL(Ax,p) * IMAG(X,j) ;
                asum += CABS (Ax, p) ;
            }
            anorm = MAX (anorm, asum) ;
        }
        *rnorm = 0 ;
        for (i = 0 ; i < n ; i++)
        {
            *rnorm = MAX (*rnorm, CABS (R, i)) ;
        }

        /* ------------------------------------------------------------------ */
        /* free numeric factorization */
        /* ------------------------------------------------------------------ */

        klu_zl_free_numeric (&Numeric, Common) ;
    }

    /* ---------------------------------------------------------------------- */
    /* free symbolic analysis, and residual */
    /* ---------------------------------------------------------------------- */

    klu_l_free_symbolic (&Symbolic, Common) ;
    return (1) ;
}

//.........这里部分代码省略.........
    {
	/* called from registerOne */
	pDevDesc dev;
	GPar *ddp;
	sd = dd->gesd[baseRegisterIndex];
	dev = dd->dev;
	bss = sd->systemSpecific = malloc(sizeof(baseSystemState));
        /* Bail out if necessary */
        if (!bss)
            return result;
	ddp = &(bss->dp);
	GInit(ddp);
	/* For some things, the device sets the starting value at least.
	 */
	ddp->ps = dev->startps;
	ddp->col = ddp->fg = dev->startcol;
	ddp->bg = dev->startfill;
	ddp->font = dev->startfont;
	ddp->lty = dev->startlty;
	ddp->gamma = dev->startgamma;
	/* Initialise the gp settings too: formerly in addDevice. */
	copyGPar(ddp, &(bss->gp));
	GReset(dd);
	/*
	 * The device has not yet received any base output
	 */
	bss->baseDevice = FALSE;
        /* Indicate success */
        result = R_BlankString;
	break;
    }
    case GE_CopyState:
    {
	/* called from GEcopyDisplayList */
	pGEDevDesc curdd = GEcurrentDevice();
	bss = dd->gesd[baseRegisterIndex]->systemSpecific;
	bss2 = curdd->gesd[baseRegisterIndex]->systemSpecific;
	copyGPar(&(bss->dpSaved), &(bss2->dpSaved));
	restoredpSaved(curdd);
	copyGPar(&(bss2->dp), &(bss2->gp));
	GReset(curdd);
	break;
    }
    case GE_SaveState:
	/* called from GEinitDisplayList */
	bss = dd->gesd[baseRegisterIndex]->systemSpecific;
	copyGPar(&(bss->dp), &(bss->dpSaved));
	break;
    case GE_RestoreState:
	/* called from GEplayDisplayList */
	bss = dd->gesd[baseRegisterIndex]->systemSpecific;
	restoredpSaved(dd);
	copyGPar(&(bss->dp), &(bss->gp));
	GReset(dd);
	break;
    case GE_SaveSnapshotState:
	/* called from GEcreateSnapshot */
	bss = dd->gesd[baseRegisterIndex]->systemSpecific;
	/* Changed from INTSXP in 2.7.0: but saved graphics lists
	   are protected by an R version number */
	PROTECT(result = allocVector(RAWSXP, sizeof(GPar)));
	copyGPar(&(bss->dpSaved), (GPar*) RAW(result));
	UNPROTECT(1);
	break;
    case GE_RestoreSnapshotState:
	/* called from GEplaySnapshot */
	bss = dd->gesd[baseRegisterIndex]->systemSpecific;
	copyGPar((GPar*) RAW(data), &(bss->dpSaved));
	restoredpSaved(dd);
	copyGPar(&(bss->dp), &(bss->gp));
	GReset(dd);
	break;
    case GE_CheckPlot:
	/* called from GEcheckState:
	   Check that the current plotting state is "valid"
	 */
	bss = dd->gesd[baseRegisterIndex]->systemSpecific;
	result = ScalarLogical(bss->baseDevice ?
			       (bss->gp.state == 1) && bss->gp.valid :
			       TRUE);
	break;
    case GE_ScalePS:
    {
	/* called from GEhandleEvent in devWindows.c */
	GPar *ddp, *ddpSaved;
	bss = dd->gesd[baseRegisterIndex]->systemSpecific;
	ddp = &(bss->dp);
	ddpSaved = &(bss->dpSaved);
	if (isReal(data) && LENGTH(data) == 1) {
	    double rf = REAL(data)[0];
	    ddp->scale *= rf;
	    /* Modify the saved settings so this effects display list too */
	    ddpSaved->scale *= rf;
	} else
	  error("event 'GE_ScalePS' requires a single numeric value");
	break;
    }
    }
    return result;
}

SEXP non_duplicates (SEXP x_, SEXP fromLast_) {
  int fromLast = asLogical(fromLast_),
      i, d=0,
      len   = length(x_);
  
  int *x_int;
  double *x_real;

  SEXP duplicates;
  int *duplicates_int;
  PROTECT(duplicates = allocVector(INTSXP, len)); /* possibly resize this */
  duplicates_int = INTEGER(duplicates);

  if(!fromLast) { /* keep first observation */
    duplicates_int[0] = ++d;
    switch(TYPEOF(x_)) {
      case INTSXP:
        x_int = INTEGER(x_);
        for(i=1; i < len-1; i++) {
          if( x_int[i-1] != x_int[i]) {
#ifdef DEBUG
            Rprintf("i=%i:  x[i-1]=%i, x[i]=%i\n",i,x_int[i-1],x_int[i]);
#endif
            duplicates_int[d++] = i+1;
          }
        }      
        break;
      case REALSXP:
        x_real = REAL(x_);
        for(i=1; i < len; i++) {
          /*
          if( x_real[i-1] == x_real[i])
            duplicates_int[d++] = (int)(-1*(i+1));
          */
          if( x_real[i-1] != x_real[i])
            duplicates_int[d++] = i+1;
        }      
        break;
      default:
        error("only numeric types supported");
        break;
    }
  } else {    /* keep last observation  */
    switch(TYPEOF(x_)) {
      case INTSXP:
        x_int = INTEGER(x_);
        for(i=1; i < len; i++) {
          if( x_int[i-1] != x_int[i])
            duplicates_int[d++] = i;
        }      
        break;
      case REALSXP:
        x_real = REAL(x_);
        for(i=1; i < len; i++) {
          if( x_real[i-1] != x_real[i])
            duplicates_int[d++] = i;
        }      
        break;
      default:
        error("only numeric types supported");
        break;
    }
    duplicates_int[d++] = len;
  }
  UNPROTECT(1);
  return(lengthgets(duplicates, d));
}

/* used in connections.c */
SEXP xlengthgets(SEXP x, R_xlen_t len)
{
    R_xlen_t lenx, i;
    SEXP rval, names, xnames, t;
    if (!isVector(x) && !isVectorizable(x))
	error(_("cannot set length of non-vector"));
    lenx = xlength(x);
    if (lenx == len)
	return (x);
    PROTECT(rval = allocVector(TYPEOF(x), len));
    PROTECT(xnames = getAttrib(x, R_NamesSymbol));
    if (xnames != R_NilValue)
	names = allocVector(STRSXP, len);
    else names = R_NilValue;	/*- just for -Wall --- should we do this ? */
    switch (TYPEOF(x)) {
    case NILSXP:
	break;
    case LGLSXP:
    case INTSXP:
	for (i = 0; i < len; i++)
	    if (i < lenx) {
		INTEGER(rval)[i] = INTEGER(x)[i];
		if (xnames != R_NilValue)
		    SET_STRING_ELT(names, i, STRING_ELT(xnames, i));
	    }
	    else
		INTEGER(rval)[i] = NA_INTEGER;
	break;
    case REALSXP:
	for (i = 0; i < len; i++)
	    if (i < lenx) {
		REAL(rval)[i] = REAL(x)[i];
		if (xnames != R_NilValue)
		    SET_STRING_ELT(names, i, STRING_ELT(xnames, i));
	    }
	    else
		REAL(rval)[i] = NA_REAL;
	break;
    case CPLXSXP:
	for (i = 0; i < len; i++)
	    if (i < lenx) {
		COMPLEX(rval)[i] = COMPLEX(x)[i];
		if (xnames != R_NilValue)
		    SET_STRING_ELT(names, i, STRING_ELT(xnames, i));
	    }
	    else {
		COMPLEX(rval)[i].r = NA_REAL;
		COMPLEX(rval)[i].i = NA_REAL;
	    }
	break;
    case STRSXP:
	for (i = 0; i < len; i++)
	    if (i < lenx) {
		SET_STRING_ELT(rval, i, STRING_ELT(x, i));
		if (xnames != R_NilValue)
		    SET_STRING_ELT(names, i, STRING_ELT(xnames, i));
	    }
	    else
		SET_STRING_ELT(rval, i, NA_STRING);
	break;
    case LISTSXP:
	for (t = rval; t != R_NilValue; t = CDR(t), x = CDR(x)) {
	    SETCAR(t, CAR(x));
	    SET_TAG(t, TAG(x));
	}
    case VECSXP:
	for (i = 0; i < len; i++)
	    if (i < lenx) {
		SET_VECTOR_ELT(rval, i, VECTOR_ELT(x, i));
		if (xnames != R_NilValue)
		    SET_STRING_ELT(names, i, STRING_ELT(xnames, i));
	    }
	break;
    case RAWSXP:
	for (i = 0; i < len; i++)
	    if (i < lenx) {
		RAW(rval)[i] = RAW(x)[i];
		if (xnames != R_NilValue)
		    SET_STRING_ELT(names, i, STRING_ELT(xnames, i));
	    }
	    else
		RAW(rval)[i] = (Rbyte) 0;
	break;
    default:
	UNIMPLEMENTED_TYPE("length<-", x);
    }
    if (isVector(x) && xnames != R_NilValue)
	setAttrib(rval, R_NamesSymbol, names);
    UNPROTECT(2);
    return rval;
}

SEXP cliques_R(SEXP net, SEXP sn, SEXP sm, SEXP stabulatebyvert, SEXP scomembership, SEXP senumerate)
/*Maximal clique enumeration as an R-callable (.Call) function.  net should be an sna edgelist (w/n vertices and m/2 edges), and must be pre-symmetrized.  stabulatebyvert should be 0 if no tabulation is to be performed, or 1 for vertex-level tabulation of clique membership.  scomembership should be 0 for no co-membership tabulation, 1 for aggregate vertex-by-vertex tabulation, and 2 for size-by-vertex-by-vertex tabulation.  Finally, senumerate should be 1 iff the enumerated clique list should be returned.  (The current algorithm enumerates them internally, regardless.  This is b/c I am lazy, and didn't fold all of the tabulation tasks into the recursion process.  Life is hard.)*/
{
  int n,tabulate,comemb,enumerate,*gotcomp,*compmemb,i,j,k,maxcsize,pc=0;
  double *ccount,*pccountvec,*pcocliquevec=NULL;
  snaNet *g;
  slelement *sep,*sep2,*k0;
  element **clist,*ep;
  SEXP smaxcsize,ccountvec,outlist,cliquevec=R_NilValue;
  SEXP temp=R_NilValue,sp=R_NilValue,cocliquevec=R_NilValue;

  /*Coerce what needs coercin'*/
  PROTECT(sn=coerceVector(sn,INTSXP)); pc++;
  PROTECT(net=coerceVector(net,REALSXP)); pc++;
  PROTECT(stabulatebyvert=coerceVector(stabulatebyvert,INTSXP)); pc++;
  PROTECT(scomembership=coerceVector(scomembership,INTSXP)); pc++;
  PROTECT(senumerate=coerceVector(senumerate,INTSXP)); pc++;
  n=INTEGER(sn)[0];
  tabulate=INTEGER(stabulatebyvert)[0];
  comemb=INTEGER(scomembership)[0];
  enumerate=INTEGER(senumerate)[0];

  /*Pre-allocate what needs pre-allocatin'*/
  ccount=(double *)R_alloc(n,sizeof(double));
  PROTECT(smaxcsize=allocVector(INTSXP,1)); pc++;
  clist=(element **)R_alloc(n,sizeof(element *));
  for(i=0;i<n;i++){
    ccount[i]=0.0;
    clist[i]=NULL;
  }
    
  /*Initialize sna internal network*/
  GetRNGstate();
  g=elMatTosnaNet(REAL(net),INTEGER(sn),INTEGER(sm));

  /*Calculate the components of g*/
  compmemb=undirComponents(g);

  /*Accumulate cliques across components*/
  gotcomp=(int *)R_alloc(compmemb[0],sizeof(int));
  for(i=0;i<compmemb[0];i++)
    gotcomp[i]=0;
  for(i=0;i<n;i++)                   /*Move through vertices in order*/
    if(!gotcomp[compmemb[i+1]-1]){   /*Take first vertex of each component*/
      gotcomp[compmemb[i+1]-1]++;              /*Mark component as visited*/
      /*Get the first maximal clique in this component*/
      k0=slistInsert(NULL,(double)i,NULL);
      k0=cliqueFirstChild(g,k0);
      /*Recursively enumerate all cliques within the component*/
      cliqueRecurse(g,k0,i,clist,ccount,compmemb);
    }
  PutRNGstate();
  
  /*Find the maximum clique size (to cut down on subsequent memory usage)*/
  INTEGER(smaxcsize)[0]=n+1;
  for(i=n-1;(i>=0)&(INTEGER(smaxcsize)[0]==n+1);i--)
    if(ccount[i]>0.0)
      INTEGER(smaxcsize)[0]=i+1;
  maxcsize=INTEGER(smaxcsize)[0];

  /*Allocate memory for R return value objects*/
  if(tabulate){
    PROTECT(ccountvec=allocVector(REALSXP,maxcsize*(1+n))); pc++;
    for(i=0;i<maxcsize*(1+n);i++)
      REAL(ccountvec)[i]=0.0;
  }else{
    PROTECT(ccountvec=allocVector(REALSXP,maxcsize)); pc++;
    for(i=0;i<maxcsize;i++)
      REAL(ccountvec)[i]=0.0;
  }
  pccountvec=REAL(ccountvec);
  switch(comemb){
    case 0:
      cocliquevec=R_NilValue;
      pcocliquevec=NULL;
      break;
    case 1:
      PROTECT(cocliquevec=allocVector(REALSXP,n*n)); pc++;
      for(i=0;i<n*n;i++)
        REAL(cocliquevec)[i]=0.0;
      pcocliquevec=REAL(cocliquevec);
      break;
    case 2:
      PROTECT(cocliquevec=allocVector(REALSXP,maxcsize*n*n)); pc++;
      for(i=0;i<maxcsize*n*n;i++)
        REAL(cocliquevec)[i]=0.0;
      pcocliquevec=REAL(cocliquevec);
      break;
  }
  if(enumerate){
    PROTECT(cliquevec=allocVector(VECSXP,maxcsize)); pc++;
    for(i=0;i<maxcsize;i++){
      if(ccount[i]==0.0)
        SET_VECTOR_ELT(cliquevec,i,R_NilValue);
      else{
        PROTECT(temp=allocVector(VECSXP,(int)(ccount[i])));
        SET_VECTOR_ELT(cliquevec,i,temp);
        UNPROTECT(1);
      }
    }
//.........这里部分代码省略.........

SEXP bicomponents_R(SEXP net, SEXP sn, SEXP sm)
{
  snaNet *g;
  int *parent,*num,*back,*dfn,i,j,n,count,pc=0;
  element *complist,*ep,*ep2,*es;
  SEXP bicomps,bcl,memb,outlist;

  /*Coerce what needs coercin'*/
  //Rprintf("Initial coercion\n");
  PROTECT(sn=coerceVector(sn,INTSXP)); pc++;
  PROTECT(sm=coerceVector(sm,INTSXP)); pc++;
  PROTECT(net=coerceVector(net,REALSXP)); pc++;
  n=INTEGER(sn)[0];

  /*Initialize sna internal network*/
  GetRNGstate();
  g=elMatTosnaNet(REAL(net),INTEGER(sn),INTEGER(sm));

  /*Calculate the sorting stat*/
  parent=(int *)R_alloc(n,sizeof(int));
  num=(int *)R_alloc(n,sizeof(int));
  back=(int *)R_alloc(n,sizeof(int));
  dfn=(int *)R_alloc(1,sizeof(int));
  for(i=0;i<n;i++){
    parent[i]=-1;
    num[i]=0;
    back[i]=n+1;
  }
  *dfn=0;

  /*Initialize the component list*/
  complist=(element *)R_alloc(1,sizeof(element));
  complist->val=0.0;
  complist->next=NULL;
  complist->dp=NULL;

  /*Walk the graph components, finding bicomponents*/
  es=(element *)R_alloc(1,sizeof(element));
  for(i=0;i<n;i++)
    if(num[i]==0){
      es->next=NULL;
      bicomponentRecurse(g,complist,es,parent,num,back,dfn,i);
    }

  /*Transfer information from complist to output vector*/
  //Rprintf("Gathering components...\n");
  count=(int)complist->val;
  PROTECT(outlist=allocVector(VECSXP,2)); pc++;
  PROTECT(bicomps=allocVector(VECSXP,count)); pc++;
  PROTECT(memb=allocVector(INTSXP,n)); pc++;
  for(i=0;i<n;i++)  /*Initialize memberships, since some have none*/
    INTEGER(memb)[i]=-1;
  ep=complist->next;
  for(i=0;i<count;i++){
    PROTECT(bcl=allocVector(INTSXP,(int)ep->val));
    j=0;
    for(ep2=(element *)ep->dp;ep2!=NULL;ep2=ep2->next){
      INTEGER(bcl)[j++]=(int)ep2->val+1;
      INTEGER(memb)[(int)ep2->val]=i+1;
    }
    SET_VECTOR_ELT(bicomps,i,bcl);
    UNPROTECT(1);
    ep=ep->next;
  }
  SET_VECTOR_ELT(outlist,0,bicomps); 
  SET_VECTOR_ELT(outlist,1,memb); 

  /*Unprotect and return*/
  PutRNGstate();
  UNPROTECT(pc);
  return outlist;
}

SEXP do_subset_xts(SEXP x, SEXP sr, SEXP sc, SEXP drop) //SEXP s, SEXP call, int drop)
{
    SEXP attr, result, dim;
    int nr, nc, nrs, ncs;
    int i, j, ii, jj, ij, iijj;
    int mode;
    int *int_x=NULL, *int_result=NULL, *int_newindex=NULL, *int_index=NULL;
    double *real_x=NULL, *real_result=NULL, *real_newindex=NULL, *real_index=NULL;

    nr = nrows(x);
    nc = ncols(x);

    if( length(x)==0 )
      return x;

    dim = getAttrib(x, R_DimSymbol);

    nrs = LENGTH(sr);
    ncs = LENGTH(sc);
    int *int_sr=NULL, *int_sc=NULL;
    int_sr = INTEGER(sr);
    int_sc = INTEGER(sc);

    mode = TYPEOF(x);

    result = allocVector(mode, nrs*ncs);
    PROTECT(result);


    if( mode==INTSXP ) {
      int_x = INTEGER(x);
      int_result = INTEGER(result);
    } else
    if( mode==REALSXP ) {
      real_x = REAL(x);
      real_result = REAL(result);
    }

    /* code to handle index of xts object efficiently */
    SEXP index, newindex;
    int indx;

    index = getAttrib(x, install("index"));
    PROTECT(index);

    if(TYPEOF(index) == INTSXP) {
      newindex = allocVector(INTSXP, LENGTH(sr));
      PROTECT(newindex);
      int_newindex = INTEGER(newindex);
      int_index = INTEGER(index);
      for(indx = 0; indx < nrs; indx++) {
        int_newindex[indx] = int_index[ (int_sr[indx])-1];
      }
      copyAttributes(index, newindex);
      setAttrib(result, install("index"), newindex);
      UNPROTECT(1);
    }
    if(TYPEOF(index) == REALSXP) {
      newindex = allocVector(REALSXP, LENGTH(sr));
      PROTECT(newindex);
      real_newindex = REAL(newindex);
      real_index = REAL(index);
      for(indx = 0; indx < nrs; indx++) {
        real_newindex[indx] = real_index[ (int_sr[indx])-1 ];
      }
      copyAttributes(index, newindex);
      setAttrib(result, install("index"), newindex);
      UNPROTECT(1);
    }

    for (i = 0; i < nrs; i++) {
      ii = int_sr[i];
      if (ii != NA_INTEGER) {
        if (ii < 1 || ii > nr)
          error("i is out of range\n");
        ii--;
      }
      /* Begin column loop */
      for (j = 0; j < ncs; j++) {
        //jj = INTEGER(sc)[j];
        jj = int_sc[j];
        if (jj != NA_INTEGER) {
        if (jj < 1 || jj > nc)
          error("j is out of range\n");
        jj--;
        }
        ij = i + j * nrs;
        if (ii == NA_INTEGER || jj == NA_INTEGER) {
          switch ( mode ) {
            case REALSXP:
                 real_result[ij] = NA_REAL;
                 break;
            case LGLSXP:
            case INTSXP:
                 int_result[ij] = NA_INTEGER;
                 break;
            case CPLXSXP:
                 COMPLEX(result)[ij].r = NA_REAL;
                 COMPLEX(result)[ij].i = NA_REAL;
                 break;
//.........这里部分代码省略.........

SEXP minc2_apply(SEXP filenames, SEXP fn, SEXP have_mask, 
		 SEXP mask, SEXP mask_value, SEXP rho) {
  int                result;
  mihandle_t         *hvol, hmask;
  int                i, v0, v1, v2, output_index, buffer_index;
  unsigned long     start[3], count[3];
  //unsigned long      location[3];
  int                num_files;
  double             *xbuffer, *xoutput, **full_buffer;
  double             *xhave_mask, *xmask_value;
  double             *mask_buffer;
  midimhandle_t      dimensions[3];
  misize_t            sizes[3];
  SEXP               output, buffer;
  //SEXP               R_fcall;
  

  /* allocate memory for volume handles */
  num_files = LENGTH(filenames);
  hvol = malloc(num_files * sizeof(mihandle_t));

  Rprintf("Number of volumes: %i\n", num_files);

  /* open the mask - if so desired */
  xhave_mask = REAL(have_mask);
  if (xhave_mask[0] == 1) {
    result = miopen_volume(CHAR(STRING_ELT(mask, 0)),
			   MI2_OPEN_READ, &hmask);
    if (result != MI_NOERROR) {
      error("Error opening mask: %s.\n", CHAR(STRING_ELT(mask, 0)));
    }
  }
  
  /* get the value inside that mask */
  xmask_value = REAL(mask_value);

  /* open each volume */
  for(i=0; i < num_files; i++) {
    result = miopen_volume(CHAR(STRING_ELT(filenames, i)),
      MI2_OPEN_READ, &hvol[i]);
    if (result != MI_NOERROR) {
      error("Error opening input file: %s.\n", CHAR(STRING_ELT(filenames,i)));
    }
  }

  /* get the file dimensions and their sizes - assume they are the same*/
  miget_volume_dimensions( hvol[0], MI_DIMCLASS_SPATIAL,
			   MI_DIMATTR_ALL, MI_DIMORDER_FILE,
			   3, dimensions);
  result = miget_dimension_sizes( dimensions, 3, sizes );
  Rprintf("Volume sizes: %i %i %i\n", sizes[0], sizes[1], sizes[2]);

  /* allocate the output buffer */
  PROTECT(output=allocVector(REALSXP, (sizes[0] * sizes[1] * sizes[2])));
  xoutput = REAL(output);

  /* allocate the local buffer that will be passed to the function */
  PROTECT(buffer=allocVector(REALSXP, num_files));
  xbuffer = REAL(buffer); 

  //PROTECT(R_fcall = lang2(fn, R_NilValue));


  /* allocate first dimension of the buffer */
  full_buffer = malloc(num_files * sizeof(double));

  /* allocate second dimension of the buffer 
     - big enough to hold one slice per subject at a time */
  for (i=0; i < num_files; i++) {
    full_buffer[i] = malloc(sizes[1] * sizes[2] * sizeof(double));
  }
  
  /* allocate buffer for mask - if necessary */
  if (xhave_mask[0] == 1) {