mpb.c

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

/* Copyright (C) 1999, 2000, 2001, 2002, Massachusetts Institute of Technology.
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */

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

/* Here, we define the external functions callable from Guile, as defined
   by mpb.scm. */

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

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <math.h>
#include <ctype.h>
#include <stddef.h>

/* GNU Guile library header file: */
#include <guile/gh.h>

/* Header files for my eigensolver routines: */
#include "../src/config.h"
#include <mpiglue.h>
#include <mpi_utils.h>
#include <check.h>
#include <blasglue.h>
#include <matrices.h>
#include <eigensolver.h>
#include <maxwell.h>

/* Header file for the ctl-file (Guile) interface; automatically
   generated from mpb.scm */
#include <ctl-io.h>

/* Routines from libctl/utils/geom.c: */
#include <ctlgeom.h>

/* shared declarations for the files in mpb-ctl: */
#include "mpb.h"

/*  matrix/field "smobs" (Scheme objects) */
#include "matrix-smob.h"
#include "field-smob.h"

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

/* use Guile malloc function instead of plain malloc, to force use
   of the garbage collector when we run low */
static void *malloc_hook(size_t n)
{
     return (void*) scm_must_malloc(n, "stuff for MPB");
}

/* The following are hook functions called from main() when
   starting the program and just before exiting.   We use
   them to initialize MPI. */

void ctl_start_hook(int *argc, char ***argv)
{
     my_malloc_hook = malloc_hook;
     MPI_Init(argc, argv);
}

void ctl_stop_hook(void)
{
     MPI_Finalize();
}

/* The following is a hook function called from main() when initializing
   Guile, which can export any additional symbols to Guile: */
void ctl_export_hook(void)
{
     register_matrix_smobs();
     register_field_smobs();
}

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

/* Some Guile-callable functions so that ctl files can know something
   about MPI. */

boolean mpi_is_masterp(void)
{
     return mpi_is_master();
}

boolean using_mpip(void)
{
#ifdef HAVE_MPI
     return 1;
#else
     return 0;
#endif
}

integer mpi_num_procs(void)
{
     int num_procs;
     MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
     return num_procs;
}

integer mpi_proc_index(void)
{
     int proc_num;
     MPI_Comm_rank(MPI_COMM_WORLD, &proc_num);
     return proc_num;
}

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

/* a couple of utilities to convert libctl data types to the data
   types of the eigensolver & maxwell routines: */

void vector3_to_arr(real arr[3], vector3 v)
{
     arr[0] = v.x;
     arr[1] = v.y;
     arr[2] = v.z;
}

void matrix3x3_to_arr(real arr[3][3], matrix3x3 m)
{
     vector3_to_arr(arr[0], m.c0);
     vector3_to_arr(arr[1], m.c1);
     vector3_to_arr(arr[2], m.c2);
}

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

/* global variables for retaining data about the eigenvectors between
   calls from Guile: */

int nwork_alloc = 0;

maxwell_data *mdata = NULL;
maxwell_target_data *mtdata = NULL;
evectmatrix H, W[MAX_NWORK], Hblock;

vector3 cur_kvector;
scalar_complex *curfield = NULL;
int curfield_band;
char curfield_type = '-';

void curfield_reset(void) { curfield = NULL; curfield_type = '-'; }

/* R[i]/G[i] are lattice/reciprocal-lattice vectors */
real R[3][3], G[3][3];
matrix3x3 Rm, Gm; /* same thing, but matrix3x3 */
real Vol; /* computational cell volume = |det Rm| */

/* index of current kpoint, for labeling output */
int kpoint_index = 0;

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

scalar_complex cnumber2cscalar(cnumber c)
{
     scalar_complex cs;
     CASSIGN_SCALAR(cs, cnumber_re(c), cnumber_im(c));
     return cs;
}

cnumber cscalar2cnumber(scalar_complex cs)
{
     return make_cnumber(CSCALAR_RE(cs), CSCALAR_IM(cs));
}

cvector3 cscalar32cvector3(const scalar_complex *cs)
{
     cvector3 v;
     v.x = cscalar2cnumber(cs[0]);
     v.y = cscalar2cnumber(cs[1]);
     v.z = cscalar2cnumber(cs[2]);
     return v;
}

void cvector32cscalar3(scalar_complex *cs, cvector3 v)
{
     cs[0] = cnumber2cscalar(v.x);
     cs[1] = cnumber2cscalar(v.y);
     cs[2] = cnumber2cscalar(v.z);
}

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

/* initialize the field to random numbers; should only be called
   after init-params.  (Guile-callable.) */
void randomize_fields(void)
{
     int i;

     if (!mdata)
	  return;
     mpi_one_printf("Initializing fields to random numbers...\n");
     for (i = 0; i < H.n * H.p; ++i) {
	  ASSIGN_SCALAR(H.data[i], rand() * 1.0 / RAND_MAX,
			rand() * 1.0 / RAND_MAX);
     }
}

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

/* Guile-callable functions for getting/setting the kpoint index. */

integer get_kpoint_index(void)
{
     return kpoint_index;
}

void set_kpoint_index(integer i)
{
     kpoint_index = i;
}

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

/* return a string describing the current parity, used for frequency
   and filename prefixes */
const char *parity_string(maxwell_data *d)
{
     static char s[128];
     strcpy(s, "");
     if (d->parity & EVEN_Z_PARITY)
	  strcat(s, (d->nz == 1) ? "te" : "zeven");
     else if (d->parity & ODD_Z_PARITY)
	  strcat(s, (d->nz == 1) ? "tm" : "zodd");
     if (d->parity & EVEN_Y_PARITY)
	  strcat(s, "yeven");
     else if (d->parity & ODD_Y_PARITY)
	  strcat(s, "yodd");
     return s;
}

/* Set the current parity to solve for. (init-params should have
   already been called.  (Guile-callable; see mpb.scm.) 

   p >= 0 means a bitwise OR of the various parity constants from
   maxwell.h (NO_PARITY, EVEN_Z_PARITY, etcetera).

   p = -1 means the parity of the previous call, 
       or NO_PARITY if this is the first call */

void set_parity(integer p)
{
     static int last_p = -2;  /* initialize to some non-value */

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

     if (p == -1)
	  p = last_p < 0 ? NO_PARITY : last_p;

     set_maxwell_data_parity(mdata, p);
     CHECK(mdata->parity == p, "k vector incompatible with parity");
     mpi_one_printf("Solving for band polarization: %s.\n",
		    parity_string(mdata));

     last_p = p;
     set_kpoint_index(0);  /* reset index */
}

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

/* Guile-callable function: init-params, which initializes any data
   that we need for the eigenvalue calculation.  When this function
   is called, the input variables (the geometry, etcetera) have already
   been read into the global variables defined in ctl-io.h.  
   
   p is the parity to use for the coming calculation, although
   this can be changed by calling set-parity.  p is interpreted
   in the same way as for set-parity.

   If reset_fields is false, then any fields from a previous run are
   retained if they are of the same dimensions.  Otherwise, new
   fields are allocated and initialized to random numbers. */
void init_params(integer p, boolean reset_fields)
{
     int i, local_N, N_start, alloc_N;
     int nx, ny, nz;
     int have_old_fields = 0;
     int block_size;
     
     /* Output a bunch of stuff so that the user can see what we're
	doing and what we've read in. */
     
     mpi_one_printf("init-params: initializing eigensolver data\n");
#ifndef SCALAR_COMPLEX
     mpi_one_printf("  -- assuming INVERSION SYMMETRY in the geometry.\n");
#endif
     
     mpi_one_printf("Computing %d bands with %e tolerance.\n",
		    num_bands, tolerance);
     if (target_freq != 0.0)
	  mpi_one_printf("Target frequency is %g\n", target_freq);
     
     get_grid_size_n(&nx, &ny, &nz);

     {
	  int true_rank = nz > 1 ? 3 : (ny > 1 ? 2 : 1);
	  if (true_rank < dimensions)
	       dimensions = true_rank;
	  else if (true_rank > dimensions) {
	       mpi_one_fprintf(stderr, 
			       "WARNING: rank of grid is > dimensions.\n"
			       "         setting extra grid dims. to 1.\n");
	       /* force extra dims to be 1: */
	       if (dimensions <= 2)
		    nz = 1;
	       if (dimensions <= 1)
		    ny = 1;
	  }
     }

     mpi_one_printf("Working in %d dimensions.\n", dimensions);
     mpi_one_printf("Grid size is %d x %d x %d.\n", nx, ny, nz);

     if (eigensolver_block_size != 0 && eigensolver_block_size < num_bands) {
	  block_size = eigensolver_block_size;
	  if (block_size < 0) {
	       /* Guess a block_size near -block_size, chosen so that
		  all blocks are nearly equal in size: */
	       block_size = (num_bands - block_size - 1) / (-block_size);
	       block_size = (num_bands + block_size - 1) / block_size;
	  }
	  mpi_one_printf("Solving for %d bands at a time.\n", block_size);
     }
     else
	  block_size = num_bands;

     if (mdata) {  /* need to clean up from previous init_params call */
	  if (nx == mdata->nx && ny == mdata->ny && nz == mdata->nz &&
	      block_size == Hblock.alloc_p && num_bands == H.p &&
	      eigensolver_nwork == nwork_alloc)
	       have_old_fields = 1; /* don't need to reallocate */
	  else {
	       destroy_evectmatrix(H);
	       for (i = 0; i < nwork_alloc; ++i)
		    destroy_evectmatrix(W[i]);
	       if (Hblock.data != H.data)
		    destroy_evectmatrix(Hblock);
	  }
	  destroy_maxwell_target_data(mtdata); mtdata = NULL;
	  destroy_maxwell_data(mdata); mdata = NULL;
	  curfield_reset();
     }
     else
	  srand(time(NULL)); /* init random seed for field initialization */
   
     if (deterministicp) {  /* check input variable "deterministic?" */
	  /* seed should be the same for each run, although
	     it should be different for each process: */
	  int rank;
	  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
	  srand(314159 * (rank + 1));
     }

     mpi_one_printf("Creating Maxwell data...\n");
     mdata = create_maxwell_data(nx, ny, nz, &local_N, &N_start, &alloc_N,
                                 block_size, NUM_FFT_BANDS);
     CHECK(mdata, "NULL mdata");

     if (target_freq != 0.0)
	  mtdata = create_maxwell_target_data(mdata, target_freq);
     else
	  mtdata = NULL;

     if (!have_old_fields) {
	  mpi_one_printf("Allocating fields...\n");
	  H = create_evectmatrix(nx * ny * nz, 2, num_bands,
				 local_N, N_start, alloc_N);
	  nwork_alloc = eigensolver_nwork;
	  for (i = 0; i < nwork_alloc; ++i)
	       W[i] = create_evectmatrix(nx * ny * nz, 2, block_size,
					 local_N, N_start, alloc_N);
	  if (block_size < num_bands)
	       Hblock = create_evectmatrix(nx * ny * nz, 2, block_size,
					   local_N, N_start, alloc_N);
	  else
	       Hblock = H;
     }

     init_epsilon();

     mpi_one_printf("%d k-points:\n", k_points.num_items);
     for (i = 0; i < k_points.num_items; ++i)
	  mpi_one_printf("     (%g,%g,%g)\n", k_points.items[i].x,
			 k_points.items[i].y, k_points.items[i].z);

     set_parity(p);