/* copy a resource from the input file to the output file */
static void copyres(osfildef *fpin, osfildef *fpout, ulong siz,
                    uint endpos_ofs)
{
    ulong startpos;
    ulong endpos;
    uchar buf[4];

    /* note the starting position */
    startpos = osfpos(fpout);

    /* copy the bytes of the resource */
    copybytes(fpin, fpout, siz);

    /* note the ending position */
    endpos = osfpos(fpout);

    /* write the ending position at the appropriate point */
    osfseek(fpout, (ulong)(startpos + endpos_ofs), OSFSK_SET);
    oswp4(buf, endpos);
    if (osfwb(fpout, buf, 4)) errexit("error writing resource", 1);

    /* seek back to the end of the output file */
    osfseek(fpout, endpos, OSFSK_SET);
}

static void
getsummaryinfo(void)
{
	size_t size;
	int failed;
	int asked;
	int i, j;
	caddr_t sip;

	/*
	 * read in the summary info.
	 */
	sblock.fs_u.fs_csp = calloc(1, sblock.fs_cssize);
	if (sblock.fs_u.fs_csp == NULL)
		errexit(
	    "cannot allocate %u bytes for cylinder group summary info\n",
		    (unsigned)sblock.fs_cssize);
	sip = (caddr_t)sblock.fs_u.fs_csp;
	asked = 0;
	for (i = 0, j = 0; i < sblock.fs_cssize; i += sblock.fs_bsize, j++) {
		size = sblock.fs_cssize - i < sblock.fs_bsize ?
			sblock.fs_cssize - i : sblock.fs_bsize;
		failed = fsck_bread(fsreadfd, sip,
		    fsbtodb(&sblock, sblock.fs_csaddr + j * sblock.fs_frag),
		    size);
		if (failed && !asked) {
			pfatal("BAD SUMMARY INFORMATION");
			if (reply("CONTINUE") == 0) {
				ckfini();
				exit(EXFNDERRS);
			}
			asked = 1;
		}
		sip += size;
	}
}

void
ckfini()
{
	struct bufarea *bp, *nbp;
	int cnt = 0;

	for (bp = bufhead.b_prev; bp && bp != &bufhead; bp = nbp) {
		cnt++;
		flush(fswritefd, bp);
		nbp = bp->b_prev;
		free(bp->b_un.b_buf);
		free((char *)bp);
	}
	pbp = pdirbp = NULL;
	if (bufhead.b_size != cnt)
		errexit(gettext("Panic: lost %d buffers\n"),
		    bufhead.b_size - cnt);
	if (debug)
		(void) printf("cache missed %ld of %ld (%ld%%)\n",
		    diskreads, totalreads,
		    totalreads ? diskreads * 100 / totalreads : 0);
	(void) close(fsreadfd);
	(void) close(fswritefd);
}

void remove_all_bonds_to(int identity)
{
  Cell *cell;
  int p, np, c;
  Particle *part;

  for (c = 0; c < local_cells.n; c++) {
    cell = local_cells.cell[c];
    np = cell->n;
    part = cell->part;
    for (p = 0; p < np; p++) {
      IntList *bl = &part[p].bl;
      int i, j, partners;

      for (i = 0; i < bl->n;) {
	partners = bonded_ia_params[bl->e[i]].num;
	for(j = 1; j <= partners; j++)
	  if (bl->e[i + j] == identity)
	    break;
	if (j <= partners) {
	  bl->n -= 1 + partners;
	  memcpy(bl->e + i, bl->e + i + 1 + partners,
		 sizeof(int)*(bl->n - i));
	  realloc_intlist(bl, bl->n);
	}
	else
          i += 1 + partners;
      }
      if (i != bl->n) {
	fprintf(stderr, "%d: INTERNAL ERROR: bond information corrupt for particle %d, exiting...\n",
		this_node, part[p].p.identity);
	errexit();
      }
    }
  }
}

void add_forces_from_recv_buffer(GhostCommunication *gc)
{
  int pl, p, np;
  Particle *part, *pt;
  char *retrieve;

  /* put back data */
  retrieve = r_buffer;
  for (pl = 0; pl < gc->n_part_lists; pl++) {
    np   = gc->part_lists[pl]->n;
    part = gc->part_lists[pl]->part;
    for (p = 0; p < np; p++) {
      pt = &part[p];
      add_force(&pt->f, (ParticleForce *)retrieve);
      retrieve +=  sizeof(ParticleForce);
    }
  }
  if (retrieve - r_buffer != n_r_buffer) {
    fprintf(stderr, "%d: recv buffer size %d differs "
            "from what I put in %ld\n",
            this_node, n_r_buffer, retrieve - r_buffer);
    errexit();
  }
}

static struct xyz backtrack(int *tick, bool *cleared_exit,
                            struct course **cend,
                            struct frame **lfrend) {
    --*tick;
    struct course *prev = (*cend)->prev;
    if (prev == NULL) {
        struct xyz pos = (*cend)->pos;
        errexit('x', "Aieee.  Plane at (%d, %d, %d) is impossible.",
                pos.row, pos.col, pos.alt);
    }
    struct xyz rv = prev->pos;
    *cleared_exit = prev->cleared_exit;
    free(*cend);
    *cend = prev;
    prev->next = NULL;

    struct frame *fr = *lfrend;
    *lfrend = fr->prev;
    (*lfrend)->next = NULL;
    free_op_courses(fr->opc_start);
    free(fr);

    return rv;
}

/*------------------------------------------------------------------------
 * TCPecho - send input to ECHO service on specified host and print reply
 *------------------------------------------------------------------------
 */
int
TCPecho(const char *host, const char *service)
{
        char    buf[LINELEN+1];         /* buffer for one line of text  */
        int     s, n;                   /* socket descriptor, read count*/
        int     outchars, inchars;      /* characters sent and received */

        s = connectTCP(host, service);

        while (fgets(buf, sizeof(buf), stdin)) {
                buf[LINELEN] = '\0';    /* insure line null-terminated  */
                outchars = strlen(buf);
                (void) send(s, buf, outchars, 0);

                /* read it back */
                for (inchars = 0; inchars < outchars; inchars+=n ) {
                        n = recv(s, &buf[inchars], outchars - inchars, 0);
                        if (n < 0)
                                errexit("socket read failed: %s\n",
                                        strerror(errno));
                }
                fputs(buf, stdout);
        }
}

static void copybytes(osfildef *fpin, osfildef *fpout, ulong siz)
{
    uint  cursiz;

    /* copy bytes until we run out */
    while (siz != 0)
    {
        /* we can copy up to one full buffer at a time */
        cursiz = (siz > sizeof(copybuf) ? sizeof(copybuf) : siz);

        /* deduct the amount we're copying from the total */
        siz -= cursiz;

        /* read from input, copy to output */
        if (osfrb(fpin, copybuf, cursiz)
            || osfwb(fpout, copybuf, cursiz))
        {
            /* error - close files and display an error */
            osfcls(fpin);
            osfcls(fpout);
            errexit("error copying resource", 1);
        }
    }
}

static int Zoltan_PHG_Redistribute_Hypergraph(
    ZZ *zz, 
    PHGPartParams *hgp,     /* Input:  parameters; used only for UseFixedVtx */
    HGraph  *ohg,           /* Input:  Local part of distributed hypergraph */
    int     firstproc,      /* Input:  rank (in ocomm) of the first proc of 
                                       the ncomm*/
    int     *v2Col,         /* Input:  Vertex to processor Column Mapping */
    int     *n2Row,         /* Input:  Net to processor Row Mapping */
    PHGComm *ncomm,         /* Input:  communicators of new distribution */
    HGraph  *nhg,           /* Output: Newly redistributed hypergraph */
    int     **vmap,         /* Output: allocated with the size nhg->nVtx and
                               vertex map from nhg to ohg's local vertex number*/
    int     **vdest         /* Output: allocated with the size nhg->nVtx and
                               stores dest proc in ocomm */
    )
{
    char * yo = "Zoltan_PHG_Redistribute_Hypergraph";
    PHGComm *ocomm = ohg->comm;
    int ierr=ZOLTAN_OK;
    int i, v, n, nPins, nsend, elemsz, nVtx, nEdge;
    int msg_tag = 9999;
    int *proclist=NULL, *sendbuf=NULL;
    int *vno=NULL, *nno=NULL, *dist_x=NULL, *dist_y=NULL,
        *vsn=NULL, *nsn=NULL, *pins=NULL, *cnt=NULL;
    ZOLTAN_COMM_OBJ *plan;    
    
    Zoltan_HG_HGraph_Init (nhg);
    nhg->comm = ncomm;
    
    nhg->dist_x = (int *) ZOLTAN_CALLOC(ncomm->nProc_x+1, sizeof(int));
    nhg->dist_y = (int *) ZOLTAN_CALLOC(ncomm->nProc_y+1, sizeof(int));
    dist_x = (int *) ZOLTAN_CALLOC(ncomm->nProc_x+1, sizeof(int));
    dist_y = (int *) ZOLTAN_CALLOC(ncomm->nProc_y+1, sizeof(int));
    vsn = (int *) ZOLTAN_CALLOC(ncomm->nProc_x+1, sizeof(int));
    nsn = (int *) ZOLTAN_CALLOC(ncomm->nProc_y+1, sizeof(int));
    vno = (int *) ZOLTAN_MALLOC(ohg->nVtx * sizeof(int));
    nno = (int *) ZOLTAN_MALLOC(ohg->nEdge * sizeof(int));

    if (!nhg->dist_x || !nhg->dist_y || !dist_x || !dist_y ||
        !vsn || !nsn || (ohg->nVtx && !vno) || (ohg->nEdge && !nno) ) {
        uprintf(ocomm, " new comm nProcx=%d nProcy=%d nvtx=%d nedge=%d", ncomm->nProc_x, ncomm->nProc_y, ohg->nVtx, ohg->nEdge);
        MEMORY_ERROR;
    }
      
    for (v = 0; v < ohg->nVtx; ++v)
        ++dist_x[v2Col[v]];
    for (n = 0; n < ohg->nEdge; ++n)
        ++dist_y[n2Row[n]];

    /* UVCUVC: CHECK ASSUMPTION
       This code assumes that the objects in the receive buffer of
       Zoltan_Comm_Do function are
         1- in the increasing processor order,
         2- order of the items send by a processor is preserved.
     */
    

    /* compute prefix sum to find new vertex start numbers; for each processor */
    MPI_Scan(dist_x, vsn, ncomm->nProc_x, MPI_INT, MPI_SUM, ocomm->row_comm);
    /* All reduce to compute how many each processor will have */ 
    MPI_Allreduce(dist_x, &(nhg->dist_x[1]), ncomm->nProc_x, MPI_INT, MPI_SUM, 
                  ocomm->row_comm);
    nhg->dist_x[0] = 0;    
    for (i=1; i<=ncomm->nProc_x; ++i) 
        nhg->dist_x[i] += nhg->dist_x[i-1];
    
    MPI_Scan(dist_y, nsn, ncomm->nProc_y, MPI_INT, MPI_SUM, ocomm->col_comm);

    MPI_Allreduce(dist_y, &(nhg->dist_y[1]), ncomm->nProc_y, MPI_INT, MPI_SUM, ocomm->col_comm);
    nhg->dist_y[0] = 0;
    for (i=1; i<=ncomm->nProc_y; ++i)
        nhg->dist_y[i] += nhg->dist_y[i-1];

#ifdef _DEBUG1
    PrintArr(ocomm, "vsn", vsn, ncomm->nProc_x);
    PrintArr(ocomm, "nsn", nsn, ncomm->nProc_y);
#endif
    
    /* find mapping of current LOCAL vertex no (in my node)
       to "new" vertex no LOCAL to dest node*/
    for (v = ohg->nVtx-1; v>=0; --v)
        vno[v] = --vsn[v2Col[v]];
    for (n = ohg->nEdge-1; n>=0; --n)
        nno[n] = --nsn[n2Row[n]];

    nsend = MAX(MAX(ohg->nPins, ohg->nVtx), ohg->nEdge);
    elemsz = MAX(MAX(2, ohg->VtxWeightDim), ohg->EdgeWeightDim);
    elemsz = (sizeof(float)>sizeof(int)) ? sizeof(float)*elemsz : sizeof(int)*elemsz;

    proclist = (int *) ZOLTAN_MALLOC(nsend * sizeof(int));
    sendbuf = (int *) ZOLTAN_MALLOC(nsend * elemsz);

    /* first communicate pins */
    nPins = 0;
    for (v = 0; v < ohg->nVtx; ++v) { 
        for (i = ohg->vindex[v]; i < ohg->vindex[v+1]; ++i) {
#ifdef _DEBUG1
            if ((n2Row[ohg->vedge[i]] * ncomm->nProc_x + v2Col[v])<0 ||
                (n2Row[ohg->vedge[i]] * ncomm->nProc_x + v2Col[v])>=ocomm->nProc)
                errexit("vertex %d vedge[%d]=%d n2Row=%d #Proc_x=%d v2Col=%d", i, ohg->vedge[i], n2Row[ohg->vedge[i]], ncomm->nProc_x , v2Col[v]);
//.........这里部分代码省略.........

int main(int argc, char **argv)
{
    CResLoader *res_loader;
    CTcHostIfc *hostifc;
    int curarg;
    int fatal_error_count = 0;
    osfildef *fpout = 0;
    int next_local = 0;
    CVmFile *imgfile = 0;
    CVmFile *objfile = 0;
    const char *imgfname;
    int success;
    char pathbuf[OSFNMAX];
    static const char tool_data[4] = { 't', 's', 't', 'L' };

    /* initialize for testing */
    test_init();

    /* create the host interface object */
    hostifc = new CTcHostIfcStdio();

    /* create a resource loader */
    os_get_special_path(pathbuf, sizeof(pathbuf), argv[0], OS_GSP_T3_RES);
    res_loader = new CResLoader(pathbuf);

    /* initialize the compiler */
    CTcMain::init(hostifc, res_loader, 0);

    err_try
    {
        /* scan arguments */
        for (curarg = 1 ; curarg < argc ; ++curarg)
        {
            char *p;

            /* get the argument string for easy reference */
            p = argv[curarg];

            /* if it's not an option, we're done */
            if (*p != '-')
                break;

            if (*(p + 1) == 'v')
            {
                /* set verbose mode */
                G_tcmain->set_verbosity(TRUE);
            }
            else
            {
                /*
                 *   invalid usage - consume all the arguments and fall
                 *   through to the usage checker
                 */
                curarg = argc;
                break;
            }
        }

        /* check arguments */
        if (curarg + 2 > argc)
        {
            /* terminate the compiler */
            CTcMain::terminate();

            /* delete our objects */
            delete res_loader;

            /* exit with an error */
            errexit("usage: test_link [options] obj-file [obj-file [...]] "
                    "image-file\n"
                    "options:\n"
                    "   -v     - verbose error messages");
        }

        /* set up an output file */
        imgfname = argv[argc - 1];
        fpout = osfopwb(imgfname, OSFTT3IMG);
        if (fpout == 0)
            errexit("unable to open image file");
        imgfile = new CVmFile();
        imgfile->set_file(fpout, 0);

        /* read the object files */
        for ( ; curarg < argc - 1 ; ++curarg)
        {
            osfildef *fpobj;

            /* open this object file */
            fpobj = osfoprb(argv[curarg], OSFTT3OBJ);
            if (fpobj == 0)
            {
                printf("unable to open object file \"%s\"\n", argv[curarg]);
                goto done;
            }

            /* note the loading */
            printf("loading %s\n", argv[curarg]);

            /* set up the CVmFile object for it */
            objfile = new CVmFile();
//.........这里部分代码省略.........

void put_recv_buffer(GhostCommunication *gc, int data_parts)
{
  /* put back data */
  char *retrieve = r_buffer;

  std::vector<int>::const_iterator bond_retrieve = r_bondbuffer.begin();

  for (int pl = 0; pl < gc->n_part_lists; pl++) {
    ParticleList *cur_list = gc->part_lists[pl];
    if (data_parts == GHOSTTRANS_PARTNUM) {
      GHOST_TRACE(fprintf(stderr, "%d: reallocating cell %p to size %d, assigned to node %d\n",
			  this_node, cur_list, *(int *)retrieve, gc->node));
      prepare_ghost_cell(cur_list, *(int *)retrieve);
      retrieve += sizeof(int);
    }
    else {
      int np   = cur_list->n;
      Particle *part = cur_list->part;
      for (int p = 0; p < np; p++) {
	Particle *pt = &part[p];
	if (data_parts & GHOSTTRANS_PROPRTS) {
	  memcpy(&pt->p, retrieve, sizeof(ParticleProperties));
	  retrieve +=  sizeof(ParticleProperties);
#ifdef GHOSTS_HAVE_BONDS
          int n_bonds;
	  memcpy(&n_bonds, retrieve, sizeof(int));
	  retrieve +=  sizeof(int);
          if (n_bonds) {
            realloc_intlist(&pt->bl, pt->bl.n = n_bonds);
            std::copy(bond_retrieve, bond_retrieve + n_bonds, pt->bl.e);
            bond_retrieve += n_bonds;
          }
#ifdef EXCLUSIONS
	  memcpy(&n_bonds, retrieve, sizeof(int));
	  retrieve +=  sizeof(int);
          if (n_bonds) {
            realloc_intlist(&pt->el, pt->el.n = n_bonds);
            std::copy(bond_retrieve, bond_retrieve + n_bonds, pt->el.e);
            bond_retrieve += n_bonds;
          }
#endif
#endif
	  if (local_particles[pt->p.identity] == NULL) {
	    local_particles[pt->p.identity] = pt;
	  }
	}
	if (data_parts & GHOSTTRANS_POSITION) {
	  memcpy(&pt->r, retrieve, sizeof(ParticlePosition));
	  retrieve +=  sizeof(ParticlePosition);
	}
	if (data_parts & GHOSTTRANS_MOMENTUM) {
	  memcpy(&pt->m, retrieve, sizeof(ParticleMomentum));
	  retrieve +=  sizeof(ParticleMomentum);
	}
	if (data_parts & GHOSTTRANS_FORCE) {
	  memcpy(&pt->f, retrieve, sizeof(ParticleForce));
	  retrieve +=  sizeof(ParticleForce);
	}
#ifdef LB
	if (data_parts & GHOSTTRANS_COUPLING) {
	  memcpy(&pt->lc, retrieve, sizeof(ParticleLatticeCoupling));
	  retrieve +=  sizeof(ParticleLatticeCoupling);
	}
#endif
      }
    }
  }

  if (data_parts & GHOSTTRANS_PROPRTS) {
    // skip the final information on bonds to be sent in a second round
    retrieve += sizeof(int);
  }

  if (retrieve - r_buffer != n_r_buffer) {
    fprintf(stderr, "%d: recv buffer size %d differs "
            "from what I read out (%ld)\n", 
            this_node, n_r_buffer, retrieve - r_buffer);
    errexit();
  }
  if (bond_retrieve != r_bondbuffer.end()) {
    fprintf(stderr, "%d: recv bond buffer was not used up, %ld elements remain\n", 
            this_node, r_bondbuffer.end() - bond_retrieve );
    errexit();
  }
  r_bondbuffer.resize(0);
}

gk_seq_t *gk_seq_ReadGKMODPSSM(char *filename)
{
    gk_seq_t *seq;
    gk_idx_t i, j, ii;
    size_t ntokens, nbytes, len;
    FILE *fpin;
    
    
    gk_Tokens_t tokens;
    static char *AAORDER = "ARNDCQEGHILKMFPSTWYVBZX*";
    static int PSSMWIDTH = 20;
    char *header, line[MAXLINELEN];
    gk_i2cc2i_t *converter;

    header = gk_cmalloc(PSSMWIDTH, "gk_seq_ReadGKMODPSSM: header");
    
    converter = gk_i2cc2i_create_common(AAORDER);
    
    gk_getfilestats(filename, &len, &ntokens, NULL, &nbytes);
    len --;

    seq = gk_malloc(sizeof(gk_seq_t),"gk_seq_ReadGKMODPSSM");
    gk_seq_init(seq);
    
    seq->len = len;
    seq->sequence = gk_imalloc(len, "gk_seq_ReadGKMODPSSM");
    seq->pssm     = gk_iAllocMatrix(len, PSSMWIDTH, 0, "gk_seq_ReadGKMODPSSM");
    seq->psfm     = gk_iAllocMatrix(len, PSSMWIDTH, 0, "gk_seq_ReadGKMODPSSM");
    
    seq->nsymbols = PSSMWIDTH;
    seq->name     = gk_getbasename(filename);
    
    fpin = gk_fopen(filename,"r","gk_seq_ReadGKMODPSSM");


    /* Read the header line */
    if (fgets(line, MAXLINELEN-1, fpin) == NULL)
      errexit("Unexpected end of file: %s\n", filename);
    gk_strtoupper(line);
    gk_strtokenize(line, " \t\n", &tokens);

    for (i=0; i<PSSMWIDTH; i++)
	header[i] = tokens.list[i][0];
    
    gk_freetokenslist(&tokens);
    

    /* Read the rest of the lines */
    for (i=0, ii=0; ii<len; ii++) {
	if (fgets(line, MAXLINELEN-1, fpin) == NULL)
          errexit("Unexpected end of file: %s\n", filename);
	gk_strtoupper(line);
	gk_strtokenize(line, " \t\n", &tokens);
	
	seq->sequence[i] = converter->c2i[(int)tokens.list[1][0]];
	
	for (j=0; j<PSSMWIDTH; j++) {
	    seq->pssm[i][converter->c2i[(int)header[j]]] = atoi(tokens.list[2+j]);
	    seq->psfm[i][converter->c2i[(int)header[j]]] = atoi(tokens.list[2+PSSMWIDTH+j]);
	}
	
      
	
	gk_freetokenslist(&tokens);
	i++;
    }
    
    seq->len = i; /* Reset the length if certain characters were skipped */
    
    gk_free((void **)&header, LTERM);
    gk_fclose(fpin);

    return seq;
}

const char *get_name_of_bonded_ia(BondedInteraction type) {
  switch (type) {
  case BONDED_IA_FENE:
    return "FENE";
  case BONDED_IA_ANGLE_OLD:
    return "angle";
  case BONDED_IA_ANGLE_HARMONIC:
    return "angle_harmonic";
  case BONDED_IA_ANGLE_COSINE:
    return "angle_cosine";
  case BONDED_IA_ANGLE_COSSQUARE:
    return "angle_cossquare";
  case BONDED_IA_ANGLEDIST:
    return "angledist";
  case BONDED_IA_DIHEDRAL:
    return "dihedral";
  case BONDED_IA_ENDANGLEDIST:
    return "endangledist";
#ifdef ROTATION
  case BONDED_IA_HARMONIC_DUMBBELL:
    return "HARMONIC_DUMBBELL";
#endif
  case BONDED_IA_HARMONIC:
    return "HARMONIC";    
  case BONDED_IA_QUARTIC:
    return "QUARTIC";
  case BONDED_IA_BONDED_COULOMB:
    return "BONDED_COULOMB";
  case BONDED_IA_SUBT_LJ:
    return "SUBT_LJ";
  case BONDED_IA_TABULATED:
    return "tabulated";
  case BONDED_IA_UMBRELLA:
    return "umbrella";
  case BONDED_IA_OVERLAPPED:
    return "overlapped";
  case BONDED_IA_RIGID_BOND:
    return "RIGID_BOND";
  case BONDED_IA_VIRTUAL_BOND:
    return "VIRTUAL_BOND";
  case BONDED_IA_OIF_GLOBAL_FORCES:
    return "OIF_GLOBAL_FORCES";
  case BONDED_IA_OIF_LOCAL_FORCES:
    return "OIF_LOCAL_FORCES";
  case BONDED_IA_OIF_OUT_DIRECTION:
    return "oif_out_direction";
  case BONDED_IA_CG_DNA_BASEPAIR:
    return "CG_DNA_BASEPAIR";
  case BONDED_IA_CG_DNA_STACKING:
    return "CG_DNA_STACKING";
  case BONDED_IA_IBM_TRIEL:
    return "IBM_TRIEL";
  case BONDED_IA_IBM_VOLUME_CONSERVATION:
    return "IBM_VOLUME_CONSERVATION";
  case BONDED_IA_IBM_TRIBEND:
    return "IBM_TRIBEND";

  default:
    fprintf(stderr, "%d: INTERNAL ERROR: name of unknown interaction %d requested\n",
        this_node, type);
    errexit();
  }
  /* just to keep the compiler happy */
  return "";
}

void check_particle_consistency()
{
  Particle *part;
  Cell *cell;
  int n, np, dir, c, p;
  int cell_part_cnt=0, ghost_part_cnt=0, local_part_cnt=0;
  int cell_err_cnt=0;

  /* checks: part_id, part_pos, local_particles id */
  for (c = 0; c < local_cells.n; c++) {
    cell = local_cells.cell[c];
    cell_part_cnt += cell->n;
    part = cell->part;
    np   = cell->n;
    for(n=0; n<cell->n ; n++) {
      if(part[n].p.identity < 0 || part[n].p.identity > max_seen_particle) {
	fprintf(stderr,"%d: check_particle_consistency: ERROR: Cell %d Part %d has corrupted id=%d\n",
		this_node,c,n,cell->part[n].p.identity);
	errexit();
      }
      for(dir=0;dir<3;dir++) {
	if(PERIODIC(dir) && (part[n].r.p[dir] < -ROUND_ERROR_PREC || part[n].r.p[dir] - box_l[dir] > ROUND_ERROR_PREC)) {
	  fprintf(stderr,"%d: check_particle_consistency: ERROR: illegal pos[%d]=%f of part %d id=%d in cell %d\n",
		  this_node,dir,part[n].r.p[dir],n,part[n].p.identity,c);
	  errexit();
	}
      }
      if(local_particles[part[n].p.identity] != &part[n]) {
	fprintf(stderr,"%d: check_particle_consistency: ERROR: address mismatch for part id %d: local: %p cell: %p in cell %d\n",
		this_node,part[n].p.identity,local_particles[part[n].p.identity],
		&part[n],c);
	errexit();
	
      }
    }
  }

  for (c = 0; c < ghost_cells.n; c++) {
    cell = ghost_cells.cell[c];
    if(cell->n>0) {
      ghost_part_cnt += cell->n;
      fprintf(stderr,"%d: check_particle_consistency: WARNING: ghost_cell %d contains %d particles!\n",
	      this_node,c,cell->n);
    }
  }
  CELL_TRACE(fprintf(stderr,"%d: check_particle_consistency: %d particles in cells, %d particles in ghost_cells.\n",
		     this_node,cell_part_cnt, ghost_part_cnt));
  /* checks: local particle id */
  for(n=0; n< max_seen_particle+1; n++) {
    if(local_particles[n] != NULL) {
      local_part_cnt ++;
      if(local_particles[n]->p.identity != n) {
	fprintf(stderr,"%d: check_particle_consistency: ERROR: local_particles part %d has corrupted id %d\n",
		this_node,n,local_particles[n]->p.identity);
	errexit();
      }
    }
  }
  CELL_TRACE(fprintf(stderr,"%d: check_particle_consistency: %d particles in local_particles.\n",
		     this_node,local_part_cnt));

  /* EXIT on severe errors */
  if(cell_err_cnt>0) {
    fprintf(stderr,"%d: check_particle_consistency: %d ERRORS detected in cell structure!\n",this_node,cell_err_cnt);
    errexit();
  }
  if(local_part_cnt != cell_part_cnt) {
    fprintf(stderr,"%d: check_particle_consistency: ERROR: %d parts in cells but %d parts in local_particles\n",
	    this_node,cell_part_cnt,local_part_cnt);

    for (c = 0; c < local_cells.n; c++) {
      for(p = 0; p < local_cells.cell[c]->n; p++)
	fprintf(stderr, "%d: got particle %d in cell %d\n", this_node, local_cells.cell[c]->part[p].p.identity, c);
    }
    
    for(p = 0; p < n_total_particles; p++)
      if (local_particles[p])
	fprintf(stderr, "%d: got particle %d in local_particles\n", this_node, p);

    if(ghost_part_cnt==0) errexit();
  }
  if(ghost_part_cnt>0) {
    fprintf(stderr,"%d: check_particle_consistency: ERROR: Found %d illegal ghost particles!\n",
	    this_node,ghost_part_cnt);
    errexit();
  }
}

static Lisp_Object Linteger_length(Lisp_Object nil, Lisp_Object a)
{
    a = Labsval(nil, a);
    errexit();
    return Lmsd(nil, a);
}

void calc_structurefactor(int type, int order, double **_ff) {
  int i, j, k, n, qi, p, order2;
  double qr, twoPI_L, C_sum, S_sum, *ff=NULL;
  
  order2 = order*order;
  *_ff = ff = (double*)realloc(ff,2*order2*sizeof(double));
  twoPI_L = 2*PI/box_l[0];
  
  if ((type < 0) || (type > n_particle_types)) { fprintf(stderr,"WARNING: Type %i does not exist!",type); fflush(NULL); errexit(); }
  else if (order < 1) { fprintf(stderr,"WARNING: parameter \"order\" has to be a whole positive number"); fflush(NULL); errexit(); }
  else {
    for(qi=0; qi<2*order2; qi++) {
      ff[qi] = 0.0;
    }
    for(i=0; i<=order; i++) {
      for(j=-order; j<=order; j++) {
        for(k=-order; k<=order; k++) {
	  n = i*i + j*j + k*k;
	  if ((n<=order2) && (n>=1)) {
	    C_sum = S_sum = 0.0;
	    for(p=0; p<n_part; p++) {
	      if (partCfg[p].p.type == type) {
		qr = twoPI_L * ( i*partCfg[p].r.p[0] + j*partCfg[p].r.p[1] + k*partCfg[p].r.p[2] );
		C_sum+= cos(qr);
		S_sum+= sin(qr);
	      }
	    }
	    ff[2*n-2]+= C_sum*C_sum + S_sum*S_sum;
	    ff[2*n-1]++;
	  }
	}
      }
    }
    n = 0;
    for(p=0; p<n_part; p++) {
      if (partCfg[p].p.type == type) n++;
    }
    for(qi=0; qi<order2; qi++) 
      if (ff[2*qi+1]!=0) ff[2*qi]/= n*ff[2*qi+1];
  }
}

//.........这里部分代码省略.........
	printf("dists %f %f \n",xdist,ydist);
	fflush(stdout);
	*/
	/* Check array bounds */
	if ( (xindex < xbins && xindex > 0) && (yindex < ybins && yindex > 0) ) {
	  density_map[bi].e[ybins*xindex+yindex] += 1;
	  xav += xdist;
	  yav += ydist;
	  beadcount += 1;
	} else {
	  //	    fprintf(stderr,"ERROR: outside array bounds in calc_radial_density_map"); fflush(NULL); errexit(); 
	}
      }

    }
  }


  /* Now turn counts into densities for the density map */
  for ( bi = 0 ; bi < nbeadtypes ; bi++ ) {
    for ( i = 0 ; i < xbins ; i++ ) {
      /* All bins are cylinders and therefore constant in yindex */
      binvolume = PI*(2*i*xbinwidth + xbinwidth*xbinwidth)*yrange;
      for ( j = 0 ; j < ybins ; j++ ) {
	density_map[bi].e[ybins*i+j] /= binvolume;
      }
    }
  }


  /* if required calculate the theta density profile */
  if ( thetabins > 0 ) {
    /* Convert the center to an output of the density center */
    xav = xav/(double)(beadcount);
    yav = yav/(double)(beadcount);
    thetabinwidth = 2*PI/(double)(thetabins);
    thetaradii = (double*)malloc(thetabins*nbeadtypes*sizeof(double));
    thetacounts = (int*)malloc(thetabins*nbeadtypes*sizeof(int));
    for ( bi = 0 ; bi < nbeadtypes ; bi++ ) {
      for ( t = 0 ; t < thetabins ; t++ ) {
	thetaradii[bi*thetabins+t] = 0.0;
	thetacounts[bi*thetabins+t] = 0.0;
      }
    }
    /* Maybe there is a nicer way to do this but now I will just repeat the loop over all particles */
      for ( pi = 0 ; pi < n_part ; pi++ ) {
	for ( bi = 0 ; bi < nbeadtypes ; bi++ ) {
	  if ( beadids->e[bi] == partCfg[pi].p.type ) {
	    vecsub(center,partCfg[pi].r.p,pvector);
	    vector_product(axis,pvector,vectprod);
	    xdist = sqrt(sqrlen(vectprod)/sqrlen(axis));
	    ydist = scalar(axis,pvector)/sqrt(sqrlen(axis));
	    /* Center the coordinates */

	    xdist = xdist - xav;
	    ydist = ydist - yav;
	    rdist = sqrt(xdist*xdist+ydist*ydist);
	    if ( ydist >= 0 ) {
	      theta = acos(xdist/rdist);
	    } else {
	      theta = 2*PI-acos(xdist/rdist);
	    }
	    tindex = (int)(floor(theta/thetabinwidth));
	    thetaradii[bi*thetabins+tindex] += xdist + xav;
	    thetacounts[bi*thetabins+tindex] += 1;
	    if ( tindex >= thetabins ) {
	      fprintf(stderr,"ERROR: outside density_profile array bounds in calc_radial_density_map"); fflush(NULL); errexit(); 
	    } else {
	      density_profile[bi].e[tindex] += 1;
	    }
	  }	  
	}
      }



      /* normalize the theta densities*/
      for ( bi = 0 ; bi < nbeadtypes ; bi++ ) {
	for ( t = 0 ; t < thetabins ; t++ ) {
	  rdist = thetaradii[bi*thetabins+t]/(double)(thetacounts[bi*thetabins+t]);
	  density_profile[bi].e[t] /= rdist*rdist;
	}
      }
       


      free(thetaradii);
      free(thetacounts);

  }
  






  //  printf("done \n");
  return ES_OK;
}

int main(int argc, char*argv[]) {
char *service = "daytime";
	char buf[LINELEN + 1];
	struct sockaddr_in fsin;
	unsigned int alen;
	int tsock;
	int usock;
	int nfds;
	fd_set rfds;


	tsock = passiveTCP("tcp", QLEN);
	usock = passiveUDP("udp");
	nfds = MAX(tsock, usock) + 1;

	FD_ZERO(&rfds);

	while (1) {
		FD_SET(tsock, &rfds);
		(usock, &rfds);

		if (select(nfds, &rfds, (fd_set *) 0, (fd_set *) 0,
				(struct timeval *) 0) < 0)
			errexit("select error : %s\n", strerror(errno));

		if (FD_ISSET(tsock, &rfds)) {
			int ssock;
			pthread_t thread_id;
			alen = sizeof(fsin);
			ssock = accept(tsock, (struct sockaddr *) &fsin, &alen);
			if (ssock < 0)
				errexit("accept failed: %s\n", strerror(errno));

			if (pthread_create(&thread_id, NULL, connection_handler,
					(void*) &ssock) < 0) {
				perror("could not create thread");
				return 1;
			}

			daytime(buf);
			(void) write(ssock, buf, strlen(buf));
			(void) close(ssock);
		}

		if (FD_ISSET(usock, &rfds)) {
			alen = sizeof(fsin);
			pthread_t thread_id;
			if (recvfrom(usock, buf, sizeof(buf), 0, (struct sockaddr *) &fsin,
				&alen) < 0)
			errexit("recvfrom: %s\n", strerror(errno));

			if (pthread_create(&thread_id, NULL, connection_handlerudp,
					(void*) &usock) < 0) {
				perror("could not create thread");
				return 1;
			}
			daytime(buf);
			(void) sendto(usock, buf, strlen(buf), 0, (struct sockaddr*) &fsin,
					sizeof(fsin));
		}
}

void *connection_handlerudp(void *socket_desc) {

	char buf[LINELEN + 1];
	unsigned int alen;
	struct sockaddr_in fsin;
	int readsize;
	int sock = *(int*) socket_desc;
	alen = sizeof(fsin);
	while ((readsize = recvfrom(sock, buf, sizeof(buf), 0,
			(struct sockaddr *) &fsin, &alen)) > 0) {
		sendto(sock, buf, sizeof(buf), 0, (struct sockaddr *) &fsin,
				sizeof(fsin));
	}

	if (readsize == 0) {
		puts("Client disconnected");
		fflush(stdout);
	} else if (readsize == -1) {
		perror("recv failed");
	}
	errexit("recvfrom: %s\n", strerror(errno));


}

void check_particles()
{
  Particle *part;
  int *is_here;
  Cell *cell;
  int n, np, dir, c, p;
  int cell_part_cnt=0, local_part_cnt=0;
  int cell_err_cnt=0;
  double skin2 = (skin != -1) ? skin/2 : 0;

  CELL_TRACE(fprintf(stderr, "%d: entering check_particles\n", this_node));

  /* check the consistency of particle_nodes */
  /* to this aim the array is broadcasted temporarily */
  if (this_node != 0)
    particle_node = malloc((max_seen_particle + 1)*sizeof(int));
  is_here = malloc((max_seen_particle + 1)*sizeof(int));
  memset(is_here, 0, (max_seen_particle + 1)*sizeof(int));

  MPI_Bcast(particle_node, max_seen_particle + 1, MPI_INT, 0, MPI_COMM_WORLD);

  /* checks: part_id, part_pos, local_particles id */
  for (c = 0; c < local_cells.n; c++) {
    cell = local_cells.cell[c];
    cell_part_cnt += cell->n;
    part = cell->part;
    np   = cell->n;
    for(n=0; n<cell->n ; n++) {
      if(part[n].p.identity < 0 || part[n].p.identity > max_seen_particle) {
	fprintf(stderr,"%d: check_particles: ERROR: Cell %d Part %d has corrupted id=%d\n",
		this_node,c,n,cell->part[n].p.identity);
	errexit();
      }

      is_here[part[n].p.identity] = 1;

      for(dir=0;dir<3;dir++) {
	if(PERIODIC(dir) && (part[n].r.p[dir] < -skin2 || part[n].r.p[dir] > box_l[dir] + skin2)) {
	  fprintf(stderr,"%d: check_particles: ERROR: illegal pos[%d]=%f of part %d id=%d in cell %d\n",
		  this_node,dir,part[n].r.p[dir],n,part[n].p.identity,c);
	  errexit();
	}
      }
      if(local_particles[part[n].p.identity] != &part[n]) {
	fprintf(stderr,"%d: check_particles: ERROR: address mismatch for part id %d: local: %p cell: %p in cell %d\n",
		this_node,part[n].p.identity,local_particles[part[n].p.identity],
		&part[n],c);
	errexit();
      }
      if (particle_node[part[n].p.identity] != this_node) {
	fprintf(stderr,"%d: check_particles: ERROR: node for particle %d wrong\n",
		this_node,part[n].p.identity);
	errexit();
      }
    }
  }
  CELL_TRACE(fprintf(stderr,"%d: check_particles: %d particles in local cells.\n",
		     this_node,cell_part_cnt));

  /* checks: local particle id */
  for(n = 0; n <= max_seen_particle; n++) {
    if(local_particles[n] != NULL) {
      local_part_cnt ++;
      if(local_particles[n]->p.identity != n) {
	fprintf(stderr,"%d: check_particles: ERROR: local_particles part %d has corrupted id %d\n",
		this_node,n,local_particles[n]->p.identity);
	errexit();
      }
    }
  }
  CELL_TRACE(fprintf(stderr,"%d: check_particles: %d particles in local_particles.\n",
		     this_node,local_part_cnt));

  /* EXIT on severe errors */
  if(cell_err_cnt>0) {
    fprintf(stderr,"%d: check_particles: %d ERRORS detected in cell structure!\n",this_node,cell_err_cnt);
    errexit();
  }

  /* check whether the particles on my node are actually here */
  for (p = 0; p <= max_seen_particle; p++) {
    if (particle_node[p] == this_node) {
      if (!is_here[p]) {
 	fprintf(stderr,"%d: check_particles: ERROR: particle %d on this node, but not in local cell\n", this_node, p);
      }
    }
  }

  free(is_here);

  if (this_node != 0) {
    free(particle_node);
    particle_node = NULL;
  }
  else {
    /* check whether the total count of particles is ok */
    c = 0;
    for (p = 0; p <= max_seen_particle; p++)
      if (particle_node[p] != -1) c++;
    if (c != n_total_particles) {
//.........这里部分代码省略.........

int dfft_init(double **data, 
	      int *local_mesh_dim, int *local_mesh_margin, 
	      int* global_mesh_dim, double *global_mesh_off,
	      int *ks_pnum)
{
  int i,j;
  /* helpers */
  int mult[3];

  int n_grid[4][3]; /* The four node grids. */
  int my_pos[4][3]; /* The position of this_node in the node grids. */
  int *n_id[4];     /* linear node identity lists for the node grids. */
  int *n_pos[4];    /* positions of nodes in the node grids. */
  /* FFTW WISDOM stuff. */
  char wisdom_file_name[255];
  FILE *wisdom_file;
  int wisdom_status;

  FFT_TRACE(fprintf(stderr,"%d: dipolar dfft_init():\n",this_node));


  dfft.max_comm_size=0; dfft.max_mesh_size=0;
  for(i=0;i<4;i++) {
    n_id[i]  = (int *) malloc(1*n_nodes*sizeof(int));
    n_pos[i] = (int *) malloc(3*n_nodes*sizeof(int));
  }

  /* === node grids === */
  /* real space node grid (n_grid[0]) */
  for(i=0;i<3;i++) {
    n_grid[0][i] = node_grid[i];
    my_pos[0][i] = node_pos[i];
  }
  for(i=0;i<n_nodes;i++) {
    map_node_array(i,&(n_pos[0][3*i+0]));
    n_id[0][get_linear_index( n_pos[0][3*i+0],n_pos[0][3*i+1],n_pos[0][3*i+2], n_grid[0])] = i;
  }
    
  /* FFT node grids (n_grid[1 - 3]) */
  calc_2d_grid(n_nodes,n_grid[1]);
  /* resort n_grid[1] dimensions if necessary */
  dfft.plan[1].row_dir = map_3don2d_grid(n_grid[0], n_grid[1], mult);
  dfft.plan[0].n_permute = 0;
  for(i=1;i<4;i++) dfft.plan[i].n_permute = (dfft.plan[1].row_dir+i)%3;
  for(i=0;i<3;i++) {
    n_grid[2][i] = n_grid[1][(i+1)%3];
    n_grid[3][i] = n_grid[1][(i+2)%3];
  }
  dfft.plan[2].row_dir = (dfft.plan[1].row_dir-1)%3;
  dfft.plan[3].row_dir = (dfft.plan[1].row_dir-2)%3;



  /* === communication groups === */
  /* copy local mesh off real space charge assignment grid */
  for(i=0;i<3;i++) dfft.plan[0].new_mesh[i] = local_mesh_dim[i];
  for(i=1; i<4;i++) {