capsolve.c

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

#endif
}



/* 
Preconditioned(possibly) Generalized Conjugate Residuals.
*/
int gcr(sys, q, p, r, ap, bp, bap, size, maxiter, tol, chglist)
ssystem *sys;
double *q, *p, *ap, *r, tol;
double **bp, **bap;
int size, maxiter;
charge *chglist;
{
  int iter, i, j;
  double norm, beta, alpha, maxnorm;
  extern double conjtime;
#if EXPGCR == ON
  extern double *sqrmat;
#endif

  /* NOTE ON EFFICIENCY: all the loops of length "size" could have */
  /*   if(sys->is_dummy[i]) continue; as their first line to save some ops */
  /* currently the entries corresponding to dummy panels are set to zero */

  for(iter = 0; iter < maxiter; iter++) {

    /* allocate the back vectors if they haven't been already (22OCT90) */
    if(bp[iter] == NULL) {
      CALLOC(bp[iter], size+1, double, ON, AMSC);
      CALLOC(bap[iter], size+1, double, ON, AMSC);
    }

    for(i=1; i <= size; i++) {
      bp[iter][i] = p[i] = r[i];
    }

    computePsi(sys, p, ap, size, chglist);
    
    starttimer;
    for(i=1; i <= size; i++) {
      bap[iter][i] = ap[i];
    }
    
    /* Subtract the backward projections from p and ap. */
    for(j= 0; j < iter; j++) {
      INNER(beta, ap, bap[j], size);
      for(i=1; i <= size; i++) {
	bp[iter][i] -= beta * bp[j][i];
	bap[iter][i] -= beta * bap[j][i];
      }
    }
    
    /* Normalize the p and ap vectors so that ap*ap = 1. */
    INNER(norm, bap[iter], bap[iter], size);
    norm = sqrt(norm);
    for(i=1; i <= size; i++) {
      bap[iter][i] /= norm;
      bp[iter][i] /= norm;
    }
    
    /* Compute the projection in the p direction and get the next p. */
    INNER(alpha, r, bap[iter], size);
    for(i=1; i <= size; i++) {
      q[i] += alpha * bp[iter][i];
      r[i] -= alpha * bap[iter][i];
    }

    /* Check convergence. */
    for(maxnorm = 0.0, i=1; i <= size; i++) maxnorm = MAX(ABS(r[i]),maxnorm);
#if ITRDAT == ON
    INNER(norm, r, r, size);
    fprintf(stdout, "max res = %g ||res|| = %g\n", maxnorm, sqrt(norm));
#else
    fprintf(stdout, "%d ", iter+1);
    if((iter+1) % 15 == 0) fprintf(stdout, "\n");
#endif
    fflush(stdout);
    stoptimer;
    conjtime += dtime;
    if(maxnorm < tol) break;
  }
/*  
#if ITRDAT == ON
  printf("Total iterative loop iterations = %d\n", iter); 
#endif 
*/
  
#if PRECOND != NONE
  /* Undo the preconditioning to get the real q. */
  for(i=1; i <= size; i++) {
    p[i] = q[i];
    ap[i] = 0.0;
  }
  mulPrecond(sys, PRECOND);
  for(i=1; i <= size; i++) {
    q[i] = p[i];
  }
#endif
  
  if(maxnorm >= tol) {
    fprintf(stdout, "\ngcr: WARNING exiting without converging\n");
  }
  return(iter+1);
}


/* 
  Preconditioned(possibly) Generalized Minimum Residual. 
  */
int gmres(sys, q, p, r, ap, bv, bh, size, maxiter, tol, chglist)
ssystem *sys;
double *q, *p, *ap, *r, tol;
double **bv, **bh;
int size, maxiter;
charge *chglist;
{
  int iter, i, j;
  double rnorm, norm, maxnorm=10.0;
  double beta, hi, hip1, length;
  extern double conjtime, prectime;
#if EXPGCR == ON
  extern double *sqrmat;
#endif

  static double *c=NULL, *s=NULL, *g=NULL, *y=NULL;
  
  starttimer;

  /* Allocation first time through. */
  if(c == NULL) {
    CALLOC(c, size+1, double, ON, AMSC);
    CALLOC(s, size+1, double, ON, AMSC);
    CALLOC(g, size+1, double, ON, AMSC);
    CALLOC(y, size+1, double, ON, AMSC);
  }
  
  /* Set up v^1 and g^0. */
  INNER(rnorm, r, r, size);
  rnorm = sqrt(rnorm);
  for(i=1; i <= size; i++) {
    p[i] = r[i] / rnorm;
    g[i] = 0.0;
  }
  g[1] = rnorm;

  stoptimer;
  conjtime += dtime;

#if ITRDAT == ON
  fprintf(stdout, "||res|| = %g\n", rnorm); /* initial guess residual norm */
#endif
  
  for(iter = 1; (iter <= maxiter) && (rnorm > tol); iter++) {
    
    starttimer;
    /* allocate the back vectors if they haven't been already */
    if(bv[iter] == NULL) {
      CALLOC(bv[iter], size+1, double, ON, AMSC);
      CALLOC(bh[iter], iter+2, double, ON, AMSC);
    }
    
    /* Save p as the v{iter}. */
    for(i=1; i <= size; i++) bv[iter][i] = p[i];
    
    stoptimer;
    conjtime += dtime;

    /* Form Av{iter}. */
    computePsi(sys, p, ap, size, chglist);

    starttimer;
    
    /* Initialize v^{iter+1} to Av^{iter}. */
    for(i=1; i <= size; i++) p[i] = ap[i];
    
    /* Make v^{iter+1} orthogonal to v^{i}, i <= iter. */
    for(j=1; j <= iter; j++) {
      INNER(hi, ap, bv[j], size);
      for(i=1; i <= size; i++) p[i] -= hi * bv[j][i];
      bh[iter][j] = hi;
    }
    
    /* Normalize v^{iter+1}. */
    INNER(norm, p, p, size);    
    norm = sqrt(norm);
    for(i=1; i <= size; i++) p[i] /= norm;
    bh[iter][iter+1] = norm;
    
    /* Apply rotations to new h column. */
    for(i=1; i < iter; i++) {
      hi = bh[iter][i];
      hip1 = bh[iter][i+1];
      bh[iter][i] = c[i] * hi - s[i] * hip1;
      bh[iter][i+1] = c[i] * hip1 + s[i] * hi;
    }
    
    /* Compute new rotations. */
    hi = bh[iter][iter];
    hip1 = bh[iter][iter+1];
    length = sqrt(hi * hi + hip1 * hip1);
    c[iter] = hi/length;
    s[iter] = -hip1/length;
    
    /* Apply new rotations. */
    bh[iter][iter] = c[iter] * hi - s[iter] * hip1;
    bh[iter][iter+1] = c[iter] * hip1 + s[iter] * hi;
    /* ASSERT(g[iter+1] == 0); WHY IS THIS HERE ??? */
    hi = g[iter];
    g[iter] = c[iter] * hi;
    g[iter+1] = s[iter] * hi;    
    
    rnorm = ABS(g[iter+1]);

    stoptimer;
    conjtime += dtime;
    
#if ITRDAT == ON
    fprintf(stdout, "||res|| = %g\n", rnorm);
#else
    fprintf(stdout, "%d ", iter);
    if((iter) % 15 == 0 && iter != 0) fprintf(stdout, "\n");
#endif
    fflush(stdout);
  }
  /* Decrement from the last increment. */
  iter--;

/*  
#if ITRDAT == ON
  printf("Total iterative loop iterations = %d\n", iter); 
#endif 
*/
  
  starttimer;
  
  /* Compute solution, note, bh is bh[col][row]. */
  for(i=1; i <= iter; i++) y[i] = g[i];
  for(i = iter; i > 0; i--) {
    y[i] /=  bh[i][i];
    for(j = i-1; j > 0; j--) {
      y[j] -= bh[i][j]*y[i];
    }
  }
  for(i=1; i <= size; i++) {
    q[i] = 0.0;
    for(j=1; j <= iter; j++) {
      q[i] += y[j] * bv[j][i];
    }
  }

  stoptimer;
  conjtime += dtime;
  
#if PRECOND != NONE
  /* Undo the preconditioning to get the real q. */
  starttimer;
  for(i=1; i <= size; i++) {
    p[i] = q[i];
    ap[i] = 0.0;
  }
  mulPrecond(sys, PRECOND);
  for(i=1; i <= size; i++) {
    q[i] = p[i];
  }
  stoptimer;
  prectime += dtime;
#endif

  if(rnorm > tol) {
    fprintf(stdout, "\ngmres: WARNING exiting without converging\n");
  }
  return(iter);
}

/* 
ComputePsi computes the potential from the charge vector, or may
include a preconditioner.  It is assumed that the vectors for the
charge and potential have already been set up and that the potential
vector has been zeroed.  ARBITRARY VECTORS CAN NOT BE USED.
*/

computePsi(sys, q, p, size, chglist)
ssystem *sys;
double *q, *p;
int size;
charge *chglist;
{
  extern double dirtime, uptime, downtime, evaltime, prectime;
  extern int real_size;
  int i;

  ASSERT(p == sys->p);
  ASSERT(q == sys->q);

  for(i=1; i <= size; i++) p[i] = 0;

#if PRECOND != NONE
  starttimer;
  mulPrecond(sys, PRECOND);
  stoptimer;
  prectime += dtime;
#endif

#if EXPGCR == ON
  blkCompressVector(q+1, size, real_size, sys->is_dummy+1);
  blkAqprod(p+1, q+1, real_size, sqrmat);	/* offset since index from 1 */
  blkExpandVector(p+1, size, real_size); /* ap changed to p, r chged to q */
  blkExpandVector(q+1, size, real_size); /*    7 Oct 91 */
#else
  starttimer;
  mulDirect(sys);
  stoptimer;
  dirtime += dtime;

  starttimer;
  mulUp(sys);
  stoptimer;
  uptime += dtime;

#if DUPVEC == ON
  dumpLevOneUpVecs(sys);
#endif

#if DNTYPE == NOSHFT
  mulDown(sys);		/* do downward pass without local exp shifts */
#endif

#if DNTYPE == GRENGD
  mulDown(sys);	       	/* do heirarchical local shift dwnwd pass */
#endif
  stoptimer;
  downtime += dtime;

  starttimer;
#if MULTI == ON
  mulEval(sys);		/* evaluate either locals or multis or both */
#endif
  stoptimer;
  evaltime += dtime;

#if DMPCHG == LAST
  fprintf(stdout, "\nPanel potentials divided by areas\n");
  dumpChgDen(stdout, p, chglist);
  fprintf(stdout, "End panel potentials\n");
#endif

  /* convert the voltage vec entries on dielectric i/f's into eps1E1-eps2E2 */
  compute_electric_fields(sys, chglist);

#if OPCNT == ON
  printops();
  exit(0);
#endif				/* OPCNT == ON */
#endif				/* EXPGCR == ON */
}