operator.cc

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

//   Operator.cc -- Operator algebra in Hilbert space
//     
//   Copyright (C) 1995  Todd Brun and Ruediger Schack
//   
//   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., 675 Mass Ave, Cambridge, MA 02139, USA.
//
//   ----------------------------------------------------------------------
//   If you improve the code or make additions to it, or if you have
//   comments or suggestions, please contact us:
//
//   Dr. Todd Brun			        Tel    [phone]
//   Department of Physics                      FAX    [phone]
//   Queen Mary and Westfield College           email  [e-mail]
//   Mile End Road, London E1 4NS, UK
//
//   Dr. Ruediger Schack                        Tel    [phone]
//   Department of Mathematics                  FAX    [phone]
//   Royal Holloway, University of London       email  [e-mail]
//   Egham, Surrey TW20 0EX, UK
/////////////////////////////////////////////////////////////////////////////

#include <stdlib.h>
#include <iostream.h>

#include "Operator.h"
#include "PrimOp.h"

static const char rcsid[] = "$Id: Operator.cc,v 3.1 1996/11/19 10:05:08 rschack Exp $";

int Operator::flag = 0;   // Static flag, set by constructors, tested by
                          // PrimaryOperator class.

void Operator::error( const char* msg ) const
{
  cerr << "Fatal error in class Operator or derived class:\n  " << msg << endl;
  exit(1);
}

void Operator::allocate
  (int comSize,int opSize=0,int cSize=0,int rSize=0,int cfSize=0,int rfSize=0)
//
// Private function, used in object construction and assignment.
{
  if( comSize<0 || opSize<0 || cSize<0 || rSize<0 || cfSize<0 || rfSize<0 )
    error("Negative stack size in allocate.");
  stack.com = 0;
  stack.op = 0;
  stack.c = 0;
  stack.r = 0;
  stack.cf = 0;
  stack.rf = 0;
  if( (stack.comSize = comSize) > 0) stack.com = new Command[comSize];
  if( (stack.opSize = opSize) > 0) stack.op = new PrimaryOperator*[opSize];
  if( (stack.cSize = cSize) > 0) stack.c = new Complex[cSize];
  if( (stack.rSize = rSize) > 0) stack.r = new double[rSize];
  if( (stack.cfSize = cfSize) > 0) stack.cf = new ComplexFunction[cfSize];
  if( (stack.rfSize = rfSize) > 0) stack.rf = new RealFunction[rfSize];
}

void Operator::deallocate()
//
// Private function, used in object destruction and assignment.
{
  #ifndef NON_GNU_DELETE
    delete[] stack.com;
    delete[] stack.op;
    delete[] stack.c;
    delete[] stack.r;
    delete[] stack.cf;
    delete[] stack.rf;
  #else  //--------NON GNU CODE--------
    delete[stack.comSize] stack.com;
    delete[stack.opSize] stack.op;
    delete[stack.cSize] stack.c;
    delete[stack.rSize] stack.r;
    delete[stack.cfSize] stack.cf;
    delete[stack.rfSize] stack.rf;
  #endif //--------NON GNU CODE--------
}

void Operator::copy(const Operator& rhs)
//
// Private function, used in object construction and assignment.
{
  time = rhs.time;
  allocate( rhs.stack.comSize, rhs.stack.opSize, rhs.stack.cSize, 
	   rhs.stack.rSize, rhs.stack.cfSize, rhs.stack.rfSize );
  int i;
  for (i=0; i<stack.comSize; i++) stack.com[i] = rhs.stack.com[i];
  for (i=0; i<stack.opSize; i++) stack.op[i] = rhs.stack.op[i];
  for (i=0; i<stack.cSize; i++) stack.c[i] = rhs.stack.c[i];
  for (i=0; i<stack.rSize; i++) stack.r[i] = rhs.stack.r[i];
  for (i=0; i<stack.cfSize; i++) stack.cf[i] = rhs.stack.cf[i];
  for (i=0; i<stack.rfSize; i++) stack.rf[i] = rhs.stack.rf[i];
}

Operator::Operator()   // Default constructor.
{
  allocate(1);
  stack.com[0] = UNINITIALIZED;
  time = 0;
  flag = 1;
}

Operator::Operator(char*)  // Private constructor called by `PrimaryOperator'.
{                          // Does not set `flag'.               
  allocate(1,1);
  stack.com[0] = UNINITIALIZED;
  time = 0;
}

Operator::Operator(const Operator& rhs)   // Copy constructor.
{
  copy(rhs);
  flag = 1;
}

Operator::~Operator()    // Destructor.
{
  deallocate();
}

Operator& Operator::operator=(const Operator& rhs)   // Assignment.
{
  deallocate();
  copy(rhs);
  return *this;
}

void Operator::printCommandStack() const      // For debugging only.
{
  cout << "Operator::printCommandStack():" << endl;
  for( int i=0; i<stack.comSize; i++ )
    cout << stack.com[i] << endl;
}

Operator& Operator::operator()(double t)
{
  time = t;
  return *this;
}

State Operator::operator*(const State& psi) const
//
// Application of the operator (i.e., `*this') to the state `psi' using the 
// syntax `State psi0 = X * psi;' where `X' is the operator.
// Defined in terms of `*='. Creates a temporary `State' object and is
// therefore inefficient. Use `*=' wherever possible.
{
  State psi1(psi);
  return psi1 *= *this;
}

State& operator*=(State& psi, const Operator& X)
//
// Application of the operator `X' to the state `psi' using the 
// syntax `psi *= X;'. The state `psi' is modified; no extra `State'
// object is created. 
{
  #ifdef DEBUG_TRACE
    cout << "operator*=(State&,Operator&) entered." << endl;
  #endif
  #ifndef NON_GNU_SCOPE
    Operator::StackPtr stackPtr = {0,0,0,0,0,0};  // Initialize stack pointers.
  #else  //----------NON GNU CODE----
    StackPtr stackPtr = {0,0,0,0,0,0};
  #endif //----------NON GNU CODE----
  X.eval( stackPtr, psi );             // Enter the recursive `Operator::eval'.
  #ifdef DEBUG_TRACE
    cout << "Returning from operator*=(State&,Operator&)." << endl;
  #endif
  return psi;
}

void Operator::eval( StackPtr& stackPtr, State& psi ) const
/////////////////////////////////////////////////////////////////////////// 
// Private function.
// Applies the "popped" stack `stack' to the state `psi', thereby modifying
// `psi', and pops the stack by incrementing the stack pointer `stackPtr'.
// "Popped stack" means -- e.g., for the command stack `stack.com' --
// that the stack level given by the value of `stackPtr.com' is interpreted
// as the bottom level of the stack.
// 
// Example:  The assignment  `C=(A+B)*4.7;', where `A' and `B' are 
// primary operators, leads to the following contents of `C.stack':
// 
//  i    C.stack.com[i]    C.stack.op[i]     C.stack.c[i]     C.stack.r[i]
// 
//  4      OP
//  3      OP
//  2      PLUS
//  1      REAL           (pointer to A)
//  0      TIMES          (pointer to B)                          4.7
//
//  (`stack.cf' and `stack.rf' are omitted to keep the example simple.)
// 
// The statement 
// psi *= C;
// then leads to the execution of the following piece of code:
// {
//   psi *= 4.7;
//   State psi1;
//   psi1.xerox(psi);
//   psi1 *= B;
//   psi *= A;
//   psi += psi1;
// }
//////////////////////////////////////////////////////////////////////////////
{
  Command c=stack.com[stackPtr.com++];   // Take command-stack content and pop.
  switch (c) {
    
  case OPERATOR: {
    #ifdef DEBUG_TRACE
    cout << "Operator::eval for case OPERATOR entered." << endl;
    #endif
    PrimaryOperator& A = *stack.op[stackPtr.op++];
    psi.apply(A,NO_HC,A.myFreedom,A.myType,time);
    // Apply to `psi' the primary operator to which `stack.op' points.
    // Pop `stack.op'
    break;
  }

  case OPERATOR_HC: {
    #ifdef DEBUG_TRACE
    cout << "Operator::eval for case OPERATOR_HC entered." << endl;
    #endif
    PrimaryOperator& A = *stack.op[stackPtr.op++];
    psi.apply(A,HC,A.myFreedom,A.myType,time);
    // Apply to `psi' the Hermitian conjugate of the primary operator to which
    // `stack.op' points. Pop `stack.op'
    break;
  }

  case COMPLEX:                     // Multiply by the content of `stack.c'.
    psi *= stack.c[stackPtr.c++];   // Pop `stack.c'.
    break;                          // `*=' defined in class State.

  case REAL:                        // Multiply by the content of `stack.r'.
    psi *= stack.r[stackPtr.r++];   // Pop `stack.r'.
    break;                          // `*=' defined in class State.

  case IMAG:            // Multiply by the imaginary unit `IM'.
    psi *= IM;          // `*=' defined in class State.
    break;

  case M_IMAG:          // Multiply by the negative imaginary unit `M_IM'.
    psi *= M_IM;        // `*=' defined in class State.
    break;

  case CFUNC:  // Multiply by the function in `stack.cf' evaluated at time t.
    psi *= stack.cf[stackPtr.cf++](time); // Pop `stack.cf'.
    break;                                // `*=' defined in class State.

  case CCFUNC:  // Multiply by the conjugate of `stack.cf' evaluated at time t.
    psi *= conj(stack.cf[stackPtr.cf++](time)); // Pop `stack.cf'.
    break;                                      // `*=' defined in class State.

  case RFUNC:   // Multiply by the function in `stack.rf' evaluated at time t.
    psi *= stack.rf[stackPtr.rf++](time);  // Pop `stack.rf'.
    break;                                 // `*=' defined in class State.

  case PLUS: {
    #ifdef DEBUG_TRACE
    cout << "Operator::eval for case PLUS entered." << endl;
    #endif
    State psi1;                     
    psi1.xerox(psi);                // Copy psi into psi1.
    eval( stackPtr, psi1 );         // Apply (popped) stack to psi1.
    eval( stackPtr, psi );          // Apply (popped) stack to psi.
    psi += psi1;                    // Add the results.
    break;                          // `+=' defined in class State.
  }  // `psi1' goes out of scope.
    
  case MINUS: {
    #ifdef DEBUG_TRACE
    cout << "Operator::eval for case MINUS entered." << endl;
    #endif
    State psi1;                     
    psi1.xerox(psi);                // Copy psi into psi1.
    eval( stackPtr, psi1 );         // Apply (popped) stack to psi1.
    eval( stackPtr, psi );          // Apply (popped) stack to psi.
    psi -= psi1;                    // Subtract the results.
    break;                          // `-=' defined in class State.
  }  // `psi1' goes out of scope.

  case TIMES:
    #ifdef DEBUG_TRACE
    cout << "Operator::eval for case TIMES entered." << endl;
    #endif
    eval( stackPtr, psi );     // Apply (popped) stack to psi.
    eval( stackPtr, psi );     // Apply (popped) stack to (modified) psi.
    break;

  case UNINITIALIZED:    
    error("Attempt to apply an uninitialized operator to a state.");

  default:
    error("eval: Unknown object in command stack.");
  }
}

void Operator::offsetCopyStack (Operator& target,
      int n_com,int n_op=0,int n_c=0,int n_r=0,int n_cf=0,int n_rf=0) const
//
// Private function, used in `+', `-', etc. Copies `stack' into
// `target.stack', offset by integers n_com, n_op, n_c, n_r, n_cf, and n_rf.
{
  int newComSize = stack.comSize + n_com;
  int newOpSize = stack.opSize + n_op;
  int newCSize = stack.cSize + n_c;
  int newRSize = stack.rSize + n_r;
  int newCfSize = stack.cfSize + n_cf;
  int newRfSize = stack.rfSize + n_rf;
  if( target.stack.comSize < newComSize || 
     target.stack.opSize < newOpSize || 
     target.stack.cSize < newCSize ||
     target.stack.rSize < newRSize ||
     target.stack.cfSize < newCfSize ||
     target.stack.rfSize < newRfSize )
    {                        // Allocate a new stack to `target' if necessary.
      target.deallocate();
      target.allocate
	(newComSize,newOpSize,newCSize,newRSize,newCfSize,newRfSize);
    }

  int i;
  for( i=0; i<stack.comSize; i++ )          // Copy stack.com offset by n_com.
  target.stack.com[i+n_com] = stack.com[i];
  
  for( i=0; i<stack.opSize; i++ )           // Copy stack.op offset by n_op.
    target.stack.op[i+n_op] = stack.op[i];

  for( i=0; i<stack.cSize; i++ )            // Copy stack.c offset by n_c.
    target.stack.c[i+n_c] = stack.c[i];

  for( i=0; i<stack.rSize; i++ )            // Copy stack.r offset by n_r.
    target.stack.r[i+n_r] = stack.r[i];

  for( i=0; i<stack.cfSize; i++ )           // Copy stack.cf offset by n_cf.
    target.stack.cf[i+n_cf] = stack.cf[i];

  for( i=0; i<stack.rfSize; i++ )           // Copy stack.rf offset by n_rf.
    target.stack.rf[i+n_rf] = stack.rf[i];
}

Operator Operator::operator+(const Operator& X) const
////////////////////////////////////////////////////////////////////////////
// Addition of two operators.
// 
// Example: If the stacks of the operators X and Y are given by
// 
//  i    X.stack.com[i]    X.stack.op[i]     X.stack.c[i]     X.stack.r[i]
// 
//  4      OP
//  3      OP
//  2      PLUS
//  1      REAL           (pointer to A)
//  0      TIMES          (pointer to B)                          4.7
// 
// and
// 
//  i    Y.stack.com[i]    Y.stack.op[i]     Y.stack.c[i]     Y.stack.r[i]
// 
//  2      OP
//  1      REAL
//  0      TIMES          (pointer to C)                          -1.0
//  
// then the stack of the `Operator Z = X + Y;' is given by 
// 
//  i    Z.stack.com[i]    Z.stack.op[i]     Z.stack.c[i]     Z.stack.r[i]
//  
//  8      OP
//  7      OP
//  6      PLUS
//  5      REAL
//  4      TIMES
//  3      OP
//  2      REAL           (pointer to A)
//  1      TIMES          (pointer to B)                           4.7
//  0      PLUS           (pointer to C)                          -1.0
//
//  (`stack.cf' and `stack.rf' are omitted to keep the example simple.)
// 
////////////////////////////////////////////////////////////////////////////
{
  if( stack.com[0] == UNINITIALIZED || X.stack.com[0] == UNINITIALIZED )
    error("Attempt to add an uninitialized operator.");
  Operator Z;
  offsetCopyStack( Z, X.stack.comSize+1, X.stack.opSize, X.stack.cSize, 
		  X.stack.rSize, X.stack.cfSize, X.stack.rfSize );
  X.offsetCopyStack( Z, 1 );
  Z.stack.com[0] = PLUS;
  return Z;
}

Operator Operator::operator-(const Operator& X) const
//
// Subtraction of two operators. Similar to `+'.
{
  if( stack.com[0] == UNINITIALIZED || X.stack.com[0] == UNINITIALIZED )
    error("Attempt to subtract an uninitialized operator.");
  Operator Z;
  offsetCopyStack( Z, X.stack.comSize+1, X.stack.opSize, X.stack.cSize,