自学教程：C++ vmaxset函数代码示例

/* Returns: */staticenum pmatchpstrmatch(SEXP target, SEXP input, size_t slen){    const char *st = "";    const void *vmax = vmaxget();    if(target == R_NilValue)	return NO_MATCH;    switch (TYPEOF(target)) {    case SYMSXP:	st = CHAR(PRINTNAME(target));	break;    case CHARSXP:	st = translateChar(target);	break;    }    if(strncmp(st, translateChar(input), slen) == 0) {	vmaxset(vmax);	return (strlen(st) == slen) ?  EXACT_MATCH : PARTIAL_MATCH;    } else {	vmaxset(vmax);	return NO_MATCH;    }}

SEXP do_aqua_custom_print(SEXP call, SEXP op, SEXP args, SEXP env){    const void *vm;    const char *ct;    int cpr;    SEXP rv, objType, obj;    if (!ptr_Raqua_CustomPrint) return R_NilValue;    checkArity(op, args);    vm = vmaxget();    objType = CAR(args); args = CDR(args);    obj = CAR(args);    if (!isString(objType) || LENGTH(objType)<1)	errorcall(call, "invalid arguments");    ct=CHAR(STRING_ELT(objType,0));    cpr=ptr_Raqua_CustomPrint(ct, obj);    /* FIXME: trying to store a pointer in an integer is wrong */    PROTECT(rv=allocVector(INTSXP, 1));    INTEGER(rv)[0]=cpr;    vmaxset(vm);    UNPROTECT(1);    return rv;}

/*  This is the routine associated with the getNativeSymbolInfo()  function and it takes the name of a symbol and optionally an  object identifier (package usually) in which to restrict the search  for this symbol. It resolves the symbol and returns it to the caller  giving the symbol address, the package information (i.e. name and  fully qualified shared object name). If the symbol was explicitly  registered (rather than dynamically resolved by R), then we pass  back that information also, giving the number of arguments it  expects and the interface by which it should be called.  The returned object has class NativeSymbol. If the symbol was  registered, we add a class identifying the interface type  for which it is intended (i.e. .C(), .Call(), etc.) */SEXP attribute_hiddenR_getSymbolInfo(SEXP sname, SEXP spackage, SEXP withRegistrationInfo){    const void *vmax = vmaxget();    const char *package, *name;    R_RegisteredNativeSymbol symbol = {R_ANY_SYM, {NULL}, NULL};    SEXP sym = R_NilValue;    DL_FUNC f = NULL;    package = "";    name = translateChar(STRING_ELT(sname, 0));    if(length(spackage)) {	if(TYPEOF(spackage) == STRSXP)	    package = translateChar(STRING_ELT(spackage, 0));	else if(TYPEOF(spackage) == EXTPTRSXP &&		R_ExternalPtrTag(spackage) == install("DLLInfo")) {	    f = R_dlsym((DllInfo *) R_ExternalPtrAddr(spackage), name, &symbol);	    package = NULL;	} else	    error(_("must pass package name or DllInfo reference"));    }    if(package)	f = R_FindSymbol(name, package, &symbol);    if(f)	sym = createRSymbolObject(sname, f, &symbol,				  LOGICAL(withRegistrationInfo)[0]);    vmaxset(vmax);    return sym;}

static void whittle2 (Array acf, Array Aold, Array Bold, int lag,		      char *direction, Array A, Array K, Array E){    int d, i, nser=DIM(acf)[1];    const void *vmax;    Array beta, tmp, id;    d = strcmp(direction, "forward") == 0;    vmax = vmaxget();    beta = make_zero_matrix(nser,nser);    tmp = make_zero_matrix(nser, nser);    id = make_identity_matrix(nser);    set_array_to_zero(E);    copy_array(id, subarray(A,0));    for(i = 0; i < lag; i++) {       matrix_prod(subarray(acf,lag - i), subarray(Aold,i), d, 1, tmp);       array_op(beta, tmp, '+', beta);       matrix_prod(subarray(acf,i), subarray(Bold,i), d, 1, tmp);       array_op(E, tmp, '+', E);    }    qr_solve(E, beta, K);    transpose_matrix(K,K);    for (i = 1; i <= lag; i++) {	matrix_prod(K, subarray(Bold,lag - i), 0, 0, tmp);	array_op(subarray(Aold,i), tmp, '-', subarray(A,i));    }    vmaxset(vmax);}

static SEXPsummary_func_builtin (Sampler *ss, int currentPeriod, SEXP currentStreams,                      SEXP currentLogWeights){        int dspp = ss->dimSummPerPeriod, ns = ss->nStreams;        int ii, jj, start, nProtected = 0;        double *scs  = REAL(currentStreams);        double *sclw = REAL(currentLogWeights);        double *scaw = REAL(ss->SEXPCurrentAdjWeights);        double sumcaw, *summ;        SEXP SEXPSumm;        void *vmax = vmaxget( );        PROTECT(SEXPSumm = allocVector(REALSXP, dspp)); ++nProtected;        summ = REAL(SEXPSumm);                /*         * Note: here dimSummPerPeriod == dimPerPeriod and we only         * provide the weighted average of each of the dimensions         */        sumcaw = sampler_adjust_log_weights(ns, sclw, scaw);        for (ii = 0; ii < dspp; ++ii) {                summ[ii] = 0.0; start = ii * ns;                for (jj = 0; jj < ns; ++jj)                        summ[ii] += scaw[jj] * scs[start + jj];                summ[ii] /= sumcaw;        }        UNPROTECT(nProtected);        vmaxset(vmax);        return SEXPSumm;}

attribute_hiddenvoid printMatrix(SEXP x, int offset, SEXP dim, int quote, int right,                 SEXP rl, SEXP cl, const char *rn, const char *cn){    /* 'rl' and 'cl' are dimnames(.)[[1]] and dimnames(.)[[2]]  whereas     * 'rn' and 'cn' are the  names(dimnames(.))     */    const void *vmax = vmaxget();    int r = INTEGER(dim)[0];    int c = INTEGER(dim)[1], r_pr;    /* PR#850 */    if ((rl != R_NilValue) && (r > length(rl)))        error(_("too few row labels"));    if ((cl != R_NilValue) && (c > length(cl)))        error(_("too few column labels"));    if (r == 0 && c == 0) {        Rprintf("<0 x 0 matrix>/n");        return;    }    r_pr = r;    if(c > 0 && R_print.max / c < r) /* avoid integer overflow */        /* using floor(), not ceil(), since 'c' could be huge: */        r_pr = R_print.max / c;    switch (TYPEOF(x)) {    case LGLSXP:        printLogicalMatrix(x, offset, r_pr, r, c, rl, cl, rn, cn);        break;    case INTSXP:        printIntegerMatrix(x, offset, r_pr, r, c, rl, cl, rn, cn);        break;    case REALSXP:        printRealMatrix	  (x, offset, r_pr, r, c, rl, cl, rn, cn);        break;    case CPLXSXP:        printComplexMatrix(x, offset, r_pr, r, c, rl, cl, rn, cn);        break;    case STRSXP:        if (quote) quote = '"';        printStringMatrix (x, offset, r_pr, r, c, quote, right, rl, cl, rn, cn);        break;    case RAWSXP:        printRawMatrix	  (x, offset, r_pr, r, c, rl, cl, rn, cn);        break;    default:        UNIMPLEMENTED_TYPE("printMatrix", x);    }#ifdef ENABLE_NLS    if(r_pr < r) // number of formats must be consistent here        Rprintf(ngettext(" [ reached getOption(/"max.print/") -- omitted %d row ]/n",                         " [ reached getOption(/"max.print/") -- omitted %d rows ]/n",                         r - r_pr),                r - r_pr);#else    if(r_pr < r)        Rprintf(" [ reached getOption(/"max.print/") -- omitted %d rows ]/n",                r - r_pr);#endif    vmaxset(vmax);}

static SEXP MakeVariable(int k, SEXP tag){    const void *vmax = vmaxget();        char buf[64];    snprintf(buf, 64, "%s%d", translateChar(STRING_ELT(tag, 0)), k);    vmaxset(vmax);    return install(buf);}

SEXP dotTclObjv(SEXP args){    SEXP t,	avec = CADR(args),	nm = getAttrib(avec, R_NamesSymbol);    int objc, i, result;    Tcl_Obj **objv;    const void *vmax = vmaxget();    for (objc = 0, i = 0; i < length(avec); i++){	if (!isNull(VECTOR_ELT(avec, i)))	    objc++;	if (!isNull(nm) && strlen(translateChar(STRING_ELT(nm, i))))	    objc++;    }    objv = (Tcl_Obj **) R_alloc(objc, sizeof(Tcl_Obj *));    for (objc = i = 0; i < length(avec); i++){	const char *s;	char *tmp;	if (!isNull(nm) && strlen(s = translateChar(STRING_ELT(nm, i)))){	    tmp = calloc(strlen(s)+2, sizeof(char));	    *tmp = '-';	    strcpy(tmp+1, s);	    objv[objc++] = Tcl_NewStringObj(tmp, -1);	    free(tmp);	}	if (!isNull(t = VECTOR_ELT(avec, i)))	    objv[objc++] = (Tcl_Obj *) R_ExternalPtrAddr(t);    }    for (i = objc; i--; ) Tcl_IncrRefCount(objv[i]);    result = Tcl_EvalObjv(RTcl_interp, objc, objv, 0);    for (i = objc; i--; ) Tcl_DecrRefCount(objv[i]);    if (result == TCL_ERROR)    {	char p[512];	if (strlen(Tcl_GetStringResult(RTcl_interp)) > 500)	    strcpy(p, _("tcl error./n"));	else {	    char *res;	    Tcl_DString  res_ds;	    Tcl_DStringInit(&res_ds);	    res = Tcl_UtfToExternalDString(NULL,					   Tcl_GetStringResult(RTcl_interp),					   -1, &res_ds);	    snprintf(p, sizeof(p), "[tcl] %s./n", res);	    Tcl_DStringFree(&res_ds);	}	error(p);    }    SEXP res = makeRTclObject(Tcl_GetObjResult(RTcl_interp));    vmaxset(vmax);    return res;}

void CRF::Decode_Tree(){	void *vmax = vmaxget(); 	TreeBP(true);	MaxOfMarginals();	vmaxset(vmax); }

static void matrix_prod(Array mat1, Array mat2, int trans1, int trans2, Array ans)/*    General matrix product between mat1 and mat2. Put answer in ans.    trans1 and trans2 are logical flags which indicate if the matrix is    to be transposed. Normal matrix multiplication has trans1 = trans2 = 0.*/{    int i,j,k,K1,K2;    const void *vmax;    double m1, m2;    Array tmp;    /* Test whether everything is a matrix */    assert(DIM_LENGTH(mat1) == 2 &&	   DIM_LENGTH(mat2) == 2 && DIM_LENGTH(ans) == 2);    /* Test whether matrices conform. K is the dimension that is       lost by multiplication */    if (trans1) {	assert ( NCOL(mat1) == NROW(ans) );	K1 = NROW(mat1);    }    else {	assert ( NROW(mat1) == NROW(ans) );	K1 = NCOL(mat1);    }    if (trans2) {	assert ( NROW(mat2) == NCOL(ans) );	K2 = NCOL(mat2);    }    else {	assert ( NCOL(mat2) == NCOL(ans) );	K2 = NROW(mat2);    }    assert (K1 == K2);    tmp = init_array();    /* In case ans is the same as mat1 or mat2, we create a temporary       matrix to hold the answer, then copy it to ans    */    vmax = vmaxget();    tmp = make_zero_matrix(NROW(ans), NCOL(ans));    for (i = 0; i < NROW(tmp); i++) {	for (j = 0; j < NCOL(tmp); j++) {	    for(k = 0; k < K1; k++) {		    m1 = (trans1) ? MATRIX(mat1)[k][i] : MATRIX(mat1)[i][k];		    m2 = (trans2) ? MATRIX(mat2)[j][k] : MATRIX(mat2)[k][j];		    MATRIX(tmp)[i][j] += m1 * m2;	    }	}    }    copy_array(tmp, ans);    vmaxset(vmax);}

FILE *RC_fopen(const SEXP fn, const char *mode, const Rboolean expand){    const void *vmax = vmaxget();    const char *filename = translateChar(fn), *res;    if(fn == NA_STRING || !filename) return NULL;    if(expand) res = R_ExpandFileName(filename);    else res = filename;    vmaxset(vmax);    return fopen(res, mode);}

// unused now from Rdouble bessel_j(double x, double alpha){    int nb, ncalc;    double na, *bj;#ifndef MATHLIB_STANDALONE    const void *vmax;#endif#ifdef IEEE_754    /* NaNs propagated correctly */    if (ISNAN(x) || ISNAN(alpha)) return x + alpha;#endif    if (x < 0) {	ML_ERROR(ME_RANGE, "bessel_j");	return ML_NAN;    }    na = floor(alpha);    if (alpha < 0) {	/* Using Abramowitz & Stegun  9.1.2	 * this may not be quite optimal (CPU and accuracy wise) */	return(((alpha - na == 0.5) ? 0 : bessel_j(x, -alpha) * cospi(alpha)) +	       ((alpha      == na ) ? 0 : bessel_y(x, -alpha) * sinpi(alpha)));    }    else if (alpha > 1e7) {	MATHLIB_WARNING("besselJ(x, nu): nu=%g too large for bessel_j() algorithm", alpha);	return ML_NAN;    }    nb = 1 + (int)na; /* nb-1 <= alpha < nb */    alpha -= (double)(nb-1);#ifdef MATHLIB_STANDALONE    bj = (double *) calloc(nb, sizeof(double));#ifndef _RENJIN    if (!bj) MATHLIB_ERROR("%s", _("bessel_j allocation error"));#endif#else    vmax = vmaxget();    bj = (double *) R_alloc((size_t) nb, sizeof(double));#endif    J_bessel(&x, &alpha, &nb, bj, &ncalc);    if(ncalc != nb) {/* error input */      if(ncalc < 0)	MATHLIB_WARNING4(_("bessel_j(%g): ncalc (=%d) != nb (=%d); alpha=%g. Arg. out of range?/n"),			 x, ncalc, nb, alpha);      else	MATHLIB_WARNING2(_("bessel_j(%g,nu=%g): precision lost in result/n"),			 x, alpha+(double)nb-1);    }    x = bj[nb-1];#ifdef MATHLIB_STANDALONE    free(bj);#else    vmaxset(vmax);#endif    return x;}

SEXP RTcl_AssignObjToVar(SEXP args){    const void *vmax = vmaxget();    Tcl_SetVar2Ex(RTcl_interp,		  translateChar(STRING_ELT(CADR(args), 0)),		  NULL,		  (Tcl_Obj *) R_ExternalPtrAddr(CADDR(args)),		  0);    vmaxset(vmax);    return R_NilValue;}

/* Match what EncodeString does with encodings */attribute_hiddenint Rstrlen(SEXP s, int quote){    cetype_t ienc = getCharCE(s);    if (ienc == CE_UTF8 || ienc == CE_BYTES)	return Rstrwid(CHAR(s), LENGTH(s), ienc, quote);    const void *vmax = vmaxget();    const char *p = translateChar(s);    int len = Rstrwid(p, (int)strlen(p), CE_NATIVE, quote);    vmaxset(vmax);    return len;}

static SEXP cross_colon(SEXP call, SEXP s, SEXP t){    SEXP a, la, ls, lt, rs, rt;    int i, j, k, n, nls, nlt;    char *cbuf;    const void *vmax = vmaxget();    if (length(s) != length(t))	errorcall(call, _("unequal factor lengths"));    n = length(s);    ls = getLevelsAttrib(s);    lt = getLevelsAttrib(t);    nls = LENGTH(ls);    nlt = LENGTH(lt);    PROTECT(a = allocVector(INTSXP, n));    PROTECT(rs = coerceVector(s, INTSXP));    PROTECT(rt = coerceVector(t, INTSXP));    for (i = 0; i < n; i++) {	int vs = INTEGER(rs)[i];	int vt = INTEGER(rt)[i];	if ((vs == NA_INTEGER) || (vt == NA_INTEGER))	    INTEGER(a)[i] = NA_INTEGER;	else	    INTEGER(a)[i] = vt + (vs - 1) * nlt;    }    UNPROTECT(2);    if (!isNull(ls) && !isNull(lt)) {	PROTECT(la = allocVector(STRSXP, nls * nlt));	k = 0;	/* FIXME: possibly UTF-8 version */	for (i = 0; i < nls; i++) {	    const char *vi = translateChar(STRING_ELT(ls, i));	    size_t vs = strlen(vi);	    for (j = 0; j < nlt; j++) {		const char *vj = translateChar(STRING_ELT(lt, j));		size_t vt = strlen(vj), len = vs + vt + 2;		cbuf = R_AllocStringBuffer(len, &cbuff);		snprintf(cbuf, len, "%s:%s", vi, vj);		SET_STRING_ELT(la, k, mkChar(cbuf));		k++;	    }	}	setAttrib(a, R_LevelsSymbol, la);	UNPROTECT(1);    }    PROTECT(la = mkString("factor"));    setAttrib(a, R_ClassSymbol, la);    UNPROTECT(2);    R_FreeStringBufferL(&cbuff);    vmaxset(vmax);    return a;}

double bessel_y(double x, double alpha){    long nb, ncalc;    double na, *by;#ifndef MATHLIB_STANDALONE    const void *vmax;#endif#ifdef IEEE_754    /* NaNs propagated correctly */    if (ISNAN(x) || ISNAN(alpha)) return x + alpha;#endif    if (x < 0) {	ML_ERROR(ME_RANGE, "bessel_y");	return ML_NAN;    }    na = floor(alpha);    if (alpha < 0) {	/* Using Abramowitz & Stegun  9.1.2	 * this may not be quite optimal (CPU and accuracy wise) */	return(bessel_y(x, -alpha) * cos(M_PI * alpha) -	       ((alpha == na) ? 0 :		bessel_j(x, -alpha) * sin(M_PI * alpha)));    }    nb = 1+ (long)na;/* nb-1 <= alpha < nb */    alpha -= (nb-1);#ifdef MATHLIB_STANDALONE    by = (double *) calloc(nb, sizeof(double));    if (!by) MATHLIB_ERROR("%s", _("bessel_y allocation error"));#else    vmax = vmaxget();    by = (double *) R_alloc((size_t) nb, sizeof(double));#endif    Y_bessel(&x, &alpha, &nb, by, &ncalc);    if(ncalc != nb) {/* error input */	if(ncalc == -1)	    return ML_POSINF;	else if(ncalc < -1)	    MATHLIB_WARNING4(_("bessel_y(%g): ncalc (=%ld) != nb (=%ld); alpha=%g. Arg. out of range?/n"),			     x, ncalc, nb, alpha);	else /* ncalc >= 0 */	    MATHLIB_WARNING2(_("bessel_y(%g,nu=%g): precision lost in result/n"),			     x, alpha+nb-1);    }    x = by[nb-1];#ifdef MATHLIB_STANDALONE    free(by);#else    vmaxset(vmax);#endif    return x;}

void Sn0(double *x, Sint *n, Sint *is_sorted, double *res, double *a2){    char *vmax;    vmax = vmaxget();    *res = sn0(x, (int)*n, (int)*is_sorted, a2);#ifdef DEBUG_Sno    REprintf("Sn0(* -> res=%g)/n", *res);#endif    vmaxset(vmax);}

/* utils::loadRconsole */SEXP in_loadRconsole(SEXP sfile){    struct structGUI gui;    const void *vmax = vmaxget();    if (!isString(sfile) || LENGTH(sfile) < 1)	error(_("invalid '%s' argument"), "file");    getActive(&gui);  /* Will get defaults if there's no active console */    if (loadRconsole(&gui, translateChar(STRING_ELT(sfile, 0)))) applyGUI(&gui);    if (strlen(gui.warning)) warning(gui.warning);    vmaxset(vmax);    return R_NilValue;}

SEXP RTcl_ObjFromVar(SEXP args){    Tcl_Obj *tclobj;    const void *vmax = vmaxget();    tclobj = Tcl_GetVar2Ex(RTcl_interp,                           translateChar(STRING_ELT(CADR(args), 0)),                           NULL,                           0);    SEXP res = makeRTclObject(tclobj);    vmaxset(vmax);    return res;}

SEXP FFI_Call(SEXP recv,SEXP sig,SEXP args) {  void (*FUN)() = NULL;  char *vmax = vmaxget();  void *handle = NULL;  /* Get the function pointer */  switch(TYPEOF(recv)) {  case EXTPTRSXP:    FUN = R_ExternalPtrAddr(recv);    break;  case STRSXP:    if(length(recv)==1) {      FUN = dlsym(NULL,CHAR(STRING_ELT(recv,0)));      if(NULL == FUN)	error("Pointer doesn't point to a function?");    } else {      handle = dlopen(CHAR(STRING_ELT(recv,0)),RTLD_NOW|RTLD_LOCAL);      if(NULL == handle)	error("Unable to load library %s",CHAR(STRING_ELT(recv,0)));      FUN    = dlsym(handle,CHAR(STRING_ELT(recv,1)));    }    if(NULL == FUN)      error("Unable to bind symbol %s",CHAR(STRING_ELT(recv,length(recv)-1)));    break;  default:    error("No receiver function address");    break;  }    /* Get the signature for the call */  ffi_cif *cif = R_ExternalPtrAddr(sig);  if(cif->nargs != length(args))    error("this signature specifies %d arguments",cif->nargs);    /* Sort out the arguments */  int    i;  void **values = (void**)R_alloc(sizeof(void),length(args));  for(i=0;i<length(args);i++)    values[i] = ffi_make_arg(TYPE_ELT(sig,i),cif->arg_types[i]->type,VECTOR_ELT(args,i));  /* Make the call and return */  ret_value retval;  SEXP      rval;  ffi_call(cif,FUN,&retval,values);  PROTECT(rval = ffi_fromType(TYPE_RET(sig),cif->rtype->type,&retval));  setAttrib(rval,install("ffi.args"),args);  vmaxset(vmax);  if(handle != NULL) dlclose(handle);  UNPROTECT(1);  return rval;}

attribute_hiddenvoid R_GE_VText(double x, double y, const char * const s, cetype_t enc,		double x_justify, double y_justify, double rotation,		const pGEcontext gc, pGEDevDesc dd){    if(!initialized) vfonts_Init();    if(initialized > 0) {	const void *vmax = vmaxget();	const char *str = reEnc(s, enc, CE_LATIN1, 2 /* '.' */);	(*routines.GEVText)(x, y, str, x_justify, y_justify, rotation, gc, dd);	vmaxset(vmax);    } else	error(_("Hershey fonts cannot be loaded"));}

SEXP dotTcl(SEXP args){    SEXP ans;    const char *cmd;    Tcl_Obj *val;    const void *vmax = vmaxget();    if(!isValidString(CADR(args)))	error(_("invalid argument"));    cmd = translateChar(STRING_ELT(CADR(args), 0));    val = tk_eval(cmd);    ans = makeRTclObject(val);    vmaxset(vmax);    return ans;}

attribute_hiddendouble R_GE_VStrWidth(const char *s, cetype_t enc, const pGEcontext gc, pGEDevDesc dd){    double res;    if(!initialized) vfonts_Init();    if(initialized > 0) {	const void *vmax = vmaxget();	const char *str = reEnc(s, enc, CE_LATIN1, 2 /* '.' */);	res = (*routines.GEVStrWidth)(str, gc, dd);	vmaxset(vmax);	return res;    } else {	error(_("Hershey fonts cannot be loaded"));	return 0.0; /* -Wall */    }}

SEXP RTcl_RemoveArrayElem(SEXP args){    SEXP x, i;    const char *xstr, *istr;    const void *vmax = vmaxget();    x = CADR(args);    i = CADDR(args);    xstr = translateChar(STRING_ELT(x, 0));    istr = translateChar(STRING_ELT(i, 0));    Tcl_UnsetVar2(RTcl_interp, xstr, istr, 0);    vmaxset(vmax);    return R_NilValue;}

static void qr_solve(Array x, Array y, Array coef)/* Translation of the R function qr.solve into pure C   NB We have to transpose the matrices since the ordering of an array is different in Fortran   NB2 We have to copy x to avoid it being overwritten.*/{    int i, info = 0, rank, *pivot, n, p;    const void *vmax;    double tol = 1.0E-7, *qraux, *work;    Array xt, yt, coeft;    assert(NROW(x) == NROW(y));    assert(NCOL(coef) == NCOL(y));    assert(NCOL(x) == NROW(coef));    vmax = vmaxget();    qraux = (double *) R_alloc(NCOL(x), sizeof(double));    pivot = (int *) R_alloc(NCOL(x), sizeof(int));    work  = (double *) R_alloc(2*NCOL(x), sizeof(double));    for(i = 0; i < NCOL(x); i++)	pivot[i] = i+1;    xt = make_zero_matrix(NCOL(x), NROW(x));    transpose_matrix(x,xt);    n = NROW(x);    p = NCOL(x);    F77_CALL(dqrdc2)(VECTOR(xt), &n, &n, &p, &tol, &rank,		       qraux, pivot, work);    if (rank != p)	error(_("Singular matrix in qr_solve"));    yt = make_zero_matrix(NCOL(y), NROW(y));    coeft = make_zero_matrix(NCOL(coef), NROW(coef));    transpose_matrix(y, yt);    F77_CALL(dqrcf)(VECTOR(xt), &NROW(x), &rank, qraux,	yt.vec, &NCOL(y), coeft.vec, &info);    transpose_matrix(coeft,coef);    vmaxset(vmax);}

static void whittle(Array acf, int nlag, Array *A, Array *B, Array p_forward,    Array v_forward, Array p_back, Array v_back){    int lag, nser = DIM(acf)[1];    const void *vmax;    Array EA, EB;	/* prediction variance */    Array KA, KB;	/* partial correlation coefficient */    Array id, tmp;    vmax = vmaxget();    KA = make_zero_matrix(nser, nser);    EA = make_zero_matrix(nser, nser);    KB = make_zero_matrix(nser, nser);    EB = make_zero_matrix(nser, nser);    id = make_identity_matrix(nser);    copy_array(id, subarray(A[0],0));    copy_array(id, subarray(B[0],0));    copy_array(id, subarray(p_forward,0));    copy_array(id, subarray(p_back,0));    for (lag = 1; lag <= nlag; lag++) {	whittle2(acf, A[lag-1], B[lag-1], lag, "forward", A[lag], KA, EB);	whittle2(acf, B[lag-1], A[lag-1], lag, "back", B[lag], KB, EA);	copy_array(EA, subarray(v_forward,lag-1));	copy_array(EB, subarray(v_back,lag-1));	copy_array(KA, subarray(p_forward,lag));	copy_array(KB, subarray(p_back,lag));    }    tmp = make_zero_matrix(nser,nser);    matrix_prod(KB,KA, 1, 1, tmp);    array_op(id, tmp, '-', tmp);    matrix_prod(EA, tmp, 0, 0, subarray(v_forward, nlag));    vmaxset(vmax);}

double bessel_k(double x, double alpha, double expo){    long nb, ncalc, ize;    double *bk;#ifndef MATHLIB_STANDALONE    const void *vmax;#endif#ifdef IEEE_754    /* NaNs propagated correctly */    if (ISNAN(x) || ISNAN(alpha)) return x + alpha;#endif    if (x < 0) {	ML_ERROR(ME_RANGE, "bessel_k");	return ML_NAN;    }    ize = (long)expo;    if(alpha < 0)	alpha = -alpha;    nb = 1+ (long)floor(alpha);/* nb-1 <= |alpha| < nb */    alpha -= (double)(nb-1);#ifdef MATHLIB_STANDALONE    bk = (double *) calloc(nb, sizeof(double));    if (!bk) MATHLIB_ERROR("%s", _("bessel_k allocation error"));#else    vmax = vmaxget();    bk = (double *) R_alloc((size_t) nb, sizeof(double));#endif    K_bessel(&x, &alpha, &nb, &ize, bk, &ncalc);    if(ncalc != nb) {/* error input */      if(ncalc < 0)	MATHLIB_WARNING4(_("bessel_k(%g): ncalc (=%ld) != nb (=%ld); alpha=%g. Arg. out of range?/n"),			 x, ncalc, nb, alpha);      else	MATHLIB_WARNING2(_("bessel_k(%g,nu=%g): precision lost in result/n"),			 x, alpha+(double)nb-1);    }    x = bk[nb-1];#ifdef MATHLIB_STANDALONE    free(bk);#else    vmaxset(vmax);#endif    return x;}

SEXP _gregexpr(SEXP _pattern, SEXP _text, SEXP _subpattern){    SEXP ans;    pcre *re_pcre;    pcre_extra *re_pe;    int erroffset;    const char *errorptr;    //int options = 0;    const unsigned char* tables = pcre_maketables();        const char* spat = as_string(_pattern);    re_pcre = pcre_compile(spat, 0, &errorptr, &erroffset, tables);    if (!re_pcre)    {        if (errorptr)        Rprintf("PCRE pattern compilation error/n/t'%s'/n/tat '%s'/n",            errorptr, spat+erroffset);        Rprintf("invalid regular expression '%s'/n", spat);    }    re_pe = pcre_study(re_pcre, 0, &errorptr);        int n = LENGTH(_text);    SEXP elt;    PROTECT(ans = allocVector(VECSXP, n));    const void *vmax = vmaxget();        for (int i = 0 ; i < n ; i++) {        if (STRING_ELT(_text, i) == NA_STRING) {        elt = gregexpr_NAInputAns();        } else {            const char* s = as_string(_text,i);            elt = _pcre(spat, s, re_pcre, re_pe, *REAL(_subpattern));        }        SET_VECTOR_ELT(ans, i, elt);        vmaxset(vmax);    }        if (re_pe) pcre_free(re_pe);    pcre_free(re_pcre);    pcre_free((void *)tables);        UNPROTECT(1);    return ans;    }

SEXP RTcl_SetArrayElem(SEXP args){    SEXP x, i;    const char *xstr, *istr;    Tcl_Obj *value;    const void *vmax = vmaxget();    x = CADR(args);    i = CADDR(args);    value = (Tcl_Obj *) R_ExternalPtrAddr(CADDDR(args));    xstr = translateChar(STRING_ELT(x, 0));    istr = translateChar(STRING_ELT(i, 0));    Tcl_SetVar2Ex(RTcl_interp, xstr, istr, value, 0);    vmaxset(vmax);    return R_NilValue;}

/* based on EncodeEnvironment in  printutils.c */static void PrintEnvironment(SEXP x){    const void *vmax = vmaxget();    if (x == R_GlobalEnv)	Rprintf("<R_GlobalEnv>");    else if (x == R_BaseEnv)	Rprintf("<base>");    else if (x == R_EmptyEnv)	Rprintf("<R_EmptyEnv>");    else if (R_IsPackageEnv(x))	Rprintf("<%s>",		translateChar(STRING_ELT(R_PackageEnvName(x), 0)));    else if (R_IsNamespaceEnv(x))	Rprintf("<namespace:%s>",		translateChar(STRING_ELT(R_NamespaceEnvSpec(x), 0)));    else Rprintf("<%p>", (void *)x);    vmaxset(vmax);}