state.cc

使用量子轨道方法计算量子主方程的C++库

#endif
  free();				// recycle memory
  totalDim = maxSize = mySize = n;	// set up state
  nFreedoms = 1;
  nSkip = 1;
  coord = 0;
  myType = FIELD;
  data = new Complex[totalDim];		// allocate memory
  myPointer = data;
  nDims = 0;
  sizes = 0;
  nSkips = 0;
  partDims = 0;
  freedomTypes = 0;
  coords = 0;
  for (int i=0; i<totalDim; i++) data[i] = 0;
  data[nstate] = 1;					// Fock state
}

void State::coherent(int n, Complex alpha)	// Create a coherent state
{
#ifndef OPTIMIZE_QSD
  if( n < 2 )
    error("Invalid state size in State::coherent!");	 // exit on error
  if( norm(alpha) > n )
    error("Coherent state too big for storage in State::coherent!");
#endif
  Complex z;
  free();				// recycle memory
  totalDim = maxSize = mySize = n;	// set up state
  nFreedoms = 1;
  nSkip = 1;
  coord = 0;
  myType = FIELD;
  data = new Complex[totalDim];		// allocate memory
  myPointer = data;
  nDims = 0;
  sizes = 0;
  nSkips = 0;
  partDims = 0;
  freedomTypes = 0;
  coords = 0;
  z = exp( - 0.5*norm(alpha));		// normalization factor
  for (int i=0; i<totalDim; i++) {
    data[i] = z;			// coherent state elements
    z *= alpha/sqrt(i+1);
  }
}

void State::productState(int n, State* stateList)	// Create a product state
{
  int i,j,k,l;
  Complex element;
#ifndef OPTIMIZE_QSD
  if( n < 1 ) error("Invalid number of degrees of freedom in productState!");
#endif
  free();
  nFreedoms = n;
  nDims = new int[n];			// allocate memory
  sizes = new int[n];
  nSkips = new int[n];
  partDims = new int[n];
  freedomTypes = new FreedomType[n];
  coords = new Complex[n];
  totalDim = 1;
  for (i=0; i<n; i++) {			// set up memory structure for data
#ifndef OPTIMIZE_QSD
    if( stateList[i].nFreedoms != 1 )
      error("Can only produce products of 1D states!");
    if( stateList[i].totalDim < 2 )
      error("Invalid dimension size in productState!");
#endif
    nSkips[i] = totalDim;
    totalDim *= (nDims[i] = stateList[i].totalDim);
    sizes[i] = nDims[i];
    partDims[i] = totalDim;
    freedomTypes[i] = stateList[i].myType;
    coords[i] = stateList[i].coord;
  }
  mySize = totalDim;
  maxSize = totalDim;
  nSkip = 1;				// initialize state variables
  coord = coords[0];
  myType = freedomTypes[0];
  data = new Complex[totalDim];
  myPointer = data;
  for (i=0; i<totalDim; i++) {
    k = i;
    element = 1;
    for (j=n-1; j>=0; j--) {
      l = k/nSkips[j];			// index i_j of jth state
      element *= stateList[j][l];	// product state
      if( element == 0 ) break;		// if zero, leave this inner loop
      k %= nSkips[j];
    }
    data[i] = element;	// psi_1[i_1]*psi_2[i_2]*psi_3[i_3]*...*psi_n[i_n]
  }
}

State& State::operator=(const State& a)		// assignment operator
{
  if (this != &a) {	// otherwise a=a would lead to trouble
    if( (totalDim != a.totalDim) || (nFreedoms != a.nFreedoms) )
      free();		// recycle any memory
    copy(a);		// copy a into b
  }
  return *this;
}

State& State::operator=(int n)		// assign zero operator
{
  if (n != 0) error("Cannot assign a State a value other than zero!");
  for (int i=0; i<maxSize; i++) data[i] = 0;
  return *this;
}

Complex& State::operator[](int* indexList)
//
// Multiple degree of freedom subscript operator
// DO NOT USE THIS WITH UNNAMED TEMPORARIES
//
{
  int i,j=0;
  if (nFreedoms == 1)
    return data[indexList[0]];		// 1 degree of freedom case...
  for (i=0; i<nFreedoms; i++) j += indexList[i]*nSkips[i];
#ifndef OPTIMIZE_QSD
  if( (j>maxSize) || (j < 0) )
    error("Illegal arguments passed to function State::operator[]!");
#endif
  return data[j];
}

Complex State::elem(const int* indexList) const
//
// Multiple degree of freedom subscript operator.  Use this function
// when writing to the state might cause an error (e.g., when reading
// an element from an unnamed temporary).
//
{
  int i,j=0;
  if (nFreedoms == 1)
    return data[indexList[0]];		// 1 degree of freedom case...
  for (i=0; i<nFreedoms; i++) j += indexList[i]*nSkips[i];
#ifndef OPTIMIZE_QSD
  if( (j>maxSize) || (j < 0) )
    error("Illegal arguments passed to function elem!");
#endif
  return data[j];
}

Complex State::operator*(const State& a) const		// inner product
{
#ifndef OPTIMIZE_QSD
  if( mySize != a.mySize )
    error("Incompatible state sizes in State::operator*!");
  if( myType != a.myType )
    error("Incompatible state types in State::operator*!");
#endif
  int i,j,k;
  Complex sum = 0.0;
  for (i=0,j=0,k=0; i<mySize; i++,j+=nSkip,k+=a.nSkip)
    sum += conj(myPointer[j])*a.myPointer[k];
  return sum;
}

State& State::operator*=(const Complex& z)	// multiplication by Complex
{
  int i,j;
  for (i=0,j=0; i<mySize; i++,j+=nSkip)
    data[j] *= z;
  return *this;
}

State& State::operator*=(double x)		// multiplication by real
{
  int i,j;
  for (i=0,j=0; i<mySize; i++,j+=nSkip)
    data[j] *= x;
  return *this;
}

State& State::operator*=(ImaginaryUnit im)	// multiplication by i or -i
{
  int i,j;
  if( im == IM )
    for (i=0,j=0; i<mySize; i++,j+=nSkip) data[j].timesI();
  else
    for (i=0,j=0; i<mySize; i++,j+=nSkip) data[j].timesMinusI();
  return *this;
}

State& State::operator+=(const State& a)	// add state
{
#ifndef OPTIMIZE_QSD
  if( mySize != a.mySize )
    error("Incompatible state sizes in State::operator+=!");
  if( myType != a.myType )
    error("Incompatible state types in State::operator+=!");
#endif
  int i,j,k;
  for (i=0,j=0,k=0; i<mySize; i++,j+=nSkip,k+=a.nSkip)
    data[j] += a.data[k];
  return *this;
}

State& State::operator-=(const State& a)	// subtract state
{
#ifndef OPTIMIZE_QSD
  if( mySize != a.mySize )
    error("Incompatible state sizes in State::operator+=!");
  if( myType != a.myType )
    error("Incompatible state types in State::operator+=!");
#endif
  int i,j,k;
  for (i=0,j=0,k=0; i<mySize; i++,j+=nSkip,k+=a.nSkip)
    data[j] -= a.data[k];
  return *this;
}

void State::normalize()			// normalize state to 1
{
  int i;
  double sum=0.0;
  double epsLimit;
  epsLimit = smallSize/maxSize;
  for (i=0; i<maxSize; i++) {
    sum += norm(data[i]);
  }
  if(sum > epsLimit ) {			// check if zero state...
    sum = 1.0/sqrt(sum);
    for (i=0; i<maxSize; i++)
      data[i] *= sum;
  }
}

double State::checkBounds(int theFreedom, int numChecks)
//
// Check how much probability is in the top numChecks levels
// of the degree of freedom theFreedom.
//
{
#ifndef OPTIMIZE_QSD
  if( (theFreedom < 0) || (theFreedom >= nFreedoms) )
    error("Illegal freedom passed to State::checkBounds!");
  if( (nFreedoms == 1) && (numChecks > maxSize) )
    error("Too many checks requested in State::checkBounds!");
  if( (nFreedoms > 1) && (numChecks > sizes[theFreedom]) )
    error("Too many checks requested in State::checkBounds!");
#endif
  double sum = 0.0;
  int i,j,k,l;
  if( nFreedoms == 1) {				// 1 freedom case
    for(i=1; i<=numChecks; i++) {
      sum += norm(data[maxSize - i]);		// add up probability
    }
  }
  else {					// multi-freedom case
    int theSkip = nSkips[theFreedom];
    int theBound = (sizes[theFreedom] - 1)*theSkip;
    int theDim = partDims[theFreedom];
    for(i=0; i<theSkip; i++)			// loop over other indices
    for(j=0; j<maxSize; j+=theDim)
    for(k=0,l=theBound; k<numChecks; k++,l-=theSkip) {
      sum += norm(data[i+j+l]);			// add up probability
    }
  }
  return sum;
}

void State::fullSize()
//
// Makes the cutoffs of the state match the physical size in memory
//
{
  if( nFreedoms == 1 ) {
    for (int i=maxSize; i<totalDim; i++)
      data[i] = 0.0;
    maxSize = totalDim;
    mySize = totalDim;
  }
  else {
    for(int i=0; i<nFreedoms; i++)
      if( sizes[i] < nDims[i] ) stretchFreedom(i,nDims[i]);
  }
}

void State::adjustCutoff(int theFreedom, double epsilon, int padSize)
//
// This member function adjusts the amount of memory actually being
// used by a State.  As very frequently only the lowest few states have
// significant probability, it is more efficient not to loop over the
// top states at all.  This routine eliminates the top N states (which
// contain total probability less than epsilon) with a pad of size padSize
// for safety.
//
{
#ifndef OPTIMIZE_QSD
  if( (theFreedom < 0) || (theFreedom >= nFreedoms) )
    error("Illegal freedom passed to State::checkBounds!");
  if( (nFreedoms == 1) && (padSize > totalDim) )
    error("Too large a padSize requested in State::adjustCutoff!");
  if( (nFreedoms > 1) && (padSize > nDims[theFreedom]) )
    error("Too large a padSize requested in State::adjustCutoff!");
  if( padSize < 1 ) error("padSize must be positive in State::adjustCutoff!");
  if( (epsilon < 0.0) || (epsilon > 1.0) )
    error("Illegal value of epsilon in State::adjustCutoff!");
  if( (nFreedoms == 1) && (myType != FIELD) )
    error("Can only adjust cutoff size for fields!");
  if( (nFreedoms > 1) && (freedomTypes[theFreedom] != FIELD) )
    error("Can only adjust cutoff size for fields!");
#endif
  double sum = 0.0;
  int i,j,l;
  normalize();				// make sure total prob. is 1
  if( nFreedoms == 1) {			// 1 freedom case
    for(i=maxSize; sum < epsilon;)
      sum += norm(data[--i]);		// add up probability
    i += (padSize + 1);
    if( i > totalDim ) {		// enough states allocated?
      cerr << "Warning!  Significant probability in top states!" << endl;
      i = totalDim;
    }
    if( i > maxSize )				// stretch # of states;
      for(j=maxSize; j<i; j++) data[j] = 0;	// zero top states
    maxSize = i;
    mySize = maxSize;
  }
  else {				// multi-freedom case
    int theSkip = nSkips[theFreedom];
    int theBound = (sizes[theFreedom] - 1)*theSkip;
    int theDim = partDims[theFreedom];
    int k;
// sum over other indices
    for(k=sizes[theFreedom],l=theBound; sum < epsilon; k--,l-=theSkip)
    for(i=0; i<theSkip; i++)
    for(j=0; j<maxSize; j+=theDim) {
      sum += norm(data[i+j+l]);		// add up probability
    }
    k += (padSize + 1);
    if( k > nDims[theFreedom] ) {	// enough states allocated?
      cerr << "Warning!  Significant probabilty in top states!" << endl;
      cerr << "k = " << k << endl;
      k = nDims[theFreedom];
    }
    if( k > sizes[theFreedom] ) stretchFreedom(theFreedom,k);
    if( k < sizes[theFreedom] ) shrinkFreedom(theFreedom,k);
  }
}

void State::stretchFreedom(int theFreedom, int theSize)
//
// If the space being used for a degree of freedom is not sufficient to
// avoid truncation errors, stretch the available memory.
//
{
  int i,j,k;
  int newDim,newSize;

// Stretch the data into a larger space
  newDim = nSkips[theFreedom]*theSize;		// new partial dimension
  newSize = theSize*maxSize/sizes[theFreedom];	// new total size
  for(i=maxSize-partDims[theFreedom],j=newSize-newDim; i>=0;
       i-=partDims[theFreedom],j-=newDim) {
    for(k=partDims[theFreedom]-1; k>=0; k--)
      data[j+k] = data[i+k];			// stretch data
    for(k=newDim-1; k>=partDims[theFreedom]; k--)
      data[j+k] = 0;				// fill empty spaces w/zero
  }
// Recompute nSkips and partDims
  sizes[theFreedom] = theSize;			// size of freedom
  partDims[theFreedom] = newDim;
  maxSize = newSize;				// total size used
  for(i=theFreedom+1; i<nFreedoms; i++) {
    nSkips[i] = partDims[i-1];			// new nSkips
    partDims[i] = nSkips[i]*sizes[i];		// new partDims
  }
  mySize = maxSize;
}

void State::shrinkFreedom(int theFreedom, int theSize)
//
// If more memory is being used than is necessary to avoid truncation
// error, shrink memory used by freedom.
//
{
  int i,j,k;
  int newDim,newSize;

// Compact the data into a smaller space
  newDim = nSkips[theFreedom]*theSize;		// new partial dimension
  newSize = theSize*maxSize/sizes[theFreedom];	// new total size
  for(i=0,j=0; i<maxSize; i+=partDims[theFreedom],j+=newDim) //old and new outer indices
  for(k=0; k<newDim; k++)			// inner indices
    data[j+k] = data[i+k];			// compact data
// Pad out the empty space with zeros
  for(i=newSize; i<maxSize; i++) data[i] = 0;
// Recompute nSkips and partDims
  sizes[theFreedom] = theSize;			// size of freedom
  partDims[theFreedom] = newDim;
  maxSize = newSize;				// total size used
  for(i=theFreedom+1; i<nFreedoms; i++) {
    nSkips[i] = partDims[i-1];			// new nSkips
    partDims[i] = nSkips[i]*sizes[i];		// new partDims
  }
  mySize = maxSize;
}

void State::diagnostic()	// for debugging
{
  int i;
  cout << "Number of Freedoms " << nFreedoms << ".\n";
  cout << "Total Dimensions " << totalDim << ".\n";
  cout << "MaxSize " << maxSize << ".\n";
  cout << "My Type " << myType << ".\n";
  cout << "My Size " << mySize << ".\n";
  cout << "nSkip " << nSkip << ".\n";
  cout << "coord " << coord << ".\n";