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

}

static void
loess_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) ); )