/**
 * Get head of repository
 *
 * @param repo S4 class git_repository
 * @return R_NilValue if unborn branch or not found. S4 class
 * git_branch if not a detached head. S4 class git_commit if detached
 * head
 */
SEXP git2r_repository_head(SEXP repo)
{
    int err;
    SEXP result = R_NilValue;
    git_commit *commit = NULL;
    git_reference *reference = NULL;
    git_repository *repository = NULL;

    repository= git2r_repository_open(repo);
    if (!repository)
        git2r_error(git2r_err_invalid_repository, __func__, NULL);

    err = git_repository_head(&reference, repository);
    if (GIT_OK != err) {
        if (GIT_EUNBORNBRANCH == err || GIT_ENOTFOUND == err)
            err = GIT_OK;
        goto cleanup;
    }

    if (git_reference_is_branch(reference)) {
        git_branch_t type = GIT_BRANCH_LOCAL;
        if (git_reference_is_remote(reference))
            type = GIT_BRANCH_REMOTE;
        PROTECT(result = NEW_OBJECT(MAKE_CLASS("git_branch")));
        err = git2r_branch_init(reference, type, repo, result);
    } else {
        err = git_commit_lookup(
            &commit,
            repository,
            git_reference_target(reference));
        if (GIT_OK != err)
            goto cleanup;
        PROTECT(result = NEW_OBJECT(MAKE_CLASS("git_commit")));
        git2r_commit_init(commit, repo, result);
    }

cleanup:
    if (commit)
        git_commit_free(commit);

    if (reference)
        git_reference_free(reference);

    if (repository)
        git_repository_free(repository);

    if (R_NilValue != result)
        UNPROTECT(1);

    if (GIT_OK != err)
        git2r_error(git2r_err_from_libgit2, __func__, giterr_last()->message);

    return result;
}

SEXP ddense_band(SEXP x, SEXP k1P, SEXP k2P)
/* Always returns a full matrix with entries outside the band zeroed
 * Class of the value can be dtrMatrix or dgeMatrix
 */
{
    SEXP aa, ans = PROTECT(dup_mMatrix_as_dgeMatrix(x));
    int *adims = INTEGER(GET_SLOT(ans, Matrix_DimSym)),
	i, j, k1 = asInteger(k1P), k2 = asInteger(k2P);
    int m = adims[0], n = adims[1], sqr = (adims[0] == adims[1]),
	tru = (k1 >= 0), trl = (k2 <= 0);
    double *ax = REAL(GET_SLOT(ans, Matrix_xSym));

    if (k1 > k2)
	error(_("Lower band %d > upper band %d"), k1, k2);
    for (j = 0; j < n; j++) {
	int i1 = j - k2, i2 = j + 1 - k1;
	for (i = 0; i < i1; i++) ax[i + j * m] = 0.;
	for (i = i2; i < m; i++) ax[i + j * m] = 0.;
    }
    if (!sqr || (!tru && !trl)) { /* return the dgeMatrix */
	UNPROTECT(1);
	return ans;
    }
    /* Copy ans to a dtrMatrix object (must be square) */
    /* Because slots of ans are freshly allocated and ans will not be
     * used, we use the slots themselves and don't duplicate */
    aa = PROTECT(NEW_OBJECT(MAKE_CLASS("dtrMatrix")));
    SET_SLOT(aa, Matrix_xSym, GET_SLOT(ans, Matrix_xSym));
    SET_SLOT(aa, Matrix_DimSym, GET_SLOT(ans, Matrix_DimSym));
    SET_SLOT(aa, Matrix_DimNamesSym, GET_SLOT(ans, Matrix_DimNamesSym));
    SET_SLOT(aa, Matrix_diagSym, mkString("N"));
    SET_SLOT(aa, Matrix_uploSym, mkString(tru ? "U" : "L"));
    UNPROTECT(2);
    return aa;
}

SEXP dppMatrix_chol(SEXP x)
{
    SEXP val = get_factors(x, "pCholesky"),
	dimP = GET_SLOT(x, Matrix_DimSym),
	uploP = GET_SLOT(x, Matrix_uploSym);
    const char *uplo = CHAR(STRING_ELT(uploP, 0));
    int *dims = INTEGER(dimP), info;

    if (val != R_NilValue) return val;
    dims = INTEGER(dimP);
    val = PROTECT(NEW_OBJECT(MAKE_CLASS("pCholesky")));
    SET_SLOT(val, Matrix_uploSym, duplicate(uploP));
    SET_SLOT(val, Matrix_diagSym, mkString("N"));
    SET_SLOT(val, Matrix_DimSym, duplicate(dimP));
    slot_dup(val, x, Matrix_xSym);
    F77_CALL(dpptrf)(uplo, dims, REAL(GET_SLOT(val, Matrix_xSym)), &info);
    if (info) {
	if(info > 0) /* e.g. x singular */
	    error(_("the leading minor of order %d is not positive definite"),
		    info);
	else /* should never happen! */
	    error(_("Lapack routine %s returned error code %d"), "dpptrf", info);
    }
    UNPROTECT(1);
    return set_factors(x, val, "pCholesky");
}

SEXP rgeos_geospoint2SpatialPoints(SEXP env, GEOSGeom geom, SEXP p4s, SEXP id, int n) {
    
    GEOSContextHandle_t GEOShandle = getContextHandle(env);
        
    int type = GEOSGeomTypeId_r(GEOShandle, geom);    
    if ( type != GEOS_POINT && type != GEOS_MULTIPOINT && type != GEOS_GEOMETRYCOLLECTION )
        error("rgeos_geospoint2SpatialPoints: invalid geometry type");
    
    int pc=0;
    SEXP bbox, crdmat;
    if (GEOSisEmpty_r(GEOShandle, geom))
        error("rgeos_geospoint2SpatialPoints: empty point found");
    //if (GEOSisEmpty_r(GEOShandle, geom)==0) {
        PROTECT(bbox = rgeos_geom2bbox(env, geom)); pc++;
        PROTECT(crdmat = rgeos_geospoint2crdMat(env, geom, id, n, type)); pc++;
    //} else {
    //    bbox = R_NilValue;
    //    crdmat = R_NilValue;
    //}
    
    SEXP ans;
    PROTECT(ans = NEW_OBJECT(MAKE_CLASS("SpatialPoints"))); pc++;    
    SET_SLOT(ans, install("coords"), crdmat);
    SET_SLOT(ans, install("bbox"), bbox);
    SET_SLOT(ans, install("proj4string"), p4s);

    UNPROTECT(pc);
    return(ans);
}

/**
 * Callback when iterating over notes
 *
 * @param blob_id Oid of the blob containing the message
 * @param annotated_object_id Oid of the git object being annotated
 * @param payload Payload data passed to `git_note_foreach`
 * @return int 0 or error code
 */
static int git2r_note_foreach_cb(
    const git_oid *blob_id,
    const git_oid *annotated_object_id,
    void *payload)
{
    git2r_note_foreach_cb_data *cb_data = (git2r_note_foreach_cb_data*)payload;

    /* Check if we have a list to populate */
    if (R_NilValue != cb_data->list) {
        int err;
        SEXP note;

        SET_VECTOR_ELT(
            cb_data->list,
            cb_data->n,
            note = NEW_OBJECT(MAKE_CLASS("git_note")));

        err = git2r_note_init(
                  blob_id,
                  annotated_object_id,
                  cb_data->repository,
                  cb_data->notes_ref,
                  cb_data->repo,
                  note);
        if (GIT_OK != err)
            return err;
    }

    cb_data->n += 1;

    return 0;
}

/**
 * Copy the contents of N to a csn_LU or csn_QR object and,
 * optionally, free N or free both N and the pointers to its contents.
 *
 * @param a csn object to be converted
 * @param cl the name of the S4 class of the object to be generated
 * @param dofree 0 - don't free a; > 0 cs_free a; < 0 Free a
 *
 * @return SEXP containing a copy of S
 */
SEXP Matrix_csn_to_SEXP(csn *N, char *cl, int dofree)
{
    SEXP ans;
    char *valid[] = {"csn_LU", "csn_QR", ""};
    int ctype = Matrix_check_class(cl, valid), n = (N->U)->n;

    if (ctype < 0)
	error("Inappropriate class '%s' for Matrix_csn_to_SEXP", cl);
    ans = PROTECT(NEW_OBJECT(MAKE_CLASS(cl)));
				/* allocate and copy common slots */
    /* FIXME: Use the triangular matrix classes for csn_LU */
    SET_SLOT(ans, install("L"),	/* these are free'd later if requested */
	     Matrix_cs_to_SEXP(N->L, "dgCMatrix", 0));
    SET_SLOT(ans, install("U"),
	     Matrix_cs_to_SEXP(N->U, "dgCMatrix", 0));
    switch(ctype) {
    case 0:
	Memcpy(INTEGER(ALLOC_SLOT(ans, install("Pinv"), INTSXP, n)),
	       N->pinv, n);
	break;
    case 1:
	Memcpy(REAL(ALLOC_SLOT(ans, install("beta"), REALSXP, n)),
	       N->B, n);
	break;
    default:
	error("Inappropriate class '%s' for Matrix_csn_to_SEXP", cl);
    }
    if (dofree > 0) cs_nfree(N);
    if (dofree < 0) {
	Free(N->L); Free(N->U); Free(N);
    }
    UNPROTECT(1);
    return ans;
}

/**
 * Copy the contents of a to an appropriate CsparseMatrix object and,
 * optionally, free a or free both a and the pointers to its contents.
 *
 * @param a matrix to be converted
 * @param cl the name of the S4 class of the object to be generated
 * @param dofree 0 - don't free a; > 0 cs_free a; < 0 Free a
 *
 * @return SEXP containing a copy of a
 */
SEXP Matrix_cs_to_SEXP(cs *a, char *cl, int dofree)
{
    SEXP ans;
    char *valid[] = {"dgCMatrix", "dsCMatrix", "dtCMatrix", ""};
    int *dims, ctype = Matrix_check_class(cl, valid), nz;

    if (ctype < 0)
	error("invalid class of object to Matrix_cs_to_SEXP");
    ans = PROTECT(NEW_OBJECT(MAKE_CLASS(cl)));
				/* allocate and copy common slots */
    dims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));
    dims[0] = a->m; dims[1] = a->n;
    Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, a->n + 1)),
	   a->p, a->n + 1);
    nz = a->p[a->n];
    Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)), a->i, nz);
    Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), a->x, nz);
    if (ctype > 0) {
	int uplo = is_sym(a);
	if (!uplo) error("cs matrix not compatible with class");
	SET_SLOT(ans, Matrix_diagSym, mkString("N"));
	SET_SLOT(ans, Matrix_uploSym, mkString(uplo < 0 ? "L" : "U"));
    }
    if (dofree > 0) cs_spfree(a);
    if (dofree < 0) Free(a);
    UNPROTECT(1);
    return ans;
}

/**
 * Copy the contents of S to a css_LU or css_QR object and,
 * optionally, free S or free both S and the pointers to its contents.
 *
 * @param a css object to be converted
 * @param cl the name of the S4 class of the object to be generated
 * @param dofree 0 - don't free a; > 0 cs_free a; < 0 Free a
 * @param m number of rows in original matrix
 * @param n number of columns in original matrix
 *
 * @return SEXP containing a copy of S
 */
SEXP Matrix_css_to_SEXP(css *S, char *cl, int dofree, int m, int n)
{
    SEXP ans;
    char *valid[] = {"css_LU", "css_QR", ""};
    int *nz, ctype = Matrix_check_class(cl, valid);

    if (ctype < 0)
	error("Inappropriate class '%s' for Matrix_css_to_SEXP", cl);
    ans = PROTECT(NEW_OBJECT(MAKE_CLASS(cl)));
				/* allocate and copy common slots */
    Memcpy(INTEGER(ALLOC_SLOT(ans, install("Q"), INTSXP, n)), S->q, n);
    nz = INTEGER(ALLOC_SLOT(ans, install("nz"), INTSXP, 3));
    nz[0] = S->m2; nz[1] = S->lnz; nz[2] = S->unz;
    switch(ctype) {
    case 0:
	break;
    case 1:
	Memcpy(INTEGER(ALLOC_SLOT(ans, install("Pinv"), INTSXP, m)),
	       S->pinv, m);
	Memcpy(INTEGER(ALLOC_SLOT(ans, install("parent"), INTSXP, n)),
	       S->parent, n);
	Memcpy(INTEGER(ALLOC_SLOT(ans, install("cp"), INTSXP, n)),
	       S->cp, n);
	break;
    default:
	error("Inappropriate class '%s' for Matrix_css_to_SEXP", cl);
    }
    if (dofree > 0) cs_sfree(S);
    if (dofree < 0) Free(S);
    UNPROTECT(1);
    return ans;
}

/**
 * Fetch new data and update tips
 *
 * @param repo S4 class git_repository
 * @param name The name of the remote to fetch from
 * @param credentials The credentials for remote repository access.
 * @param msg The one line long message to be appended to the reflog
 * @return R_NilValue
 */
SEXP git2r_remote_fetch(
    SEXP repo,
    SEXP name,
    SEXP credentials,
    SEXP msg)
{
    int err;
    SEXP result = R_NilValue;
    const git_transfer_progress *stats;
    git_remote *remote = NULL;
    git_repository *repository = NULL;
    git_remote_callbacks callbacks = GIT_REMOTE_CALLBACKS_INIT;

    if (git2r_arg_check_string(name))
        git2r_error(git2r_err_string_arg, __func__, "name");
    if (git2r_arg_check_credentials(credentials))
        git2r_error(git2r_err_credentials_arg, __func__, "credentials");
    if (git2r_arg_check_string(msg))
        git2r_error(git2r_err_string_arg, __func__, "msg");

    repository = git2r_repository_open(repo);
    if (!repository)
        git2r_error(git2r_err_invalid_repository, __func__, NULL);

    err = git_remote_lookup(&remote, repository, CHAR(STRING_ELT(name, 0)));
    if (GIT_OK != err)
        goto cleanup;

    callbacks.credentials = &git2r_cred_acquire_cb;
    callbacks.payload = credentials;
    err = git_remote_set_callbacks(remote, &callbacks);
    if (GIT_OK != err)
        goto cleanup;

    err = git_remote_fetch(remote, NULL,  CHAR(STRING_ELT(msg, 0)));
    if (GIT_OK != err)
        goto cleanup;

    stats = git_remote_stats(remote);
    PROTECT(result = NEW_OBJECT(MAKE_CLASS("git_transfer_progress")));
    git2r_transfer_progress_init(stats, result);

cleanup:
    if (remote) {
        if (git_remote_connected(remote))
            git_remote_disconnect(remote);
        git_remote_free(remote);
    }

    if (repository)
        git_repository_free(repository);

    if (R_NilValue != result)
        UNPROTECT(1);

    if (GIT_OK != err)
        git2r_error(git2r_err_from_libgit2, __func__, giterr_last()->message);

    return result;
}

SEXP tsc_to_dgTMatrix(SEXP x)
{
    SEXP ans;
    if (*diag_P(x) != 'U')
	ans = compressed_to_dgTMatrix(x, ScalarLogical(1));
    else {			/* unit triangular matrix */
	SEXP islot = GET_SLOT(x, Matrix_iSym),
	    pslot = GET_SLOT(x, Matrix_pSym);
	int *ai, *aj, j,
	    n = length(pslot) - 1,
	    nod = length(islot),
	    nout = n + nod,
	    *p = INTEGER(pslot);
	double *ax;

	ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgTMatrix")));
	SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(x, Matrix_DimSym)));
	SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nout));
	ai = INTEGER(GET_SLOT(ans, Matrix_iSym));
	Memcpy(ai, INTEGER(islot), nod);
	SET_SLOT(ans, Matrix_jSym, allocVector(INTSXP, nout));
	aj = INTEGER(GET_SLOT(ans, Matrix_jSym));
	SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nout));
	ax = REAL(GET_SLOT(ans, Matrix_xSym));
	Memcpy(ax, REAL(GET_SLOT(x, Matrix_xSym)), nod);
	for (j = 0; j < n; j++) {
	    int jj, npj = nod + j,  p2 = p[j+1];
	    aj[npj] = ai[npj] = j;
	    ax[npj] = 1.;
	    for (jj = p[j]; jj < p2; jj++) aj[jj] = j;
	}
	UNPROTECT(1);
    }
    return ans;
}

static ad_conn_t *conn_new(ad_server_t *server, struct bufferevent *buffer) {
    if (server == NULL || buffer == NULL) {
        return NULL;
    }

    // Create a new connection container.
    ad_conn_t *conn = NEW_OBJECT(ad_conn_t);
    if (conn == NULL) return NULL;

    // Initialize with default values.
    conn->server = server;
    conn->buffer = buffer;
    conn->in = bufferevent_get_input(buffer);
    conn->out = bufferevent_get_output(buffer);
    conn_reset(conn);

    // Bind callback
    bufferevent_setcb(buffer, conn_read_cb, conn_write_cb, conn_event_cb, (void *)conn);
    bufferevent_setwatermark(buffer, EV_WRITE, 0, 0);
    bufferevent_enable(buffer, EV_WRITE);
    bufferevent_enable(buffer, EV_READ);

    // Run callbacks with AD_EVENT_INIT event.
    conn->status = call_hooks(AD_EVENT_INIT | AD_EVENT_WRITE, conn);

    return conn;
}

SEXP dsyMatrix_dgeMatrix_mm_R(SEXP a, SEXP b)
{
    int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
	*bdims = INTEGER(GET_SLOT(b, Matrix_DimSym)),
	*cdims,
	m = adims[0], n = bdims[1], k = adims[1];
    SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
    double one = 1., zero = 0.;

    if (bdims[0] != k)
	error(_("Matrices are not conformable for multiplication"));
    if (m < 1 || n < 1 || k < 1)
	error(_("Matrices with zero extents cannot be multiplied"));
    SET_SLOT(val, Matrix_rcondSym, allocVector(REALSXP, 0));
    SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));
    SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, m * n));
    SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));
    cdims = INTEGER(GET_SLOT(val, Matrix_DimSym));
    cdims[0] = m; cdims[1] = n;
    F77_CALL(dsymm)("R", uplo_P(a), adims, bdims+1, &one,
		    REAL(GET_SLOT(a, Matrix_xSym)), adims,
		    REAL(GET_SLOT(b, Matrix_xSym)), bdims,
		    &zero, REAL(GET_SLOT(val, Matrix_xSym)), adims);
    UNPROTECT(1);
    return val;
}

SEXP dpoMatrix_chol(SEXP x)
{
    SEXP val = get_factors(x, "Cholesky"),
         dimP = GET_SLOT(x, Matrix_DimSym),
         uploP = GET_SLOT(x, Matrix_uploSym);
    char *uplo = CHAR(STRING_ELT(uploP, 0));
    int *dims = INTEGER(dimP), info;
    int n = dims[0];
    double *vx;

    if (val != R_NilValue) return val;
    dims = INTEGER(dimP);
    val = PROTECT(NEW_OBJECT(MAKE_CLASS("Cholesky")));
    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);
    if (n > 0) {
        F77_CALL(dpotrf)(uplo, &n, vx, &n, &info);
        if (info) {
            if(info > 0)
                error(_("the leading minor of order %d is not positive definite"),
                      info);
            else /* should never happen! */
                error(_("Lapack routine %s returned error code %d"), "dpotrf", info);
        }
    }
    UNPROTECT(1);
    return set_factors(x, val, "Cholesky");
}

void strong_generators::init_from_ascii_coding(action *A, BYTE *ascii_coding, INT verbose_level)
{
	INT f_v = (verbose_level >= 1);
	INT f_vv = (verbose_level >= 2);
	longinteger_object go;
	group *G;

	if (f_v) {
		cout << "strong_generators::init_from_ascii_coding" << endl;
		}
	G = NEW_OBJECT(group);
	G->init(A);
	if (f_vv) {
		cout << "strong_generators::init_from_ascii_coding before G->init_ascii_coding_to_sims" << endl;
		}
	G->init_ascii_coding_to_sims(ascii_coding);
	if (f_vv) {
		cout << "strong_generators::init_from_ascii_coding after G->init_ascii_coding_to_sims" << endl;
		}
		

	G->S->group_order(go);

	if (f_vv) {
		cout << "strong_generators::init_from_ascii_coding Group order=" << go << endl;
		}

	init_from_sims(G->S, 0 /* verbose_level */);

	FREE_OBJECT(G);
	if (f_v) {
		cout << "strong_generators::init_from_ascii_coding done" << endl;
		}
}

SEXP csc_to_dgeMatrix(SEXP x)
{
    SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix"))),
	Dimslot = GET_SLOT(x, Matrix_DimSym);
    int *dims = INTEGER(Dimslot),
	*xp = INTEGER(GET_SLOT(x, Matrix_pSym)),
	*xi = INTEGER(GET_SLOT(x, Matrix_iSym));
    double *xx = REAL(GET_SLOT(x, Matrix_xSym)), *ax;
    int j, nrow = dims[0], ncol = dims[1];

    SET_SLOT(ans, Matrix_DimSym, duplicate(Dimslot));
    SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nrow*ncol));
    SET_SLOT(ans, Matrix_rcondSym, allocVector(REALSXP, 0));
    SET_SLOT(ans, Matrix_factorSym, allocVector(VECSXP, 0));
    ax = REAL(GET_SLOT(ans, Matrix_xSym));
    for (j = 0; j < (nrow * ncol); j++) ax[j] = 0.;
    for (j = 0; j < ncol; j++) {
	int ind;
	for (ind = xp[j]; ind < xp[j+1]; ind++) {
	    ax[j * nrow + xi[ind]] = xx[ind];
	}
    }
    UNPROTECT(1);
    return ans;
}

/**
 * Invoked 'callback' for each entry in the given FETCH_HEAD file.
 *
 * @param ref_name The name of the ref.
 * @param remote_url The url of the remote.
 * @param oid The id of the remote head that were updated during the
 * last fetch.
 * @param is_merge Is head for merge.
 * @return 0
 */
static int git2r_repository_fetchhead_foreach_cb(
    const char *ref_name,
    const char *remote_url,
    const git_oid *oid,
    unsigned int is_merge,
    void *payload)
{
    git2r_fetch_head_cb_data *cb_data = (git2r_fetch_head_cb_data*)payload;

    /* Check if we have a list to populate */
    if (R_NilValue != cb_data->list) {
        char sha[GIT_OID_HEXSZ + 1];
        SEXP fetch_head;

        SET_VECTOR_ELT(
            cb_data->list,
            cb_data->n,
            fetch_head = NEW_OBJECT(MAKE_CLASS("git_fetch_head")));

        SET_SLOT(fetch_head, Rf_install("ref_name"), mkString(ref_name));
        SET_SLOT(fetch_head, Rf_install("remote_url"), mkString(remote_url));
        git_oid_tostr(sha, sizeof(sha), oid);
        SET_SLOT(fetch_head, Rf_install("sha"), mkString(sha));
        SET_SLOT(fetch_head, Rf_install("is_merge"), ScalarLogical(is_merge));
        SET_SLOT(fetch_head, Rf_install("repo"), cb_data->repo);
    }

    cb_data->n += 1;

    return 0;
}

void flag_orbit_node::read_file(ifstream &fp, int verbose_level)
{
	int f_v = (verbose_level >= 1);
	
	if (f_v) {
		cout << "flag_orbit_node::read_file" << endl;
		}
	
	fp.read((char *) &downstep_primary_orbit, sizeof(int));
	fp.read((char *) &downstep_secondary_orbit, sizeof(int));
	fp.read((char *) &downstep_orbit_len, sizeof(int));
	fp.read((char *) &upstep_primary_orbit, sizeof(int));
	fp.read((char *) &upstep_secondary_orbit, sizeof(int));
	fp.read((char *) &f_fusion_node, sizeof(int));
	if (f_fusion_node) {
		fp.read((char *) &fusion_with, sizeof(int));
		fusion_elt = NEW_int(Flag_orbits->A->elt_size_in_int);
		Flag_orbits->A->element_read_from_file_binary(fusion_elt, fp, 0);
		}

	if (FALSE) {
		cout << "flag_orbit_node::read_file "
				"before gens->read_from_file_binary" << endl;
		}
	gens = NEW_OBJECT(strong_generators);
	gens->read_from_file_binary(Flag_orbits->A, fp, verbose_level);

	if (f_v) {
		cout << "flag_orbit_node::read_file finished" << endl;
		}
}

SEXP csc_matrix_crossprod(SEXP x, SEXP y, SEXP classed)
{
    int cl = asLogical(classed);
    SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
    int *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
	*ydims = INTEGER(cl ? GET_SLOT(y, Matrix_DimSym) :
			 getAttrib(y, R_DimSymbol)),
	*vdims = INTEGER(ALLOC_SLOT(val, Matrix_DimSym, INTSXP, 2));
    int *xi = INTEGER(GET_SLOT(x, Matrix_iSym)),
	*xp = INTEGER(GET_SLOT(x, Matrix_pSym));
    int j, k = xdims[0], m = xdims[1], n = ydims[1];
    double *vx, *xx = REAL(GET_SLOT(x, Matrix_xSym)),
	*yx = REAL(cl ? GET_SLOT(y, Matrix_xSym) : y);

    if (!cl && !(isMatrix(y) && isReal(y)))
	error(_("y must be a numeric matrix"));
    if (ydims[0] != k)
	error(_("x and y must have the same number of rows"));
    if (m < 1 || n < 1 || k < 1)
	error(_("Matrices with zero extents cannot be multiplied"));
    vdims[0] = m; vdims[1] = n;
    vx = REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, m * n));
    for (j = 0; j < n; j++) {
	int i; double *ypt = yx + j * k;
	for(i = 0; i < m; i++) {
	    int ii; double accum = 0.;
	    for (ii = xp[i]; ii < xp[i+1]; ii++) {
		accum += xx[ii] * ypt[xi[ii]];
	    }
	    vx[i + j * m] = accum;
	}
    }
    UNPROTECT(1);
    return val;
}

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

SEXP new_LinStatExpectCovarMPinv(int p, int q) {

    SEXP ans, expect, covar, linearstatistic, MPinv;

    PROTECT(ans = NEW_OBJECT(MAKE_CLASS("LinStatExpectCovarMPinv")));

    SET_SLOT(ans, PL2_expectationSym, expect = PROTECT(allocVector(REALSXP, p * q)));
    for (int i = 0; i < p * q; i++)
        REAL(expect)[i] = 0.0;

    SET_SLOT(ans, PL2_covarianceSym, covar = PROTECT(allocMatrix(REALSXP, p * q, p * q)));
    for (int i = 0; i < p * q * p * q; i++)
        REAL(covar)[i] = 0.0;

    SET_SLOT(ans, PL2_dimensionSym, PROTECT(ScalarInteger(p * q)));

    SET_SLOT(ans, PL2_linearstatisticSym, linearstatistic = PROTECT(allocVector(REALSXP, p * q)));
    for (int i = 0; i < p * q; i++)
        REAL(linearstatistic)[i] = 0.0;

    SET_SLOT(ans, PL2_MPinvSym, MPinv = PROTECT(allocMatrix(REALSXP, p * q, p * q)));
    for (int i = 0; i < p * q * p * q; i++)
        REAL(MPinv)[i] = 0.0;

    SET_SLOT(ans, PL2_rankSym, PROTECT(ScalarReal(0.0)));

    SET_SLOT(ans, PL2_svdmemSym, PROTECT(new_svd_mem(p * q)));

    SET_SLOT(ans, PL2_expcovinfSym, PROTECT(new_ExpectCovarInfluence(q)));

    UNPROTECT(9);
    return(ans);
}