mulmats_old.c

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

/*
Copyright (c) 1990 Massachusetts Institute of Technology, Cambridge, MA.
All rights reserved.

This Agreement gives you, the LICENSEE, certain rights and obligations.
By using the software, you indicate that you have read, understood, and
will comply with the terms.

Permission to use, copy and modify for internal, noncommercial purposes
is hereby granted.  Any distribution of this program or any part thereof
is strictly prohibited without prior written consent of M.I.T.

Title to copyright to this software and to any associated documentation
shall at all times remain with M.I.T. and LICENSEE agrees to preserve
same.  LICENSEE agrees not to make any copies except for LICENSEE'S
internal noncommercial use, or to use separately any portion of this
software without prior written consent of M.I.T.  LICENSEE agrees to
place the appropriate copyright notice on any such copies.

Nothing in this Agreement shall be construed as conferring rights to use
in advertising, publicity or otherwise any trademark or the name of
"Massachusetts Institute of Technology" or "M.I.T."

M.I.T. MAKES NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED.  By
way of example, but not limitation, M.I.T. MAKES NO REPRESENTATIONS OR
WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE OR
THAT THE USE OF THE LICENSED SOFTWARE COMPONENTS OR DOCUMENTATION WILL
NOT INFRINGE ANY PATENTS, COPYRIGHTS, TRADEMARKS OR OTHER RIGHTS.
M.I.T. shall not be held liable for any liability nor for any direct,
indirect or consequential damages with respect to any claim by LICENSEE
or any third party on account of or arising from this Agreement or use
of this software.
*/

#include "mulGlobal.h"

double **Q2P(), **Q2PDiag();
double **mulMulti2P(), **mulQ2Multi(), **mulMulti2Multi();
double **mulLocal2Local(), **mulLocal2P(), **mulQ2Local(), **mulMulti2Local();

int *localcnt, *multicnt, *evalcnt;	/* counts of builds done by level */

int **Q2Mcnt, **Q2Lcnt, **Q2Pcnt, **L2Lcnt; /* counts of xformation mats */
int **M2Mcnt, **M2Lcnt, **M2Pcnt, **L2Pcnt, **Q2PDcnt;

/*
MulMatDirect creates the matrices for the piece of the problem that is done
directly exactly.
*/
mulMatDirect(sys)
ssystem *sys;
{
  cube *nextc, *nextnbr;
  int i, nummats, **temp;
  extern double lutime, dirtime;

#if DIRSOL == ON || EXPGCR == ON
  double **ludecomp();
  extern double *trimat, *sqrmat; /* flattened triangular, square matrices */
  extern int up_size, eval_size;
  extern int *real_index;	/* for map btwn condensed/expanded vectors */
#endif

  /* First count the number of matrices to be done directly. */
  for(nextc=sys->directlist; nextc != NULL; nextc = nextc->dnext) {
    for(nummats=1, i=0; i < nextc->numnbrs; i++) {
      nextnbr = nextc->nbrs[i];
      ASSERT(nextnbr->upnumvects > 0);
      nummats++;
    }

  /* Allocate space for the vects and mats. */
    nextc->directnumvects = nummats;
    CALLOC(nextc->directq, nummats, double*, ON, AMSC);
    CALLOC(temp, nummats, int*, ON, AMSC);
    CALLOC(nextc->directnumeles, nummats, int, ON, AMSC);
    CALLOC(nextc->directmats, nummats, double**, ON, AMSC);
    CALLOC(nextc->precondmats, nummats, double**, ON, AMSC);

    /* initialize the pointer from this cube to its part of dummy vector
       - save the self part found in indexkid() */
    temp[0] = nextc->nbr_is_dummy[0];
    nextc->nbr_is_dummy = temp;
  }

/* Now place in the matrices. */
  for(nextc=sys->directlist; nextc != NULL; nextc = nextc->dnext) {
    nextc->directq[0] = nextc->upvects[0];
    nextc->directnumeles[0] = nextc->upnumeles[0];

    starttimer;
#if DIRSOL == ON || EXPGCR == ON
    if(nextc == sys->directlist) {
      if(eval_size < MAXSIZ) {
	fprintf(stderr, 
		"mulMatDirect: non-block direct methods not supported\n");
	exit(1);
	/* if this is going to work, need a special, condensing Q2P
	   as well as some way to use it in the framework of the GCR loop */
	nextc->directmats[0] = Q2P(nextc->chgs, eval_size,
				   nextc->nbr_is_dummy[0], nextc->chgs, 
				   eval_size, TRUE);
      }
      else blkQ2Pfull(sys->directlist, up_size, eval_size,
		      &trimat, &sqrmat, &real_index, sys->is_dummy);
    }
    else nextc->directmats[0] 
	= Q2PDiag(nextc->chgs, nextc->upnumeles[0], nextc->nbr_is_dummy[0],
		  TRUE);
#else
    nextc->directmats[0] 
	= Q2PDiag(nextc->chgs, nextc->upnumeles[0], nextc->nbr_is_dummy[0],
		  TRUE);
    nextc->precondmats[0] 
	= Q2PDiag(nextc->chgs, nextc->upnumeles[0], nextc->nbr_is_dummy[0],
		  FALSE);
    /*dumpMatCor(nextc->directmats[0], (double *)NULL, nextc->upnumeles[0]);*/

#endif
    stoptimer;
    dirtime += dtime;

#if DSQ2PD == ON
    dumpQ2PDiag(nextc);
#endif

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

#if DIRSOL == ON
    /* transform A into LU */
    if(eval_size > MAXSIZ) {
      blkLUdecomp(sqrmat, trimat, up_size);
    }
    else if(nextc == sys->directlist) {
      starttimer;
      nextc->directlu = ludecomp(nextc->directmats[0], eval_size, TRUE);
      stoptimer;
      lutime += dtime;
    }
#endif
    
    starttimer;
    for(nummats=1, i=0; i < nextc->numnbrs; i++) {
      nextnbr = nextc->nbrs[i];
      ASSERT(nextnbr->upnumvects > 0);
      nextc->directq[nummats] = nextnbr->upvects[0];
      nextc->nbr_is_dummy[nummats] = nextnbr->nbr_is_dummy[0];
      nextc->directnumeles[nummats] = nextnbr->upnumeles[0];
      nextc->directmats[nummats] = Q2P(nextnbr->chgs, 
				       nextnbr->upnumeles[0], 
				       nextnbr->nbr_is_dummy[0],
				       nextc->chgs, nextc->upnumeles[0],
				       TRUE);
      nextc->precondmats[nummats++] = Q2P(nextnbr->chgs, 
					  nextnbr->upnumeles[0], 
					  nextnbr->nbr_is_dummy[0],
					  nextc->chgs, nextc->upnumeles[0],
					  FALSE);
#if DMTCNT == ON
      Q2Pcnt[nextc->level][nextnbr->level]++;
#endif
    }
    stoptimer;
    dirtime += dtime;
  }
}


/*
MulMatPrecond creates the preconditioner matrix
*/
bdmulMatPrecond(sys)
ssystem *sys;
{
  cube *nc, *kid, *kidnbr;
  double **mat, **nbrmat;
  int i, j, k, l, kidi;
  int kidsize, nbrsize, size, row, col, first, offset;
  double **ludecomp(), factor;
  charge *pc;
  surface *surf;

  for(nc=sys->precondlist; nc != NULL; nc = nc->pnext) {

    /* find total number of charges in cube to dimension P. */
    for(size=0, i=0; i < nc->numkids; i++) {
      kid = nc->kids[i];
      if(kid != NULL) {
	ASSERT(kid->level == sys->depth);
	size += kid->directnumeles[0];  /* Equals number of charges. */
      }
    }

    /* allocate and zero a preconditioner matrix. */
    MALLOC(mat, size, double*, ON, AMSC);
    for(i = 0; i < size; i++) {
      MALLOC(mat[i], size, double, ON, AMSC);
    }
    for(i = 0; i < size; i++) {
      for(j = 0; j < size; j++) {
	mat[i][j] = 0.0;
      }
    }

    /* Chase through the kids to place in potential coeffs. */
    for(first = TRUE, row=0, kidi=0; kidi < nc->numkids; kidi++) {
      kid = nc->kids[kidi];
      if(kid != NULL) {
	/* Exploit the hierarchical charge numbering to get precond vector. */
	if(first == TRUE) {
	  first = FALSE;
	  nc->prevectq = kid->directq[0];
	  nc->prevectp = kid->eval;
	}
	/* Get the diagonal block of P^{-1}. */
	kidsize = kid->directnumeles[0];
	for(k = kidsize - 1; k >= 0; k--) {
	  for(l = kidsize - 1; l >= 0; l--) {
	    mat[row + k][row + l] = kid->directmats[0][k][l];
	  }
	}
	/* Get the off-diagonals of P^{-1}. */
	for(col = 0, i = 0; i < nc->numkids; i++) {
	  kidnbr = nc->kids[i];
	  if(kidnbr != NULL) {
	    if(kidnbr != kid) {
	      /* Chase thru list of nbrs to get matrix associated with this 
		 kidnbr.  Note, this is because the kid list and nbr matrix
		 list are in different orders, could be fixed. */
	      for(j = kid->numnbrs-1; j >= 0; j--) {
		if(kidnbr == kid->nbrs[j]) {
		  nbrmat = kid->directmats[j+1];
		  nbrsize = kidnbr->directnumeles[0];
		  for(k = kidsize - 1; k >= 0; k--) {
		    for(l = nbrsize - 1; l >= 0; l--) {
		      mat[row + k][col + l] = nbrmat[k][l];
		    }
		  }
		  break;
		}
	      }
	    }
	    col += kidnbr->directnumeles[0];
	  }
	}
	ASSERT(col == size);
	row += kidsize;
      }
    }    
    ASSERT(row == size);

    nc->precond = ludecomp(mat, size, FALSE);
    nc->presize = size;
  }
}
      

/* This near picks up only the hamming distance one cubes. */    
#define HNEAR(nbr, nj, nk, nl) \
((ABS((nbr)->j - (nj)) + ABS((nbr)->k - (nk)) + ABS((nbr)->l - (nl))) <= 1)

/* This near picks up all 27 neighboring cubes. */
#define NEAR(nbr, nj, nk, nl) \
((ABS((nbr)->j - (nj)) <= 1) && \
 (ABS((nbr)->k - (nk)) <= 1) && \
 (ABS((nbr)->l - (nl)) <= 1))

/* This near picks only the diagonal, for testing. */
#define DNEAR(nbr, nj, nk, nl) \
(((nbr)->j == (nj)) && \
 ((nbr)->k == (nk)) && \
 ((nbr)->l == (nl)) )
	
olmulMatPrecond(sys)
ssystem *sys;
{
  cube *nc, *nnbr, *nnnbr;
  double **mat, **nmat;
  int i, j, k, l, m;
  int maxsize, nsize, nnsize, nnnsize, *reorder;
  int nj, nk, nl, offset, noffset;
  int dindex, *nc_dummy, *nnbr_dummy, *nnnbr_dummy;
  static int *is_dummy;		/* local dummy flag vector, stays around */
  static int big_mat_size = 0;	/* size of previous mat */
  charge **nnnbr_pc, **nnbr_pc, **nc_pc, **mpc, *dp;
  surface *surf;
  double factor;

/* Figure out the max number of elements in any set of near cubes. */
  for(maxsize=0, nc=sys->directlist; nc != NULL; nc = nc->dnext) {
    nsize = nc->directnumeles[0];
    nj = nc->j;
    nk = nc->k;
    nl = nc->l;
    for(i=0; i < nc->numnbrs; i++) {
      nnbr = nc->nbrs[i];
      if(NEAR(nnbr, nj, nk, nl)) nsize += nnbr->directnumeles[0];
    }
    maxsize = MAX(nsize, maxsize);
  }

/* Allocate a matrix big enough for any set of 7. */
  printf("max direct size =%d\n", maxsize);
  MALLOC(reorder, maxsize,int, ON, AMSC);
  MALLOC(mat, maxsize, double*, ON, AMSC);
  for(i=0; i < maxsize; i++) {
    MALLOC(mat[i], maxsize, double, ON, AMSC);
  }

/* Now go fill-in a matrix. */
  for(maxsize=0, nc=sys->directlist; nc != NULL; nc = nc->dnext) {
    nsize = nc->directnumeles[0];
    nc_dummy = nc->nbr_is_dummy[0];
    nc_pc = nc->chgs;
#if CHKDUM == ON
    chkDummyList(nc_pc, nc_dummy, nsize);
#endif
    nj = nc->j;
    nk = nc->k;
    nl = nc->l;
    for(i = nsize - 1; i >= 0; i--) {
      if(nc_dummy[i]) continue;	/* dummy rows copied only in divived diff */
      if(nc_pc[i]->surf->type != DIELEC) {
	for(j = nsize - 1; j >= 0; j--) {
	  mat[i][j] = nc->directmats[0][i][j];
	}