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

SEXP lsq_dense_QR(SEXP X, SEXP y){    SEXP ans;    int info, n, p, k, *Xdims, *ydims, lwork;    double *work, tmp, *xvals;    if (!(isReal(X) & isMatrix(X)))	error(_("X must be a numeric (double precision) matrix"));    Xdims = INTEGER(coerceVector(getAttrib(X, R_DimSymbol), INTSXP));    n = Xdims[0];    p = Xdims[1];    if (!(isReal(y) & isMatrix(y)))	error(_("y must be a numeric (double precision) matrix"));    ydims = INTEGER(coerceVector(getAttrib(y, R_DimSymbol), INTSXP));    if (ydims[0] != n)	error(_(	    "number of rows in y (%d) does not match number of rows in X (%d)"),	    ydims[0], n);    k = ydims[1];    if (k < 1 || p < 1) return allocMatrix(REALSXP, p, k);    xvals = (double *) R_alloc(n * p, sizeof(double));    Memcpy(xvals, REAL(X), n * p);    ans = PROTECT(duplicate(y));    lwork = -1;    F77_CALL(dgels)("N", &n, &p, &k, xvals, &n, REAL(ans), &n,		    &tmp, &lwork, &info);    if (info)	error(_("First call to Lapack routine dgels returned error code %d"),	      info);    lwork = (int) tmp;    work = (double *) R_alloc(lwork, sizeof(double));    F77_CALL(dgels)("N", &n, &p, &k, xvals, &n, REAL(ans), &n,		    work, &lwork, &info);    if (info)	error(_("Second call to Lapack routine dgels returned error code %d"),	      info);    UNPROTECT(1);    return ans;}

SEXP dtrMatrix_rcond(SEXP obj, SEXP type){    char typnm[] = {'/0', '/0'};    int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym)), info;    double rcond;    typnm[0] = rcond_type(CHAR(asChar(type)));    F77_CALL(dtrcon)(typnm, uplo_P(obj), diag_P(obj), dims,                     REAL(GET_SLOT(obj, Matrix_xSym)), dims, &rcond,                     (double *) R_alloc(3*dims[0], sizeof(double)),                     (int *) R_alloc(dims[0], sizeof(int)), &info);    return ScalarReal(rcond);}

/* function called by Fortran to check for output */static void C_solout_bim (int * m, int *k, int * ord,   double * t0, double * tstep, double * y, double * f,   double *dd, double * rpar, int * ipar, int * irtrn){  *irtrn = 1;  while ((*t0 <= tt[it]) && (tt[it] < tstep[*k-1])) { 	  F77_CALL(contsolall) (&tt[it], m, k, t0, tstep, dd, ytmp);    saveOut(tt[it], ytmp);	  it++;	  if (it >= maxt) break;  }}

//x %*% t(y)void tcrossprod(double *x, int* nrx, int* ncx,		      double *y, int* nry, int* ncy, double *z){    char *transa = "N", *transb = "T";    double one = 1.0, zero = 0.0;    if (*nrx > 0 && *ncx > 0 && *nry > 0 && *ncy > 0) {	F77_CALL(dgemm)(transa, transb, nrx, nry, ncx, &one,			x, nrx, y, nry, &zero, z, nrx);    } else { /* zero-extent operations should return zeroes */	int i;	for(i = 0; i < (*nrx)*(*nry); i++) z[i] = 0;    }}

static void tcrossprod(double *x, int nrx, int ncx,		      double *y, int nry, int ncy, double *z){    char *transa = "N", *transb = "T";    double one = 1.0, zero = 0.0;    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {	F77_CALL(dgemm)(transa, transb, &nrx, &nry, &ncx, &one,			x, &nrx, y, &nry, &zero, z, &nrx);    } else { /* zero-extent operations should return zeroes */	R_xlen_t NRX = nrx;	for(R_xlen_t i = 0; i < NRX*nry; i++) z[i] = 0;    }}

SEXP dppMatrix_rcond(SEXP obj, SEXP type){    SEXP Chol = dppMatrix_chol(obj);    char typnm[] = {'O', '/0'};	/* always use the one norm */    int *dims = INTEGER(GET_SLOT(Chol, Matrix_DimSym)), info;    double anorm = get_norm_sp(obj, typnm), rcond;    F77_CALL(dppcon)(uplo_P(Chol), dims,		     REAL(GET_SLOT(Chol, Matrix_xSym)), &anorm, &rcond,		     (double *) R_alloc(3*dims[0], sizeof(double)),		     (int *) R_alloc(dims[0], sizeof(int)), &info);    return ScalarReal(rcond);}

CHM_DN Cholesky_rd::solveA(CHM_DN rhs) {    int info, nrhs = (int)rhs->ncol;    if (n != (int)rhs->nrow)	error(_("%s dimension mismatch: lhs of size %d, rhs has %d rows"),	      "Cholesky_rd::solveA", n, rhs->nrow);    CHM_DN ans = M_cholmod_copy_dense(rhs, &c);    F77_CALL(dpotrs)(uplo, &n, &nrhs, X, &n,		     (double*)ans->x, &n, &info);    if (info)	error(_("dpotrs in Cholesky_rd::solveA returned error code %d"),	      info);    return ans;}

static voidHF_fact(double *par, longint *time, longint *n, double *mat, double *logdet){    longint job = 11L, info, i, nsq = *n * (*n), np1 = *n + 1;    double *work = Calloc(*n, double), *work1 = Calloc(nsq, double);#ifndef USING_R    longint zero = 0L;#endif    HF_mat(par, time, n, mat);#ifdef USING_R    F77_CALL(chol) (mat, n, n, mat, &info);#else    F77_CALL(chol) (mat, n, work, &zero, &zero, &info);#endif    for(i = 0; i < *n; i++) {	work1[i * np1] = 1;	F77_CALL(dtrsl) (mat, n, n, work1 + i * (*n), &job, &info);	*logdet -= log(fabs(mat[i * np1]));    }    Memcpy(mat, work1, nsq);    Free(work); Free(work1);}

/* C-level function to compute Moore-Penrose Generalized Inverse of a square matrix. */void c_ginv(double *covariance, int ncols, double *mpinv) {int i = 0, j = 0, errcode = 0;double *u = NULL, *d = NULL, *vt = NULL, *backup = NULL;double sv_tol = 0, zero = 0, one = 1;char transa = 'N', transb = 'N';  c_udvt(&u, &d, &vt, ncols);  if (covariance != mpinv) {    backup = Calloc1D(ncols * ncols, sizeof(double));    memcpy(backup, covariance, ncols * ncols * sizeof(double));  }/*THEN*/  /* compute the SVD decomposition. */  c_svd(covariance, u, d, vt, &ncols, &ncols, &ncols, FALSE, &errcode);  /* if SVD fails, catch the error code and free all buffers. */  if (errcode == 0) {    /* set the threshold for the singular values as in corpcor. */    sv_tol = ncols * d[0] * MACHINE_TOL * MACHINE_TOL;    /* the first multiplication, U * D^{-1} is easy. */    for (i = 0; i < ncols; i++)      for (j = 0; j < ncols; j++)        u[CMC(i, j, ncols)] = u[CMC(i, j, ncols)] * ((d[j] > sv_tol) ? 1/d[j] : 0);    /* the second one, (U * D^{-1}) * Vt  is a real matrix multiplication. */    F77_CALL(dgemm)(&transa, &transb, &ncols, &ncols, &ncols, &one, u,      &ncols, vt, &ncols, &zero, mpinv, &ncols);  }/*THEN*/  if (covariance != mpinv) {    memcpy(covariance, backup, ncols * ncols * sizeof(double));    Free1D(backup);  }/*THEN*/  Free1D(u);  Free1D(d);  Free1D(vt);  if (errcode)    error("an error (%d) occurred in the call to c_ginv()./n", errcode);}/*C_GINV*/

SEXP dgeMatrix_svd(SEXP x, SEXP nnu, SEXP nnv){    int /* nu = asInteger(nnu),	   nv = asInteger(nnv), */	*dims = INTEGER(GET_SLOT(x, Matrix_DimSym));    double *xx = REAL(GET_SLOT(x, Matrix_xSym));    SEXP val = PROTECT(allocVector(VECSXP, 3));    if (dims[0] && dims[1]) {	int m = dims[0], n = dims[1], mm = (m < n)?m:n,	    lwork = -1, info;	double tmp, *work;	int *iwork, n_iw = 8 * mm;	C_or_Alloca_TO(iwork, n_iw, int);	SET_VECTOR_ELT(val, 0, allocVector(REALSXP, mm));	SET_VECTOR_ELT(val, 1, allocMatrix(REALSXP, m, mm));	SET_VECTOR_ELT(val, 2, allocMatrix(REALSXP, mm, n));	F77_CALL(dgesdd)("S", &m, &n, xx, &m,			 REAL(VECTOR_ELT(val, 0)),			 REAL(VECTOR_ELT(val, 1)), &m,			 REAL(VECTOR_ELT(val, 2)), &mm,			 &tmp, &lwork, iwork, &info);	lwork = (int) tmp;	C_or_Alloca_TO(work, lwork, double);	F77_CALL(dgesdd)("S", &m, &n, xx, &m,			 REAL(VECTOR_ELT(val, 0)),			 REAL(VECTOR_ELT(val, 1)), &m,			 REAL(VECTOR_ELT(val, 2)), &mm,			 work, &lwork, iwork, &info);	if(n_iw  >= SMALL_4_Alloca) Free(iwork);	if(lwork >= SMALL_4_Alloca) Free(work);    }    UNPROTECT(1);    return val;}

/* **** CRSF_chol2inv ****  * This function is a C interface to the fortran implemented  * scalapack driver function "callpdpotri" that performs * inverting a matrix from its Choleski Factorization */ int CRSF_chol2inv(int dim[], int iMyRank) {	int iMemSize = 0;	double *dpWork = NULL;	int ipZero[] = { 0, 1, 2, 3 };	int NPRow = dim[6];	int NPCol = dim[7];	int MyRow = iMyRank / NPCol;	int MyCol = iMyRank % NPCol;	int rowOfA = dim[0];	int colOfA = dim[1];	int rowBlockSize = dim[4];	int colBlockSize = dim[5];	/* Calculate required memory size */	int localRowSizeOfA = F77_CALL(numroc)(&rowOfA, &rowBlockSize, &MyRow, ipZero, &NPRow);	int localColSizeOfA = F77_CALL(numroc)(&colOfA, &colBlockSize, &MyCol, ipZero, &NPCol);		int localSizeOfA = localRowSizeOfA * localColSizeOfA;	int workSpace = max (rowBlockSize, colBlockSize);	iMemSize = localSizeOfA + workSpace;		dpWork = (double *) malloc(sizeof(double) * iMemSize);	memset(dpWork, 0xcc, sizeof(double) * iMemSize);	D_Rprintf (("After allocating memory .. /n "));		F77_CALL(callpdpotri)(dim, dpWork, &iMemSize);	D_Rprintf (("AFTER FORTRAN FUNCTION EXECUTION /n "));	free (dpWork);	return 0;}

staticdouble get_norm(SEXP obj, const char *typstr){    char typnm[] = {'/0', '/0'};    int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym));    double *work = (double *) NULL;    typnm[0] = norm_type(typstr);    if (*typnm == 'I') {        work = (double *) R_alloc(dims[0], sizeof(double));    }    return F77_CALL(dlantr)(typnm, uplo_P(obj), diag_P(obj), dims, dims+1,                            REAL(GET_SLOT(obj, Matrix_xSym)), dims, work);}

double calculateLambdaMax(int *n, int *p, double *X, double *U, double *y,                           double *D, int *degrees, int *cum_degrees, int *numcolsU,                           int *family, double gamma) {  double curr_max = 0.0;  double norm = 0.0;  double trDinv;  for(int j=0;j<*p;j++){    trDinv = 0.0;    double *Ujy = malloc(degrees[j]*sizeof(double));    // Calculate alpha norm    norm = fabs(F77_CALL(ddot)(n, X+(*n)*j, &inc_one, y, &inc_one))/gamma;    curr_max = max(curr_max, norm);    // Calculate beta norm    F77_CALL(dgemv)("T",n,degrees+j,&one,U+(*n)*(cum_degrees[j]),n,y,      &inc_one, &zero, Ujy, &inc_one);    for(int i=0; i<degrees[j];i++) {      trDinv += 1/D[cum_degrees[j] + i];    }    // Calculate norm of D^{-1/2}Ujy and scale    free(Ujy);  }  return curr_max;}

SEXP dppMatrix_solve(SEXP x){    SEXP Chol = dppMatrix_chol(x);    SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dppMatrix")));    int *dims = INTEGER(GET_SLOT(x, Matrix_DimSym)), info;    slot_dup(val, Chol, Matrix_uploSym);    slot_dup(val, Chol, Matrix_xSym);    slot_dup(val, Chol, Matrix_DimSym);    F77_CALL(dpptri)(uplo_P(val), dims,		     REAL(GET_SLOT(val, Matrix_xSym)), &info);    UNPROTECT(1);    return val;}

static void symtcrossprod(double *x, int nr, int nc, double *z){    char *trans = "N", *uplo = "U";    double one = 1.0, zero = 0.0;    if (nr > 0 && nc > 0) {	F77_CALL(dsyrk)(uplo, trans, &nr, &nc, &one, x, &nr, &zero, z, &nr);	for (int i = 1; i < nr; i++)	    for (int j = 0; j < i; j++) z[i + nr *j] = z[j + nr * i];    } else { /* zero-extent operations should return zeroes */	R_xlen_t NR = nr;	for(R_xlen_t i = 0; i < NR*NR; i++) z[i] = 0;    }}

static voidloess_grow(int *parameter, int *a, double *xi,	   double *vert, double *vval){    int d, vc, nc, nv, a1, v1, xi1, vv1, i, k;    d = parameter[0];    vc = parameter[2];    nc = parameter[3];    nv = parameter[4];    liv = parameter[5];    lv = parameter[6];    iv = Calloc(liv, int);    v = Calloc(lv, double);    iv[1] = d;    iv[2] = parameter[1];    iv[3] = vc;    iv[5] = iv[13] = nv;    iv[4] = iv[16] = nc;    iv[6] = 50;    iv[7] = iv[6] + nc;    iv[8] = iv[7] + vc * nc;    iv[9] = iv[8] + nc;    iv[10] = 50;    iv[12] = iv[10] + nv * d;    iv[11] = iv[12] + (d + 1) * nv;    iv[27] = 173;    v1 = iv[10] - 1;    xi1 = iv[11] - 1;    a1 = iv[6] - 1;    vv1 = iv[12] - 1;    for(i = 0; i < d; i++) {	k = nv * i;	v[v1 + k] = vert[i];	v[v1 + vc - 1 + k] = vert[i + d];    }    for(i = 0; i < nc; i++) {	v[xi1 + i] = xi[i];	iv[a1 + i] = a[i];    }    k = (d + 1) * nv;    for(i = 0; i < k; i++)	v[vv1 + i] = vval[i];    F77_CALL(ehg169)(&d, &vc, &nc, &nc, &nv, &nv, v+v1, iv+a1,		    v+xi1, iv+iv[7]-1, iv+iv[8]-1, iv+iv[9]-1);}

void matrix_inverse(Matrix *X, Matrix *X_inverse, Matrix *Xsamedims){  int n=numrows(X), e_code, ipiv[n];  // Need to set X_inverse to the identity matrix on input:  matrix_identity(X_inverse);  // Copy X to Xsamedims (error check for dims inside matrix_copy):  matrix_copy(X, Xsamedims);  // Compute: Solution to a real system of linear equations: A * X = B  // Where A is an N-by-N matrix and X and B are N-by-NRHS matrices.  // The LU decomposition with partial pivoting and row interchanges is  // used to factor A as A = P * L * U,  // where P is a permutation matrix, L is unit lower triangular, and U is  // upper triangular.  The factored form of A is then used to solve the  // system of equations A * X = B.  //  // N    = The number of linear equations, i.e., numrows(A)  // NRHS = The number of right hand sides, i.e., numcols(B)  //  // A    = LDA-by-N matrix, the leading N-by-N matrix of A is the   //        coefficient matrix A. On exit, the factors L and U from the  //        factorization. A = P*L*U  // LDA = The leading dimension of the array A (LDA >= max(1,N))  //  // IPIV = N-vector containing the pivot indices that define P;  //        row i of the matrix was interchanged with row IPIV(i)  //  // B    = LDB-by-NRHS matrix, the leading N-by-NRHS matrix of B is the  //        right hand side matrix. On exit, the N-by-NRHS solution X.  //  // LDB = The leading dimension of the array B (LDB >= max(1,N))  // INFO  =0 => Successful exit  //       <0 => If INFO = -i, the i-th argument had an illegal value  //       >0 => If INFO = i, U(i,i) is exactly zero.  The factorization  //             has been completed, but the factor U is exactly  //              singular, so the solution could not be computed.//dgesv(n,n,Xsamedims,n,ipiv,X_inverse,n,&e_code);               // C version  F77_CALL(dgesv)(&n,&n,Xsamedims,&n,ipiv,X_inverse,&n,&e_code); // R version  if (!e_code)    return;  if (e_code<0)    error("Singular value in mat_inverse./n");  else     error("Illegal value in mat_inverse./n");  return;}

/** * Update the fixed effects and the orthogonal random effects in an MCMC sample * from an mer object. * * @param x an mer object * @param sigma current standard deviation of the per-observation *        noise terms. * @param fvals pointer to memory in which to store the updated beta * @param rvals pointer to memory in which to store the updated b (may *              be (double*)NULL) */static void MCMC_beta_u(SEXP x, double sigma, double *fvals, double *rvals){    int *dims = DIMS_SLOT(x);    int i1 = 1, p = dims[p_POS], q = dims[q_POS];    double *V = V_SLOT(x), *fixef = FIXEF_SLOT(x), *muEta = MUETA_SLOT(x),            *u = U_SLOT(x), mone[] = {-1,0}, one[] = {1,0};    CHM_FR L = L_SLOT(x);    double *del1 = Calloc(q, double), *del2 = Alloca(p, double);    CHM_DN sol, rhs = N_AS_CHM_DN(del1, q, 1);    R_CheckStack();    if (V || muEta) {        error(_("Update not yet written"));    } else {			/* Linear mixed model */        update_L(x);        update_RX(x);        lmm_update_fixef_u(x);        /* Update beta */        for (int j = 0; j < p; j++) del2[j] = sigma * norm_rand();        F77_CALL(dtrsv)("U", "N", "N", &p, RX_SLOT(x), &p, del2, &i1);        for (int j = 0; j < p; j++) fixef[j] += del2[j];        /* Update u */        for (int j = 0; j < q; j++) del1[j] = sigma * norm_rand();        F77_CALL(dgemv)("N", &q, &p, mone, RZX_SLOT(x), &q,                        del2, &i1, one, del1, &i1);        sol = M_cholmod_solve(CHOLMOD_Lt, L, rhs, &c);        for (int j = 0; j < q; j++) u[j] += ((double*)(sol->x))[j];        M_cholmod_free_dense(&sol, &c);        update_mu(x);	     /* and parts of the deviance slot */    }    Memcpy(fvals, fixef, p);    if (rvals) {        update_ranef(x);        Memcpy(rvals, RANEF_SLOT(x), q);    }    Free(del1);}

void C_dgesvd(int* jobu,int* jobv,int* nrx,int* ncx,double* x,double* s,double* u,double* vt,int* info){	char const jobs[] = "NOSA";	char JOBU[2];JOBU[0] = jobs[*jobu];JOBU[1] = '/0';	char JOBV[2];JOBV[0] = jobs[*jobv];JOBV[1] = '/0';	// Rprintf("jobi(%i %i) jobs(%s,%s)/n",*jobu,*jobv,&JOBU[0],&JOBV[0]);		// set leading dimensions to default values no matrices are submatrices here	int ldx = MAX(1,*nrx); 	int ldu = 1;	if((JOBU[0] == 'S') || (JOBU[0] == 'A'))		ldu = *nrx;	int ldvt = 1;	if(JOBV[0] == 'S')		ldvt = MIN(*nrx,*ncx); 	else if(JOBV[0] == 'A')				ldvt = *ncx;    // Rprintf("n=%i p=%i ldx=%i ldu=%i ldvt=%i/n",*nrx,*ncx,ldx,ldu,ldvt);	// dgesvd    int lwork = -1;	double _work;    F77_CALL(dgesvd)(JOBU, JOBV, nrx,ncx,x,&ldx,s,u,&ldu, vt,&ldvt,&_work,&lwork,info);	if(*info){		Rprintf("Illegal arguments to Lapack routine '%s' returning error code %d", "dgesvd" ,*info);		return;	}	lwork = (int)_work;    double *work = (double *) malloc(lwork * sizeof(double));    F77_CALL(dgesvd)(JOBU, JOBV, nrx,ncx,x,&ldx,s,u,&ldu, vt,&ldvt,work,&lwork,info);	free(work);	if(*info){		Rprintf("error code %d from Lapack routine '%s'", *info, "dgesvd");		//return;	}		}

void hdrOutL( char *param,              char *xname,              char *item,              char *commen,              int *value,              int *status ) {  DECLARE_CHARACTER_DYN(fparam);  DECLARE_CHARACTER_DYN(fxname);  DECLARE_CHARACTER_DYN(fitem);  DECLARE_CHARACTER_DYN(fcommen);  F77_LOGICAL_TYPE *fvalue;  int i;  int nparam;  /*  Count the number of parameters and create a Fortran logical      array of the correct size */  nparam = img1CountParams( param, status );  fvalue = (F77_LOGICAL_TYPE *) malloc( nparam * sizeof(F77_LOGICAL_TYPE) );  /*  Convert the input values into Fortran logical values */  for ( i = 0; i < nparam; i++ ) {    if ( value[i] ) {      fvalue[i] = F77_TRUE;    } else {      fvalue[i] = F77_FALSE;    }  }  F77_CREATE_CHARACTER(fparam,strlen( param ));  cnf_exprt( param, fparam, fparam_length );  F77_CREATE_CHARACTER(fxname,strlen( xname ));  cnf_exprt( xname, fxname, fxname_length );  F77_CREATE_CHARACTER(fitem,strlen( item ));  cnf_exprt( item, fitem, fitem_length );  F77_CREATE_CHARACTER(fcommen,strlen( commen ));  cnf_exprt( commen, fcommen, fcommen_length );  F77_LOCK( F77_CALL(hdr_outl)( CHARACTER_ARG(fparam),                      CHARACTER_ARG(fxname),                      CHARACTER_ARG(fitem),                      CHARACTER_ARG(fcommen),                      LOGICAL_ARRAY_ARG(fvalue),                      INTEGER_ARG(status)                      TRAIL_ARG(fparam)                      TRAIL_ARG(fxname)                      TRAIL_ARG(fitem)                      TRAIL_ARG(fcommen) ); )

SEXP MulspeC(SEXP n, SEXP d, SEXP lag1, SEXP lag3, SEXP cov){    double *d1,*d2,*d3,*d4,*d5,*d6;    int *i1,*i2,*i3,*i4;    SEXP ans =  R_NilValue,  spec1 = R_NilValue, spec2 = R_NilValue, stat = R_NilValue,  coh1 = R_NilValue, coh2 = R_NilValue;    double *xspec1, *xspec2, *xstat, *xcoh1, *xcoh2 = NULL;    int   i, nd, nd2, lg1;    i1 = INTEGER_POINTER(n);    i2 = INTEGER_POINTER(d);    i3 = INTEGER_POINTER(lag1);    i4 = INTEGER_POINTER(lag3);    d1 = NUMERIC_POINTER(cov);    nd = *i2;    nd2 = nd * nd;    lg1 = *i3;    PROTECT(ans = allocVector(VECSXP, 5));    SET_VECTOR_ELT(ans, 0, spec1 = allocVector(REALSXP, lg1*nd2));    SET_VECTOR_ELT(ans, 1, spec2 = allocVector(REALSXP, lg1*nd2));    SET_VECTOR_ELT(ans, 2, stat = allocVector(REALSXP, lg1*nd));     SET_VECTOR_ELT(ans, 3, coh1 = allocVector(REALSXP, lg1*nd2));    SET_VECTOR_ELT(ans, 4, coh2 = allocVector(REALSXP, lg1*nd2));     d2 = NUMERIC_POINTER(spec1);    d3 = NUMERIC_POINTER(spec2);    d4 = NUMERIC_POINTER(stat);    d5 = NUMERIC_POINTER(coh1);    d6 = NUMERIC_POINTER(coh2);    F77_CALL(mulspef) (i1,i2,i3,i4,d1,d2,d3,d4,d5,d6);    xspec1 = REAL(spec1);    xspec2 = REAL(spec2);    xstat = REAL(stat);    xcoh1 = REAL(coh1);    xcoh2 = REAL(coh2);    for(i=0; i<lg1*nd2; i++) xspec1[i] = d2[i];    for(i=0; i<lg1*nd2; i++) xspec2[i] = d3[i];    for(i=0; i<lg1*nd; i++) xstat[i] = d4[i];    for(i=0; i<lg1*nd2; i++) xcoh1[i] = d5[i];    for(i=0; i<lg1*nd2; i++) xcoh2[i] = d6[i];    UNPROTECT(1);    return ans;}

/** * Generate the direction for "running Shake-and-Bake" according to 1.3.3 of * Boender et al. (1991) */void hitandrun_rsabDir(double *d, Matrix *constr, int index) {	const int inc1 = 1; // for BLAS	int n = constr->nCol - 1;	double c[n]; // the constraint vector	for (int i = 0; i < n; ++i) {		c[i] = *get(constr, index, i);	}  if (n == 1) {    d[0] = -c[0];    return;  }	double r = root(unif_rand(), n - 1);	hitandrun_randDir(d, n); // /~{u} in the paper	double cd = F77_CALL(ddot)(&n, c, &inc1, d, &inc1);	double fd = r / sqrt(1 - cd * cd);	double fc = -(r * cd / sqrt(1 - cd * cd) + sqrt(1 - r * r));	F77_CALL(dscal)(&n, &fd, d, &inc1); // d := fd * d	F77_CALL(daxpy)(&n, &fc, c, &inc1, d, &inc1); // d := fc * c + d}

voidloess_dfit(double *y, double *x, double *x_evaluate, double *weights,	   double *span, int *degree, int *nonparametric,	   int *drop_square, int *sum_drop_sqr,	   int *d, int *n, int *m, double *fit){    int zero = 0;    double dzero = 0.0;    loess_workspace(d, n, span, degree, nonparametric, drop_square,		    sum_drop_sqr, &zero);    F77_CALL(lowesf)(x, y, weights, iv, &liv, &lv, v, m, x_evaluate,		    &dzero, &zero, fit);    loess_free();}

voidloess_dfitse(double *y, double *x, double *x_evaluate, double *weights,	     double *robust, int *family, double *span, int *degree,	     int *nonparametric, int *drop_square,	     int *sum_drop_sqr,	     int *d, int *n, int *m, double *fit, double *L){    int zero = 0, two = 2;    double dzero = 0.0;    loess_workspace(d, n, span, degree, nonparametric, drop_square,		    sum_drop_sqr, &zero);    if(*family == GAUSSIAN)	F77_CALL(lowesf)(x, y, weights, iv, &liv, &lv, v, m,			x_evaluate, L, &two, fit);    else if(*family == SYMMETRIC)    {	F77_CALL(lowesf)(x, y, weights, iv, &liv, &lv, v, m,			x_evaluate, L, &two, fit);	F77_CALL(lowesf)(x, y, robust, iv, &liv, &lv, v, m,			x_evaluate, &dzero, &zero, fit);    }    loess_free();}

SEXP magma_dgeMatrix_svd(SEXP x, SEXP nnu, SEXP nnv) {#ifdef HIPLAR_WITH_MAGMA	int /* nu = asInteger(nnu),		   nv = asInteger(nnv), */		*dims = INTEGER(GET_SLOT(x, Matrix_DimSym));	double *xx = REAL(GET_SLOT(x, Matrix_xSym));	SEXP val = PROTECT(allocVector(VECSXP, 3));	if (dims[0] && dims[1]) {		int m = dims[0], n = dims[1], mm = (m < n)?m:n,			lwork = -1, info;		double tmp, *work;		int *iwork = Alloca(8 * mm, int);		R_CheckStack();		SET_VECTOR_ELT(val, 0, allocVector(REALSXP, mm));		SET_VECTOR_ELT(val, 1, allocMatrix(REALSXP, m, mm));		SET_VECTOR_ELT(val, 2, allocMatrix(REALSXP, mm, n));				if(GPUFlag == 0) {			F77_CALL(dgesdd)("S", &m, &n, xx, &m,					REAL(VECTOR_ELT(val, 0)),					REAL(VECTOR_ELT(val, 1)), &m,					REAL(VECTOR_ELT(val, 2)), &mm,					&tmp, &lwork, iwork, &info);			lwork = (int) tmp;			work = Alloca(lwork, double);			R_CheckStack();			F77_CALL(dgesdd)("S", &m, &n, xx, &m,					REAL(VECTOR_ELT(val, 0)),					REAL(VECTOR_ELT(val, 1)), &m,					REAL(VECTOR_ELT(val, 2)), &mm,					work, &lwork, iwork, &info);		} else {

static void tccrossprod(Rcomplex *x, int nrx, int ncx,			Rcomplex *y, int nry, int ncy, Rcomplex *z){    char *transa = "N", *transb = "T";    Rcomplex one, zero;    one.r = 1.0; one.i = zero.r = zero.i = 0.0;    if (nrx > 0 && ncx > 0 && nry > 0 && ncy > 0) {	F77_CALL(zgemm)(transa, transb, &nrx, &nry, &ncx, &one,			x, &nrx, y, &nry, &zero, z, &nrx);    } else { /* zero-extent operations should return zeroes */	R_xlen_t NRX = nrx;	for(R_xlen_t i = 0; i < NRX*nry; i++) z[i].r = z[i].i = 0;    }}

SEXP dsyMatrix_rcond(SEXP obj, SEXP type){    SEXP trf = dsyMatrix_trf(obj);    int *dims = INTEGER(GET_SLOT(obj, Matrix_DimSym)), info;    double anorm = get_norm_sy(obj, "O");    double rcond;    F77_CALL(dsycon)(uplo_P(trf), dims,		     REAL   (GET_SLOT(trf, Matrix_xSym)), dims,		     INTEGER(GET_SLOT(trf, Matrix_permSym)),		     &anorm, &rcond,		     (double *) R_alloc(2*dims[0], sizeof(double)),		     (int *) R_alloc(dims[0], sizeof(int)), &info);    return ScalarReal(rcond);}