mulmats.c

很经典的电磁计算（电容计算）软件 MIT90年代开发

#if OFF == ON			/* OFF == OFF produces the M2M error!! */
    for(i=0; i < 8; i++) multimats[i] = NULL;
#endif

  numterms = multerms(order);

  if(sys->depth < 2) {
    /* fprintf(stdout, "\nWarning: no multipole acceleration\n");*/
    return;	/* return if upward pass not possible */
  }

/* Handle the lowest level cubes first (set up Q2M's). */
  for(nextc=sys->multilist[sys->depth]; nextc != NULL; nextc = nextc->mnext) {
    nextc->multisize = numterms;
    CALLOC(nextc->multi, numterms, double, ON, AMSC);
    CALLOC(nextc->upmats, 1, double**, ON, AMSC);
    nextc->upmats[0] = mulQ2Multi(nextc->chgs, nextc->nbr_is_dummy[0],
				  nextc->upnumeles[0],
				  nextc->x, nextc->y, nextc->z, order);

#if DISSYN == ON
    multicnt[nextc->level]++;
#endif

#if DMTCNT == ON
    Q2Mcnt[nextc->level][nextc->level]++;
#endif

  }

  if(sys->locallist[sys->depth] == NULL
     && sys->multilist[sys->depth] == NULL) {
    fprintf(stdout, "No expansions at level %d (lowest)\n", sys->depth);
    /*if(sys->depth < 3) 
	fprintf(stdout, " (Warning: no multipole acceleration)\n");*/
  }
  else if(sys->locallist[sys->depth] == NULL) {
    fprintf(stdout, "No local expansions at level %d (lowest)\n", sys->depth);
  }
  else if(sys->multilist[sys->depth] == NULL) {
    fprintf(stdout, "No multipole expansions at level %d (lowest)\n", 
	    sys->depth); 
    /*if(sys->depth < 3) 
	fprintf(stdout, " (Warning: no multipole acceleration)\n");*/
  }

/* Allocate the vectors and matrices for the cubes. */
  /* no multipoles over root cube or its kids (would not be used if made) */
  for(depth = (sys->depth - 1); depth > 1; depth--) {
  /* set up M2M's and Q2M's to compute the multipoles needed for this level */
    if(sys->locallist[depth] == NULL && sys->multilist[depth] == NULL) {
      fprintf(stdout, "No expansions at level %d\n", depth);
      /*if(depth < 3) fprintf(stdout, " (Warning: no multipole acceleration)\n");
      else fprintf(stdout, "\n");*/
    }
    else if(sys->locallist[depth] == NULL) {
      fprintf(stdout, "No local expansions at level %d\n", depth);
    }
    else if(sys->multilist[depth] == NULL) {
      fprintf(stdout, "No multipole expansions at level %d\n", depth); 
      /*if(depth < 3) fprintf(stdout, " (Warning: no multipole acceleration)\n");
      else fprintf(stdout, "\n");*/
    }

#if ON == ON			/* ON == OFF produces the M2M error!! */
    /* NULL out pointers to same-geometry M2M mats for this level */
    for(i=0; i < 8; i++) multimats[i] = NULL;
#endif

  /* Hit nonempty cubes at this level assigning ptrs to precomputed   */
  /* M2M mats (for this lev), or if a kid is exact, computing Q2M matrices. */
    for(nextc=sys->multilist[depth]; nextc != NULL; nextc = nextc->mnext) {
      
#if DISSYN == ON
      multicnt[nextc->level]++;
#endif

    /* Save space for upvector sizes, upvect ptrs, and upmats. */
      nextc->multisize = numterms;
      if(numterms > 0) {
	CALLOC(nextc->multi, numterms, double, ON, AMSC);
      }
      if(nextc->upnumvects) {
	CALLOC(nextc->upnumeles, nextc->upnumvects, int, ON, AMSC);
	CALLOC(nextc->upvects, nextc->upnumvects, double*, ON, AMSC);
	CALLOC(nextc->upmats, nextc->upnumvects, double**, ON, AMSC);
      }
    /* Go through nonempty kids and fill in upvectors and upmats. */
      for(i=0, j=0; j < nextc->numkids; j++) {
	if((kid = nextc->kids[j]) != NULL) { /* NULL => empty kid cube */
	  if(kid->mul_exact == FALSE) { /* if kid has a multi */
	    nextc->upvects[i] = kid->multi;
	    nextc->upnumeles[i] = kid->multisize;
	    if(multimats[j] == NULL) { /* Build the needed matrix only once. */
	      multimats[j] = mulMulti2Multi(kid->x, kid->y, kid->z, nextc->x, 
					    nextc->y, nextc->z, order);
	    }
	    nextc->upmats[i] = multimats[j];

#if DMTCNT == ON
	    M2Mcnt[kid->level][nextc->level]++; /* cnts usage, ~computation */
#endif				/* only at most 8 mats really built/level */

	  }
	  else {		/* if kid is exact, has no multi */
	    nextc->upvects[i] = kid->upvects[0];
	    nextc->upnumeles[i] = kid->upnumeles[0];
	    nextc->upmats[i] = mulQ2Multi(kid->chgs, kid->nbr_is_dummy[0],
					  kid->upnumeles[0],
					  nextc->x, nextc->y, nextc->z, order);

#if DMTCNT == ON
	    Q2Mcnt[kid->level][nextc->level]++;
#endif

	  }
	  i++;			/* only increments if kid is not empty */
	}
      }
    }
  }
}

/*
  builds the transformation matrices for the final evaluation pass (M2P, L2P)
  for all the direct list (generated by linkcubes(), non-empty
  lowest level cubes) cubes:

  for each cube A in the direct list:
  1) if the cube is not exact (always the case if ADAPT = OFF)
     a) and if DNTYPE = GRENGD build an L2P matrix from A to A 
     b) and if DNTYPE = NOSHFT build an L2P matrix from each of A's ancestors
        with level > 1 (including A) to A
     c) and if DNTYPE = NOLOCL build an M2P matrix from each of A's fake 
        ilist entries to A (same action as 2b)
  2) if the cube is exact, find the 1st ancestor of A, cube B, 
     which either is not exact and is at level 2,3,4... or is at level 1
     a) if B is at level 2,3,4... 
        i) if DNTYPE = GRENGD, construct an L2P from B to A and M2P's
           from the cubes in the true interaction lists of A and all its
	   ancestors up to and including B (a partial fake interaction list)
	j) if DNTYPE = NOSHFT, find cube C, the ancestor of B at level 1;
	   construct L2P's from the ancestors of B (including B but not C)
	   to A and Q- or M2P's from the cubes in the true interaction lists 
	   of A and all its ancestors up to and including B (a partial fake 
	   interaction list)
	k) if DNTYPE = NOLOCL, do 2b
     b) if B is at level 1 construct M2P's for all the cubes in A's
        fake interaction list

  true interaction list - RADINTER = OFF, those sibling
  (same level) cubes of a given cube who are children of the neighbors
  of the given cube's parent and are not neighbors of the given cube 
  - ie those cubes required to cover charges well separated from the given
  cube but not accounted for in the parent's local expansion 
  - the flag NNBRS is the number of sibling cube "shells" taken as neighbors 
  
  fake interaction list - RADINTER = OFF, the combined true interaction lists
  of a given cube and all its ancestors at levels 2,3,4...

  if RADINTER = ON, any 8 siblings of the given cube which form a well 
  separated cube one level up are included in the lists as a single higher
  level cube
  
  if ADAPT = OFF, no cube is exact so step 2 is never done

  this routine is used alone if compiled with DNTYPE = NOLOCL or after
  mulMatDown, which produces M2L and L2L matrices (DNTYPE = GRENGD) or
  just M2L matrices (DNTYPE = NOSHFT) --  DNTYPE = GRENGD does the full
  Greengard hiearchical downward pass

*/
void mulMatEval(sys)
ssystem *sys;
{
  int i, j, k, ttlvects, vects;
  cube *na, *nc, *nexti;

  if(sys->depth < 2) return;	/* ret if upward pass not possible/worth it */

  for(nc = sys->directlist; nc != NULL; nc = nc->dnext) {

    ASSERT(nc->level == sys->depth);
    ASSERT(nc->upnumvects > 0);

    /* allocate space for evaluation pass vectors; check nc's ancestors */
    /* First count the number of transformations to do. */
    for(na = nc, ttlvects = 0; na->level > 1; na = na->parent) { 
      if(na->loc_exact == FALSE && DNTYPE != NOLOCL) {
	ttlvects++;  /* allow for na to na local expansion (L2P) */
	if(DNTYPE == GRENGD) break; /* Only one local expansion if shifting. */
      }
      else {
	ttlvects += na->interSize; /* room for Q2P and M2P xformations */
      }
    }
    nc->evalnumvects = ttlvects; /* save ttl # of transformations to do */
    if(ttlvects > 0) {
      CALLOC(nc->evalvects, ttlvects, double*, ON, AMSC);
      CALLOC(nc->evalnumeles, ttlvects, int, ON, AMSC);
      CALLOC(nc->evalmats, ttlvects, double**, ON, AMSC);
    }
    
#if DILIST == ON
    fprintf(stdout, "\nInteraction list (%d entries) for ", ttlvects);
    disExParsimpcube(nc);
#endif
    
    /* set up exp/charge vectors and L2P, Q2P and/or M2P matrices as req'd */
    for(j=0, na = nc, ttlvects = 0; na->level > 1; na = na->parent) { 
      if(na->loc_exact == FALSE && DNTYPE != NOLOCL) {  
	/* build matrices for local expansion evaluation */
	nc->evalmats[j] = mulLocal2P(na->x, na->y, na->z, nc->chgs,
				     nc->upnumeles[0], sys->order);
	nc->evalnumeles[j] = na->localsize;
	nc->evalvects[j] = na->local;
	j++; 
	
#if DMTCNT == ON
	L2Pcnt[na->level][nc->level]++;
#endif
	
#if DILIST == ON
	fprintf(stdout, "L2P: ");
	disExtrasimpcube(na);
#endif
	if(DNTYPE == GRENGD) break; /* Only one local expansion if shifting. */
      }
      else { /* build matrices for ancestor's (or cube's if 1st time) ilist */
	for(i=0; i < na->interSize; i++) {
	  nexti = na->interList[i];
	  if(nexti->mul_exact == TRUE) {
	    nc->evalvects[j] = nexti->upvects[0];
	    nc->evalmats[j] = Q2P(nexti->chgs, nexti->upnumeles[0], 
				  nexti->nbr_is_dummy[0], nc->chgs, 
				  nc->upnumeles[0], TRUE);
	    nc->evalnumeles[j] = nexti->upnumeles[0];
	    j++;

#if DMTCNT == ON
	    Q2Pcnt[nexti->level][nc->level]++;
#endif

#if DILIST == ON
	    fprintf(stdout, "Q2P: ");
	    disExtrasimpcube(nexti);
#endif
	  }
	  else {
	    nc->evalvects[j] = nexti->multi;
	    nc->evalmats[j] = mulMulti2P(nexti->x, nexti->y, nexti->z, 
					 nc->chgs, nc->upnumeles[0], 
					 sys->order);
	    nc->evalnumeles[j] = nexti->multisize;
	    j++;
	    
#if DMTCNT == ON
	    M2Pcnt[nexti->level][nc->level]++;
#endif

#if DILIST == ON
	    fprintf(stdout, "M2P: ");
	    disExtrasimpcube(nexti);
#endif
	  }
	}
      }
    }
  }
}


/* 
  sets up matrices for the downward pass
  For each cube in local list (parents always in list before kids):
  1) parent's local to child's local unless DNTYPE = NOSHFT or no parent local
  2) multipoles for (Parent+parent's nbrs - child nbrs) to child's local
  -eval is sum of ancestral local evals for each lowest lev cube if NOSHFT
    otherwise only lowest level local is evaluated (see mulMatEval)
  -with ADAPT = OFF no cube is exact so local list is all non-empty cube lev>1
  -mats that give potentials (M2P, L2P, Q2P) are calculated in mulMatEval()
  -this routine makes only L2L, M2L and Q2L matrices
*/
mulMatDown(sys)
ssystem *sys;
{
  int i, j, vects;
  cube *nc, *parent, *ni;
  int depth;

  ASSERT(DNTYPE != NOLOCL);	/* use mulMatEval() alone if NOLOCL */

  for(depth = 2; depth <= sys->depth; depth++) { /* no locals before level 2 */
    for(nc=sys->locallist[depth]; nc != NULL; nc = nc->lnext) {

      /* Allocate for interaction list, include one for parent if needed */
      if((depth <= 2) || (DNTYPE == NOSHFT)) vects = nc->interSize;
      else vects = nc->interSize + 1;
      nc->downnumvects = vects;
      if(vects > 0) {
	CALLOC(nc->downvects, vects, double*, ON, AMSC);
	CALLOC(nc->downnumeles, vects, int, ON, AMSC);
	CALLOC(nc->downmats, vects, double**, ON, AMSC);
      }

      parent = nc->parent;
      ASSERT(parent->loc_exact == FALSE); /* has >= #evals of any of its kids*/

#if DISSYN == ON
      localcnt[nc->level]++;
#endif

      if((depth <= 2) || (DNTYPE == NOSHFT)) i = 0; /* No parent local. */
      else { /* Create the mapping matrix for the parent to kid. */
	i = 1;

	nc->downmats[0] = mulLocal2Local(parent->x, parent->y, parent->z,
					 nc->x, nc->y, nc->z, sys->order);
	nc->downnumeles[0] = parent->localsize;
	nc->downvects[0] = parent->local;

#if DMTCNT == ON
	L2Lcnt[parent->level][nc->level]++;
#endif
      }

      /* Go through the interaction list and create mapping matrices. */
      for(j = 0; j < nc->interSize; j++, i++) {
	ni = nc->interList[j];
	if(ni->mul_exact == TRUE) {	/* ex->ex (Q2P) xforms in mulMatEval */
	  nc->downvects[i] = ni->upvects[0];
	  nc->downmats[i] = mulQ2Local(ni->chgs, ni->upnumeles[0],
				       ni->nbr_is_dummy[0],
				       nc->x, nc->y, nc->z, sys->order);
	  nc->downnumeles[i] = ni->upnumeles[0];
#if DMTCNT == ON
	  Q2Lcnt[ni->level][nc->level]++;
#endif
	}
	else {
	  nc->downvects[i] = ni->multi;
	  nc->downmats[i] = mulMulti2Local(ni->x, ni->y, ni->z, nc->x, 
					   nc->y, nc->z, sys->order);
	  nc->downnumeles[i] = ni->multisize;
#if DMTCNT == ON
	  M2Lcnt[ni->level][nc->level]++;
#endif
	}
      }
    }
  }
}