mpb.c

麻省理工的计算光子晶体的程序

     if (!have_old_fields || reset_fields)
	  randomize_fields();

     {
	  int ierr = check_maxwell_dielectric(mdata, negative_epsilon_okp);
	  if (ierr == 1)
	       mpi_one_fprintf(stderr,
			   "ERROR: non positive-definite dielectric tensor\n");
	  else if (ierr == 2)
	       mpi_one_fprintf(stderr, 
		       "ERROR: dielectric tensor must not couple xy "
		       "plane with z direction for 2D TE/TM calculations\n");
	  CHECK(!ierr, "invalid dielectric function\n");
     }

     evectmatrix_flops = eigensolver_flops; /* reset, if changed */
}

/**************************************************************************/

/* When we are solving for a few bands at a time, we solve for the
   upper bands by "deflation"--by continually orthogonalizing them
   against the already-computed lower bands.  (This constraint
   commutes with the eigen-operator, of course, so all is well.) */

typedef struct {
     evectmatrix Y;  /* the vectors to orthogonalize against; Y must
			itself be normalized (Yt Y = 1) */
     int p;  /* the number of columns of Y to orthogonalize against */
     scalar *S;  /* a matrix for storing the dot products; should have
		    at least p * X.p elements (see below for X) */
     scalar *S2; /* a scratch matrix the same size as S */
} deflation_data;

static void deflation_constraint(evectmatrix X, void *data)
{
     deflation_data *d = (deflation_data *) data;

     CHECK(X.n == d->Y.n && d->Y.p >= d->p, "invalid dimensions");

     /* (Sigh...call the BLAS functions directly since we are not
	using all the columns of Y...evectmatrix is not set up for
	this case.) */

     /* compute S = Xt Y (i.e. all the dot products): */
     blasglue_gemm('C', 'N', X.p, d->p, X.n,
		   1.0, X.data, X.p, d->Y.data, d->Y.p, 0.0, d->S2, d->p);
     mpi_allreduce(d->S2, d->S, d->p * X.p * SCALAR_NUMVALS,
		   real, SCALAR_MPI_TYPE, MPI_SUM, MPI_COMM_WORLD);

     /* compute X = X - Y*St = (1 - Y Yt) X */
     blasglue_gemm('N', 'C', X.n, X.p, d->p,
		   -1.0, d->Y.data, d->Y.p, d->S, d->p,
		   1.0, X.data, X.p);
}

/**************************************************************************/

/* Solve for the bands at a given k point.
   Must only be called after init_params! */
void solve_kpoint(vector3 kvector)
{
     int i, total_iters = 0, ib, ib0;
     real *eigvals;
     real k[3];
     int flags;
     deflation_data deflation;
     int prev_parity;

     mpi_one_printf("solve_kpoint (%g,%g,%g):\n",
		    kvector.x, kvector.y, kvector.z);
     
     curfield_reset();

     if (num_bands == 0) {
	  mpi_one_printf("  num-bands is zero, not solving for any bands\n");
	  return;
     }

     if (!mdata) {
	  mpi_one_fprintf(stderr,
			  "init-params must be called before solve-kpoint!\n");
	  return;
     }

     /* if this is the first k point, print out a header line for
	for the frequency grep data: */
     if (!kpoint_index && mpi_is_master()) {
	  printf("%sfreqs:, k index, k1, k2, k3, kmag/2pi",
		 parity_string(mdata));
	  for (i = 0; i < num_bands; ++i)
	       printf(", %s%sband %d",
		      parity_string(mdata),
		      mdata->parity == NO_PARITY ? "" : " ",
		      i + 1);
	  printf("\n");
     }

     prev_parity = mdata->parity;
     cur_kvector = kvector;
     vector3_to_arr(k, kvector);
     update_maxwell_data_k(mdata, k, G[0], G[1], G[2]);
     CHECK(mdata->parity == prev_parity,
	   "k vector is incompatible with specified parity");

     CHK_MALLOC(eigvals, real, num_bands);

     flags = eigensolver_flags; /* ctl file input variable */
     if (verbose)
	  flags |= EIGS_VERBOSE;

     /* constant (zero frequency) bands at k=0 are handled specially,
        so remove them from the solutions for the eigensolver: */
     if (mdata->zero_k && !mtdata) {
	  int in, ip;
	  ib0 = maxwell_zero_k_num_const_bands(H, mdata);
	  for (in = 0; in < H.n; ++in)
	       for (ip = 0; ip < H.p - ib0; ++ip)
		    H.data[in * H.p + ip] = H.data[in * H.p + ip + ib0];
	  evectmatrix_resize(&H, H.p - ib0, 1);
     }
     else
	  ib0 = 0; /* solve for all bands */

     /* Set up deflation data: */
     if (H.data != Hblock.data) {
	  deflation.Y = H;
	  deflation.p = 0;
	  CHK_MALLOC(deflation.S, scalar, H.p * Hblock.p);
	  CHK_MALLOC(deflation.S2, scalar, H.p * Hblock.p);
     }

     for (ib = ib0; ib < num_bands; ib += Hblock.alloc_p) {
	  evectconstraint_chain *constraints;
	  int num_iters;

	  /* don't solve for too many bands if the block size doesn't divide
	     the number of bands: */
	  if (ib + mdata->num_bands > num_bands) {
	       maxwell_set_num_bands(mdata, num_bands - ib);
	       for (i = 0; i < nwork_alloc; ++i)
		    evectmatrix_resize(&W[i], num_bands - ib, 0);
	       evectmatrix_resize(&Hblock, num_bands - ib, 0);
	  }

	  mpi_one_printf("Solving for bands %d to %d...\n",
			 ib + 1, ib + Hblock.p);

	  constraints = NULL;
	  constraints = evect_add_constraint(constraints,
					     maxwell_parity_constraint,
					     (void *) mdata);

	  if (mdata->zero_k)
	       constraints = evect_add_constraint(constraints,
						  maxwell_zero_k_constraint,
						  (void *) mdata);

	  if (Hblock.data != H.data) {  /* initialize fields of block from H */
	       int in, ip;
	       for (in = 0; in < Hblock.n; ++in)
		    for (ip = 0; ip < Hblock.p; ++ip)
			 Hblock.data[in * Hblock.p + ip] =
			      H.data[in * H.p + ip + (ib-ib0)];
	       deflation.p = ib-ib0;
	       if (deflation.p > 0)
		    constraints = evect_add_constraint(constraints,
						       deflation_constraint,
						       &deflation);
	  }

	  if (mtdata) {  /* solving for bands near a target frequency */
	       if (eigensolver_davidsonp)
		    eigensolver_davidson(
			 Hblock, eigvals + ib,
			 maxwell_target_operator, (void *) mtdata,
			 simple_preconditionerp ? 
			 maxwell_target_preconditioner :
			 maxwell_target_preconditioner2,
			 (void *) mtdata,
			 evectconstraint_chain_func,
			 (void *) constraints,
			 W, nwork_alloc, tolerance, &num_iters, flags, 0.0);
	       else
		    eigensolver(Hblock, eigvals + ib,
				maxwell_target_operator, (void *) mtdata,
				simple_preconditionerp ? 
				maxwell_target_preconditioner :
				maxwell_target_preconditioner2,
				(void *) mtdata,
				evectconstraint_chain_func,
				(void *) constraints,
				W, nwork_alloc, tolerance, &num_iters, flags);

	       /* now, diagonalize the real Maxwell operator in the
		  solution subspace to get the true eigenvalues and
		  eigenvectors: */
	       CHECK(nwork_alloc >= 2, "not enough workspace");
	       eigensolver_get_eigenvals(Hblock, eigvals + ib,
					 maxwell_operator,mdata, W[0],W[1]);
	  }
	  else {
	       if (eigensolver_davidsonp)
		    eigensolver_davidson(
			 Hblock, eigvals + ib,
			 maxwell_operator, (void *) mdata,
			 simple_preconditionerp ?
			 maxwell_preconditioner :
			 maxwell_preconditioner2,
			 (void *) mdata,
			 evectconstraint_chain_func,
			 (void *) constraints,
			 W, nwork_alloc, tolerance, &num_iters, flags, 0.0);
	       else
		    eigensolver(Hblock, eigvals + ib,
				maxwell_operator, (void *) mdata,
				simple_preconditionerp ?
				maxwell_preconditioner :
				maxwell_preconditioner2,
				(void *) mdata,
				evectconstraint_chain_func,
				(void *) constraints,
				W, nwork_alloc, tolerance, &num_iters, flags);
	  }
	  
	  if (Hblock.data != H.data) {  /* save solutions of current block */
	       int in, ip;
	       for (in = 0; in < Hblock.n; ++in)
		    for (ip = 0; ip < Hblock.p; ++ip)
			 H.data[in * H.p + ip + (ib-ib0)] =
			      Hblock.data[in * Hblock.p + ip];
	  }

	  evect_destroy_constraints(constraints);
	  
	  mpi_one_printf("Finished solving for bands %d to %d after "
			 "%d iterations.\n", ib + 1, ib + Hblock.p, num_iters);
	  total_iters += num_iters * Hblock.p;
     }

     if (num_bands - ib0 > Hblock.alloc_p)
	  mpi_one_printf("Finished k-point with %g mean iterations/band.\n",
			 total_iters * 1.0 / num_bands);

     /* Manually put in constant (zero-frequency) solutions for k=0: */
     if (mdata->zero_k && !mtdata) {
	  int in, ip;
	  evectmatrix_resize(&H, H.alloc_p, 1);
	  for (in = 0; in < H.n; ++in)
	       for (ip = H.p - ib0 - 1; ip >= 0; --ip)
		    H.data[in * H.p + ip + ib0] = H.data[in * H.p + ip];
	  maxwell_zero_k_set_const_bands(H, mdata);
	  for (ib = 0; ib < ib0; ++ib)
	       eigvals[ib] = 0;
     }

     /* Reset scratch matrix sizes: */
     evectmatrix_resize(&Hblock, Hblock.alloc_p, 0);
     for (i = 0; i < nwork_alloc; ++i)
	  evectmatrix_resize(&W[i], W[i].alloc_p, 0);
     maxwell_set_num_bands(mdata, Hblock.alloc_p);

     /* Destroy deflation data: */
     if (H.data != Hblock.data) {
	  free(deflation.S2);
	  free(deflation.S);
     }

     if (num_write_output_vars > 1) {
	  /* clean up from prev. call */
	  free(freqs.items);
	  free(parity);
     }

     CHK_MALLOC(parity, char, strlen(parity_string(mdata)) + 1);
     parity = strcpy(parity, parity_string(mdata));

     iterations = total_iters; /* iterations output variable */

     /* create freqs array for storing frequencies in a Guile list */
     freqs.num_items = num_bands;
     CHK_MALLOC(freqs.items, number, freqs.num_items);
     
     set_kpoint_index(kpoint_index + 1);

     mpi_one_printf("%sfreqs:, %d, %g, %g, %g, %g",
		    parity,
		    kpoint_index, k[0], k[1], k[2],
		    vector3_norm(matrix3x3_vector3_mult(Gm, kvector)));
     for (i = 0; i < num_bands; ++i) {
	  freqs.items[i] =
	       negative_epsilon_okp ? eigvals[i] : sqrt(eigvals[i]);
	  mpi_one_printf(", %g", freqs.items[i]);
     }
     mpi_one_printf("\n");

     eigensolver_flops = evectmatrix_flops;

     free(eigvals);
}

/**************************************************************************/

/* Return a list of the z/y parities, one for each band. */

number_list compute_zparities(void)
{
     number_list z_parity;
     z_parity.num_items = num_bands;
     z_parity.items = maxwell_zparity(H, mdata);
     return z_parity;
}

number_list compute_yparities(void)
{
     number_list y_parity;
     y_parity.num_items = num_bands;
     y_parity.items = maxwell_yparity(H, mdata);
     return y_parity;
}

/**************************************************************************/

/* Compute the group velocity dw/dk in the given direction d (where
   the length of d is ignored).  d is in the reciprocal lattice basis.
   Should only be called after solve_kpoint.  Returns a list of the
   group velocities, one for each band, in units of c. */
number_list compute_group_velocity_component(vector3 d)
{
     number_list group_v;
     number *gv_scratch;
     real u[3];
     int i, ib;

     group_v.num_items = 0;  group_v.items = (number *) NULL;

     curfield_reset(); /* has the side effect of overwriting curfield scratch */

     if (!mdata) {
	  mpi_one_fprintf(stderr, "init-params must be called first!\n");
	  return group_v;
     }
     if (!kpoint_index) {
	  mpi_one_fprintf(stderr, "solve-kpoint must be called first!\n");
	  return group_v;
     }

     /* convert d to unit vector in Cartesian coords: */
     d = unit_vector3(matrix3x3_vector3_mult(Gm, d));
     u[0] = d.x; u[1] = d.y; u[2] = d.z;

     group_v.num_items = num_bands;
     CHK_MALLOC(group_v.items, number, group_v.num_items);
     CHK_MALLOC(gv_scratch, number, group_v.num_items);
     
     /* now, compute group_v.items = diag Re <H| curl 1/eps i u x |H>: */

     /* ...we have to do this in blocks of eigensolver_block_size since
	the work matrix W[0] may not have enough space to do it all at once. */
     
     for (ib = 0; ib < num_bands; ib += Hblock.alloc_p) {
	  if (ib + mdata->num_bands > num_bands) {
	       maxwell_set_num_bands(mdata, num_bands - ib);
	       evectmatrix_resize(&W[0], num_bands - ib, 0);
	       evectmatrix_resize(&Hblock, num_bands - ib, 0);
	  }
	  if (Hblock.data != H.data) {  /* initialize fields of block from H */
	       int in, ip;
	       for (in = 0; in < Hblock.n; ++in)
		    for (ip = 0; ip < Hblock.p; ++ip)
			 Hblock.data[in * Hblock.p + ip] =
			      H.data[in * H.p + ip + ib];
	  }

	  maxwell_ucross_op(Hblock, W[0], mdata, u);
	  evectmatrix_XtY_diag_real(Hblock, W[0],
				    group_v.items + ib, gv_scratch);
     }

     free(gv_scratch);

     /* Reset scratch matrix sizes: */
     evectmatrix_resize(&Hblock, Hblock.alloc_p, 0);
     evectmatrix_resize(&W[0], W[0].alloc_p, 0);
     maxwell_set_num_bands(mdata, Hblock.alloc_p);

     /* The group velocity is given by:

	grad_k(omega)*d = grad_k(omega^2)*d / 2*omega
	   = grad_k(<H|maxwell_op|H>)*d / 2*omega
	   = Re <H| curl 1/eps i u x |H> / omega
        
        Note that our k is in units of 2*Pi/a, and omega is in
        units of 2*Pi*c/a, so the result will be in units of c. */
     for (i = 0; i < num_bands; ++i) {
	  if (freqs.items[i] == 0)  /* v is undefined in this case */
	       group_v.items[i] = 0.0;  /* just set to zero */
	  else
	       group_v.items[i] /= 
		    negative_epsilon_okp ? sqrt(fabs(freqs.items[i]))
		    : freqs.items[i];
     }
     
     return group_v;
}

/**************************************************************************/