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 void
HF_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;
}

static
double 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;
    }

}