C++ R_CheckStack函数代码示例

SEXP sparseQR_coef(SEXP qr, SEXP y)
{
    SEXP ans = PROTECT(dup_mMatrix_as_dgeMatrix(y)),
	qslot = GET_SLOT(qr, install("q"));
    CSP V = AS_CSP(GET_SLOT(qr, install("V"))),
	R = AS_CSP(GET_SLOT(qr, install("R")));
    int *ydims = INTEGER(GET_SLOT(ans, Matrix_DimSym)),
	*q = INTEGER(qslot),
	j, lq = LENGTH(qslot), m = R->m, n = R->n;
    double *ax = REAL(GET_SLOT(ans, Matrix_xSym)),
	*x = Alloca(m, double);
    R_CheckStack();
    R_CheckStack();

    /* apply row permutation and multiply by Q' */
    sparseQR_Qmult(V, REAL(GET_SLOT(qr, install("beta"))),
		   INTEGER(GET_SLOT(qr, Matrix_pSym)), 1,
		   REAL(GET_SLOT(ans, Matrix_xSym)), ydims);
    for (j = 0; j < ydims[1]; j++) {
	double *aj = ax + j * m;
	cs_usolve(R, aj);
	if (lq) {
	    cs_ipvec(q, aj, x, n);
	    Memcpy(aj, x, n);
	}
    }
    UNPROTECT(1);
    return ans;
}

SEXP dgCMatrix_matrix_solve(SEXP Ap, SEXP b, SEXP give_sparse)
{
    Rboolean sparse = asLogical(give_sparse);
    if(sparse) {
	// FIXME: implement this
	error(_("dgCMatrix_matrix_solve(.., sparse=TRUE) not yet implemented"));

	/* Idea: in the  for(j = 0; j < nrhs ..) loop below, build the *sparse* result matrix
	 * ----- *column* wise -- which is perfect for dgCMatrix
	 * --> build (i,p,x) slots "increasingly" [well, allocate in batches ..]
	 *
	 * --> maybe first a protoype in R
	 */

    }
    SEXP ans = PROTECT(dup_mMatrix_as_dgeMatrix(b)),
	lu, qslot;
    CSP L, U;
    int *bdims = INTEGER(GET_SLOT(ans, Matrix_DimSym)), *p, *q;
    int j, n = bdims[0], nrhs = bdims[1];
    double *ax = REAL(GET_SLOT(ans, Matrix_xSym)),
	*x = Alloca(n, double);
    R_CheckStack();

    if (isNull(lu = get_factors(Ap, "LU"))) {
	install_lu(Ap, /* order = */ 1, /* tol = */ 1.0, /* err_sing = */ TRUE);
	lu = get_factors(Ap, "LU");
    }
    qslot = GET_SLOT(lu, install("q"));
    L = AS_CSP__(GET_SLOT(lu, install("L")));
    U = AS_CSP__(GET_SLOT(lu, install("U")));
    R_CheckStack();

    p = INTEGER(GET_SLOT(lu, Matrix_pSym));
    q = LENGTH(qslot) ? INTEGER(qslot) : (int *) NULL;

    if (U->n != n || nrhs < 1 || n < 1)
	error(_("Dimensions of system to be solved are inconsistent"));
    for (j = 0; j < nrhs; j++) {
	cs_pvec(p, ax + j * n, x, n);  /* x = b(p) */
	cs_lsolve(L, x);	       /* x = L\x */
	cs_usolve(U, x);	       /* x = U\x */
	if (q)			       /* r(q) = x , hence r = Q' U{^-1} L{^-1} P b = A^{-1} b */
	    cs_ipvec(q, x, ax + j * n, n);
	else
	    Memcpy(ax + j * n, x, n);
    }
    UNPROTECT(1);
    return ans;
}

SEXP sparseQR_resid_fitted(SEXP qr, SEXP y, SEXP resid)
{
    SEXP ans = PROTECT(dup_mMatrix_as_dgeMatrix(y));
    CSP V = AS_CSP(GET_SLOT(qr, install("V")));
    int *ydims = INTEGER(GET_SLOT(ans, Matrix_DimSym)),
	*p = INTEGER(GET_SLOT(qr, Matrix_pSym)),
	i, j, m = V->m, n = V->n, res = asLogical(resid);
    double *ax = REAL(GET_SLOT(ans, Matrix_xSym)),
	*beta = REAL(GET_SLOT(qr, install("beta")));
    R_CheckStack();

    /* apply row permutation and multiply by Q' */
    sparseQR_Qmult(V, beta, p, 1, ax, ydims);
    for (j = 0; j < ydims[1]; j++) {
	if (res)		/* zero first n rows */
	    for (i = 0; i < n; i++) ax[i + j * m] = 0;
	else 			/* zero last m - n rows */
	    for (i = n; i < m; i++) ax[i + j * m] = 0;
    }
    /* multiply by Q and apply inverse row permutation */
    sparseQR_Qmult(V, beta, p, 0, ax, ydims);

    UNPROTECT(1);
    return ans;
}

/**
 * Apply Householder transformations and the row permutation P to y
 *
 * @param a sparse matrix containing the vectors defining the
 *        Householder transformations
 * @param beta scaling factors for the Householder transformations
 * @param y contents of a V->m by nrhs dense matrix
 * @param p 0-based permutation vector of length V->m
 * @param nrhs number of right hand sides (i.e. ncol(y))
 * @param trans logical value - if TRUE create Q'y[p] otherwise Qy[p]
 */
static
void sparseQR_Qmult(cs *V, double *beta, int *p, int trans,
		    double *y, int *ydims)
{
    int j, k, m = V->m, n = V->n;
    double *x = Alloca(m, double);	/* workspace */
    R_CheckStack();

    if (ydims[0] != m)
	error(_("Dimensions of system are inconsistent"));
    for (j = 0; j < ydims[1]; j++) {
	double *yj = y + j * m;
	if (trans) {
	    cs_pvec(p, yj, x, m);	/* x(0:m-1) = y(p(0:m-1, j)) */
	    Memcpy(yj, x, m);	/* replace it */
	    for (k = 0 ; k < n ; k++)   /* apply H[1]...H[n] */
		cs_happly(V, k, beta[k], yj);
	} else {
	    for (k = n - 1 ; k >= 0 ; k--) /* apply H[n]...H[1] */
		cs_happly(V, k, beta[k], yj);
	    cs_ipvec(p, yj, x, m); /* inverse permutation */
	    Memcpy(yj, x, m);
	}
    }
}

SEXP dsyMatrix_trf(SEXP x)
{
    SEXP val = get_factors(x, "BunchKaufman"),
	dimP = GET_SLOT(x, Matrix_DimSym),
	uploP = GET_SLOT(x, Matrix_uploSym);
    int *dims = INTEGER(dimP), *perm, info;
    int lwork = -1, n = dims[0];
    const char *uplo = CHAR(STRING_ELT(uploP, 0));
    double tmp, *vx, *work;

    if (val != R_NilValue) return val;
    dims = INTEGER(dimP);
    val = PROTECT(NEW_OBJECT(MAKE_CLASS("BunchKaufman")));
    SET_SLOT(val, Matrix_uploSym, duplicate(uploP));
    SET_SLOT(val, Matrix_diagSym, mkString("N"));
    SET_SLOT(val, Matrix_DimSym, duplicate(dimP));
    vx = REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, n * n));
    AZERO(vx, n * n);
    F77_CALL(dlacpy)(uplo, &n, &n, REAL(GET_SLOT(x, Matrix_xSym)), &n, vx, &n);
    perm = INTEGER(ALLOC_SLOT(val, Matrix_permSym, INTSXP, n));
    F77_CALL(dsytrf)(uplo, &n, vx, &n, perm, &tmp, &lwork, &info);
    lwork = (int) tmp;
    work = Alloca(lwork, double);
    R_CheckStack();
    F77_CALL(dsytrf)(uplo, &n, vx, &n, perm, work, &lwork, &info);
    if (info) error(_("Lapack routine dsytrf returned error code %d"), info);
    UNPROTECT(1);
    return set_factors(x, val, "BunchKaufman");
}

/**
 * "Indexing" aka subsetting : Compute  x[i,j], also for vectors i and j
 * Working via CHOLMOD_submatrix, see ./CHOLMOD/MatrixOps/cholmod_submatrix.c
 * @param x CsparseMatrix
 * @param i row     indices (0-origin), or NULL (R's)
 * @param j columns indices (0-origin), or NULL
 *
 * @return x[i,j]  still CsparseMatrix --- currently, this loses dimnames
 */
SEXP Csparse_submatrix(SEXP x, SEXP i, SEXP j)
{
    CHM_SP chx = AS_CHM_SP(x); /* << does diagU2N() when needed */
    int rsize = (isNull(i)) ? -1 : LENGTH(i),
	csize = (isNull(j)) ? -1 : LENGTH(j);
    int Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
    R_CheckStack();

    if (rsize >= 0 && !isInteger(i))
	error(_("Index i must be NULL or integer"));
    if (csize >= 0 && !isInteger(j))
	error(_("Index j must be NULL or integer"));

    if (!chx->stype) {/* non-symmetric Matrix */
	return chm_sparse_to_SEXP(cholmod_submatrix(chx,
						    (rsize < 0) ? NULL : INTEGER(i), rsize,
						    (csize < 0) ? NULL : INTEGER(j), csize,
						    TRUE, TRUE, &c),
				  1, 0, Rkind, "",
				  /* FIXME: drops dimnames */ R_NilValue);
    }
				/* for now, cholmod_submatrix() only accepts "generalMatrix" */
    CHM_SP tmp = cholmod_copy(chx, /* stype: */ 0, chx->xtype, &c);
    CHM_SP ans = cholmod_submatrix(tmp,
				   (rsize < 0) ? NULL : INTEGER(i), rsize,
				   (csize < 0) ? NULL : INTEGER(j), csize,
				   TRUE, TRUE, &c);
    cholmod_free_sparse(&tmp, &c);
    return chm_sparse_to_SEXP(ans, 1, 0, Rkind, "", R_NilValue);
}

SEXP Csparse_diagN2U(SEXP x)
{
    const char *cl = class_P(x);
    /* dtCMatrix, etc; [1] = the second character =?= 't' for triangular */
    if (cl[1] != 't' || *diag_P(x) != 'N') {
	/* "trivially fast" when not triangular (<==> no 'diag' slot),
	   or already *unit* triangular */
	return (x);
    }
    else { /* triangular with diag='N'): now drop the diagonal */
	/* duplicate, since chx will be modified: */
	SEXP xx = PROTECT(duplicate(x));
	CHM_SP chx = AS_CHM_SP__(xx);
	int uploT = (*uplo_P(x) == 'U') ? 1 : -1,
	    Rkind = (chx->xtype != CHOLMOD_PATTERN) ? Real_kind(x) : 0;
	R_CheckStack();

	chm_diagN2U(chx, uploT, /* do_realloc */ FALSE);

	UNPROTECT(1);
	return chm_sparse_to_SEXP(chx, /*dofree*/ 0/* or 1 ?? */,
				  uploT, Rkind, "U",
				  GET_SLOT(x, Matrix_DimNamesSym));
    }
}

// called from package MatrixModels's R code:
SEXP dgCMatrix_cholsol(SEXP x, SEXP y)
{
    /* Solve Sparse Least Squares X %*% beta ~= y  with dense RHS y,
     * where X = t(x) i.e. we pass  x = t(X)  as argument,
     * via  "Cholesky(X'X)" .. well not really:
     * cholmod_factorize("x", ..) finds L in  X'X = L'L directly */
    CHM_SP cx = AS_CHM_SP(x);
    /* FIXME: extend this to work in multivariate case, i.e. y a matrix with > 1 column ! */
    CHM_DN cy = AS_CHM_DN(coerceVector(y, REALSXP)), rhs, cAns, resid;
    CHM_FR L;
    int n = cx->ncol;/* #{obs.} {x = t(X) !} */
    double one[] = {1,0}, zero[] = {0,0}, neg1[] = {-1,0};
    const char *nms[] = {"L", "coef", "Xty", "resid", ""};
    SEXP ans = PROTECT(Rf_mkNamed(VECSXP, nms));
    R_CheckStack();

    if (n < cx->nrow || n <= 0)
	error(_("dgCMatrix_cholsol requires a 'short, wide' rectangular matrix"));
    if (cy->nrow != n)
	error(_("Dimensions of system to be solved are inconsistent"));
    rhs = cholmod_allocate_dense(cx->nrow, 1, cx->nrow, CHOLMOD_REAL, &c);
    /* cholmod_sdmult(A, transp, alpha, beta, X, Y, &c):
     *		Y := alpha*(A*X) + beta*Y or alpha*(A'*X) + beta*Y ;
     * here: rhs := 1 * x %*% y + 0 =  x %*% y =  X'y  */
    if (!(cholmod_sdmult(cx, 0 /* trans */, one, zero, cy, rhs, &c)))
	error(_("cholmod_sdmult error (rhs)"));
    L = cholmod_analyze(cx, &c);
    if (!cholmod_factorize(cx, L, &c))
	error(_("cholmod_factorize failed: status %d, minor %d from ncol %d"),
	      c.status, L->minor, L->n);
/* FIXME: Do this in stages so an "effects" vector can be calculated */
    if (!(cAns = cholmod_solve(CHOLMOD_A, L, rhs, &c)))
	error(_("cholmod_solve (CHOLMOD_A) failed: status %d, minor %d from ncol %d"),
	      c.status, L->minor, L->n);
    /* L : */
    SET_VECTOR_ELT(ans, 0, chm_factor_to_SEXP(L, 0));
    /* coef : */
    SET_VECTOR_ELT(ans, 1, allocVector(REALSXP, cx->nrow));
    Memcpy(REAL(VECTOR_ELT(ans, 1)), (double*)(cAns->x), cx->nrow);
    /* X'y : */
/* FIXME: Change this when the "effects" vector is available */
    SET_VECTOR_ELT(ans, 2, allocVector(REALSXP, cx->nrow));
    Memcpy(REAL(VECTOR_ELT(ans, 2)), (double*)(rhs->x), cx->nrow);
    /* resid := y */
    resid = cholmod_copy_dense(cy, &c);
    /* cholmod_sdmult(A, transp, alp, bet, X, Y, &c):
     *		Y := alp*(A*X) + bet*Y or alp*(A'*X) + beta*Y ;
     * here: resid := -1 * x' %*% coef + 1 * y = y - X %*% coef  */
    if (!(cholmod_sdmult(cx, 1/* trans */, neg1, one, cAns, resid, &c)))
	error(_("cholmod_sdmult error (resid)"));
    /* FIXME: for multivariate case, i.e. resid  *matrix* with > 1 column ! */
    SET_VECTOR_ELT(ans, 3, allocVector(REALSXP, n));
    Memcpy(REAL(VECTOR_ELT(ans, 3)), (double*)(resid->x), n);

    cholmod_free_factor(&L, &c);
    cholmod_free_dense(&rhs, &c);
    cholmod_free_dense(&cAns, &c);
    UNPROTECT(1);
    return ans;
}

//main function used .Call()
SEXP doCollisionTest(SEXP num, SEXP n, SEXP m)
{
        if (!isNumeric(num) || !isNumeric(n) || !isNumeric(m))
            error(_("invalid argument"));
        
        
        //temporary working variables
        int lenSample  = asInteger( n ); //length of the sample 'num'
        int nbUrns = asInteger( m ); //number of urns
        
        int * rNum = INTEGER( num ); // vector of length n with random urn numbers
        
        int * Urns = (int *) R_alloc(nbUrns, sizeof(int));
        
        //result
        int *nbCollision = (int *) R_alloc(1, sizeof(int));
        SEXP resultinR; //result in R
        PROTECT(resultinR = allocVector(INTSXP, 1)); //allocate an integer
        /*
         if(resultinR == NULL) Rprintf("zog\n");
         else Rprintf("toooo2\n");*/
        
        nbCollision = INTEGER( resultinR ); //plug the C pointer on the R type
        
        R_CheckStack();
        
        //computation step        
        collisionTest(rNum, lenSample, nbUrns, Urns, nbCollision);
        
        UNPROTECT(1);
        
        return resultinR;
}

SEXP dense_to_Csparse(SEXP x)
{
    CHM_DN chxd = AS_CHM_xDN(PROTECT(mMatrix_as_geMatrix(x)));
    /* cholmod_dense_to_sparse() in CHOLMOD/Core/ below does only work for
       "REAL" 'xtypes', i.e. *not* for "nMatrix".
       ===> need "_x" in above AS_CHM_xDN() call.

       Also it cannot keep symmetric / triangular, hence the
       as_geMatrix() above.  Note that this is already a *waste* for
       symmetric matrices; However, we could conceivably use an
       enhanced cholmod_dense_to_sparse(), with an extra boolean
       argument for symmetry.
    */
    CHM_SP chxs = cholmod_dense_to_sparse(chxd, 1, &c);
    int Rkind = (chxd->xtype == CHOLMOD_REAL) ? Real_KIND2(x) : 0;
    /* Note: when 'x' was integer Matrix, Real_KIND(x) = -1, but *_KIND2(.) = 0 */
    R_CheckStack();

    UNPROTECT(1);
    /* chm_sparse_to_SEXP() *could* deal with symmetric
     * if chxs had such an stype; and we should be able to use uplo below */
    return chm_sparse_to_SEXP(chxs, 1, 0/*TODO: uplo_P(x) if x has an uplo slot*/,
			      Rkind, "",
			      isMatrix(x) ? getAttrib(x, R_DimNamesSymbol)
			      : GET_SLOT(x, Matrix_DimNamesSym));
}

SEXP Csparse_dense_prod(SEXP a, SEXP b)
{
    CHM_SP cha = AS_CHM_SP(a);
    SEXP b_M = PROTECT(mMatrix_as_dgeMatrix(b));
    CHM_DN chb = AS_CHM_DN(b_M);
    CHM_DN chc = cholmod_l_allocate_dense(cha->nrow, chb->ncol, cha->nrow,
					chb->xtype, &c);
    SEXP dn = PROTECT(allocVector(VECSXP, 2));
    double one[] = {1,0}, zero[] = {0,0};
    int nprot = 2;
    R_CheckStack();
    /* Tim Davis, please FIXME:  currently (2010-11) *fails* when  a  is a pattern matrix:*/
    if(cha->xtype == CHOLMOD_PATTERN) {
	/* warning(_("Csparse_dense_prod(): cholmod_sdmult() not yet implemented for pattern./ ngCMatrix" */
	/* 	  " --> slightly inefficient coercion")); */

	// This *fails* to produce a CHOLMOD_REAL ..
	// CHM_SP chd = cholmod_l_copy(cha, cha->stype, CHOLMOD_REAL, &c);
	// --> use our Matrix-classes
	SEXP da = PROTECT(nz2Csparse(a, x_double)); nprot++;
	cha = AS_CHM_SP(da);
    }
    cholmod_l_sdmult(cha, 0, one, zero, chb, chc, &c);
    SET_VECTOR_ELT(dn, 0,	/* establish dimnames */
		   duplicate(VECTOR_ELT(GET_SLOT(a, Matrix_DimNamesSym), 0)));
    SET_VECTOR_ELT(dn, 1,
		   duplicate(VECTOR_ELT(GET_SLOT(b_M, Matrix_DimNamesSym), 1)));
    UNPROTECT(nprot);
    return chm_dense_to_SEXP(chc, 1, 0, dn);
}

static SEXP mkCharUcs(wchar_t *name)
{
    int n = wcslen(name), N = 3*n+1;
    char buf[N];
    R_CheckStack();
    wcstombs(buf, name, N); buf[N-1] = '\0';
    return mkCharCE(buf, CE_UTF8);
}

SEXP destructive_CHM_update(SEXP object, SEXP parent, SEXP mult)
{
    CHM_FR L = AS_CHM_FR(object);
    CHM_SP A = AS_CHM_SP__(parent);
    R_CheckStack();

    return chm_factor_to_SEXP(chm_factor_update(L, A, asReal(mult)), 0);
}

/* Modified version of Tim Davis's cs_lu_mex.c file for MATLAB */
void install_lu(SEXP Ap, int order, double tol, Rboolean err_sing)
{
    // (order, tol) == (1, 1) by default, when called from R.
    SEXP ans;
    css *S;
    csn *N;
    int n, *p, *dims;
    CSP A = AS_CSP__(Ap), D;
    R_CheckStack();

    n = A->n;
    if (A->m != n)
	error(_("LU decomposition applies only to square matrices"));
    if (order) {		/* not using natural order */
	order = (tol == 1) ? 2	/* amd(S'*S) w/dense rows or I */
	    : 1;		/* amd (A+A'), or natural */
    }
    S = cs_sqr(order, A, /*qr = */ 0);	/* symbolic ordering */
    N = cs_lu(A, S, tol);	/* numeric factorization */
    if (!N) {
	if(err_sing)
	    error(_("cs_lu(A) failed: near-singular A (or out of memory)"));
	else {
	    /* No warning: The useR should be careful :
	     * Put  NA  into  "LU" factor */
	    set_factors(Ap, ScalarLogical(NA_LOGICAL), "LU");
	    return;
	}
    }
    cs_dropzeros(N->L);		/* drop zeros from L and sort it */
    D = cs_transpose(N->L, 1);
    cs_spfree(N->L);
    N->L = cs_transpose(D, 1);
    cs_spfree(D);
    cs_dropzeros(N->U);		/* drop zeros from U and sort it */
    D = cs_transpose(N->U, 1);
    cs_spfree(N->U);
    N->U = cs_transpose(D, 1);
    cs_spfree(D);
    p = cs_pinv(N->pinv, n);	/* p=pinv' */
    ans = PROTECT(NEW_OBJECT(MAKE_CLASS("sparseLU")));
    dims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));
    dims[0] = n; dims[1] = n;
    SET_SLOT(ans, install("L"),
	     Matrix_cs_to_SEXP(N->L, "dtCMatrix", 0));
    SET_SLOT(ans, install("U"),
	     Matrix_cs_to_SEXP(N->U, "dtCMatrix", 0));
    Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_pSym, /* "p" */
			      INTSXP, n)), p, n);
    if (order)
	Memcpy(INTEGER(ALLOC_SLOT(ans, install("q"),
				  INTSXP, n)), S->q, n);
    cs_nfree(N);
    cs_sfree(S);
    cs_free(p);
    UNPROTECT(1);
    set_factors(Ap, ans, "LU");
}

SEXP Csparse_validate_(SEXP x, Rboolean maybe_modify)
{
    /* NB: we do *NOT* check a potential 'x' slot here, at all */
    SEXP pslot = GET_SLOT(x, Matrix_pSym),
	islot = GET_SLOT(x, Matrix_iSym);
    Rboolean sorted, strictly;
    int j, k,
	*dims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
	nrow = dims[0],
	ncol = dims[1],
	*xp = INTEGER(pslot),
	*xi = INTEGER(islot);

    if (length(pslot) != dims[1] + 1)
	return mkString(_("slot p must have length = ncol(.) + 1"));
    if (xp[0] != 0)
	return mkString(_("first element of slot p must be zero"));
    if (length(islot) < xp[ncol]) /* allow larger slots from over-allocation!*/
	return
	    mkString(_("last element of slot p must match length of slots i and x"));
    for (j = 0; j < xp[ncol]; j++) {
	if (xi[j] < 0 || xi[j] >= nrow)
	    return mkString(_("all row indices must be between 0 and nrow-1"));
    }
    sorted = TRUE; strictly = TRUE;
    for (j = 0; j < ncol; j++) {
	if (xp[j] > xp[j + 1])
	    return mkString(_("slot p must be non-decreasing"));
	if(sorted) /* only act if >= 2 entries in column j : */
	    for (k = xp[j] + 1; k < xp[j + 1]; k++) {
		if (xi[k] < xi[k - 1])
		    sorted = FALSE;
		else if (xi[k] == xi[k - 1])
		    strictly = FALSE;
	    }
    }
    if (!sorted) {
	if(maybe_modify) {
	    CHM_SP chx = (CHM_SP) alloca(sizeof(cholmod_sparse));
	    R_CheckStack();
	    as_cholmod_sparse(chx, x, FALSE, TRUE);/*-> cholmod_l_sort() ! */
	    /* as chx = AS_CHM_SP__(x)  but  ^^^^ sorting x in_place !!! */

	    /* Now re-check that row indices are *strictly* increasing
	     * (and not just increasing) within each column : */
	    for (j = 0; j < ncol; j++) {
		for (k = xp[j] + 1; k < xp[j + 1]; k++)
		    if (xi[k] == xi[k - 1])
			return mkString(_("slot i is not *strictly* increasing inside a column (even after cholmod_l_sort)"));
	    }
	} else { /* no modifying sorting : */
	    return mkString(_("row indices are not sorted within columns"));
	}
    } else if(!strictly) {  /* sorted, but not strictly */
	return mkString(_("slot i is not *strictly* increasing inside a column"));
    }
    return ScalarLogical(1);
}