operator.cc

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

		  X.stack.rSize, X.stack.cfSize, X.stack.rfSize );
  X.offsetCopyStack( Z, 1 );
  Z.stack.com[0] = MINUS;
  return Z;
}

Operator Operator::operator*(const Operator& X) const
//
// Multiplication of two operators. Similar to `+'.
{
  if( stack.com[0] == UNINITIALIZED || X.stack.com[0] == UNINITIALIZED )
    error("Attempt to multiply 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] = TIMES;
  return Z;
}

Operator Operator::operator*(const Complex& alpha) const
///////////////////////////////////////////////////////////////////////////
// Scalar multiplication of an operator by a complex number.
// 
// Example: If `Complex alpha (1,1);' and the stack of the operator X 
// is 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      COMPLEX        (pointer to A)
//  0      TIMES          (pointer to B)       3 + 2i
// 
// then the stack of the `Operator Z = X * alpha;' is given by
// 
//  i    Z.stack.com[i]    Z.stack.op[i]     Z.stack.c[i]     Z.stack.r[i]
//   
//  6      OP
//  5      OP
//  4      PLUS
//  3      COMPLEX
//  2      TIMES
//  1      COMPLEX        (pointer to A)       3 + 2i
//  0      TIMES          (pointer to B)       1 + i
//
//  (`stack.cf' and `stack.rf' are omitted to keep the example simple.)
// 
/////////////////////////////////////////////////////////////////////////////
{
  if( stack.com[0] == UNINITIALIZED )
    error("Attempt to multiply an uninitialized operator by a scalar.");
  Operator Z;
  offsetCopyStack( Z, 2, 0, 1 );
  Z.stack.com[1] = COMPLEX;
  Z.stack.com[0] = TIMES;
  Z.stack.c[0] = alpha;
  return Z;
}

Operator Operator::operator*(double a) const
//
// Scalar multiplication of an operator by a real number.
{
  if( stack.com[0] == UNINITIALIZED )
    error("Attempt to multiply an uninitialized operator by a scalar.");
  Operator Z;
  offsetCopyStack( Z, 2, 0, 0, 1 );
  Z.stack.com[1] = REAL;
  Z.stack.com[0] = TIMES;
  Z.stack.r[0] = a;
  return Z;
}

Operator Operator::operator*(ImaginaryUnit im) const
//
// Scalar multiplication of an operator by the imaginary unit `IM' or 
// by its negative `M_IM'. `IM' and `M_IM' are defined in "State.h".
{
  if( stack.com[0] == UNINITIALIZED )
    error("Attempt to multiply an uninitialized operator by a scalar.");
  Operator Z;
  offsetCopyStack( Z, 2 );
  if( im == IM )
    Z.stack.com[1] = IMAG;
  else
    Z.stack.com[1] = M_IMAG;
  Z.stack.com[0] = TIMES;
  return Z;
}

Operator Operator::operator*(ComplexFunction f) const
//
// Scalar multiplication of an operator by a complex-valued function.
{
  if( stack.com[0] == UNINITIALIZED )
    error("Attempt to multiply an uninitialized operator by a function.");
  Operator Z;
  offsetCopyStack( Z, 2, 0, 0, 0, 1 );
  Z.stack.com[1] = CFUNC;
  Z.stack.com[0] = TIMES;
  Z.stack.cf[0] = f;
  return Z;
}

Operator Operator::operator*(RealFunction f) const
//
// Scalar multiplication of an operator by a real-valued function.
{
  if( stack.com[0] == UNINITIALIZED )
    error("Attempt to multiply an uninitialized operator by a function.");
  Operator Z;
  offsetCopyStack( Z, 2, 0, 0, 0, 0, 1 );
  Z.stack.com[1] = RFUNC;
  Z.stack.com[0] = TIMES;
  Z.stack.rf[0] = f;
  return Z;
}

Operator Operator::pow(int n) const
///////////////////////////////////////////////////////////////////////
// Integer power (n > 0) of an operator.
//
// Example: If the stack of the `Operator X = A + B;' is given by
//
//  i    X.stack.com[i]    X.stack.op[i]     X.stack.c[i]     X.stack.r[i]
// 
//  2      OP
//  1      OP              (pointer to A)
//  0      PLUS            (pointer to B)
//  
// then the stack of the `Operator Z = X.pow(3);' is given by
//
//  i    Z.stack.com[i]    Z.stack.op[i]     Z.stack.c[i]     Z.stack.r[i]
//
//  10     OP
//  9      OP
//  8      PLUS
//  7      OP
//  6      OP
//  5      PLUS            (pointer to A)
//  4      TIMES           (pointer to B)
//  3      OP              (pointer to A)
//  2      OP              (pointer to B)
//  1      PLUS            (pointer to A)
//  0      TIMES           (pointer to B)
//
//  (`stack.cf' and `stack.rf' are omitted to keep the example simple.)
// 
// Each of the statements
// psi *= Z;
// or
// psi *= (A+B).pow(3);
// leads to the execution of the following code:
// {
//   State psi1(psi);
//   psi1 *= B;
//   psi *= A;
//   psi += psi1;
// }     // Here psi1 goes out of scope and is destroyed.
// {
//   State psi1(psi);
//   psi1 *= B;
//   psi *= A;
//   psi += psi1;
// }     // Here psi1 goes out of scope and is destroyed.
// {
//   State psi1(psi);
//   psi1 *= B;
//   psi *= A;
//   psi += psi1;
// }     // Here psi1 goes out of scope and is destroyed.
////////////////////////////////////////////////////////////////////
{
  if(n <= 0) error("pow(n) not defined for n <= 0.");
  if( stack.com[0] == UNINITIALIZED )
    error("Attempt to compute a power of an uninitialized operator.");
  Operator Z;
  offsetCopyStack (Z, (n-1)*(stack.comSize+1), (n-1)*stack.opSize,
    (n-1)*stack.cSize,(n-1)*stack.rSize,(n-1)*stack.cfSize,(n-1)*stack.rfSize);
  for( int i=0; i<n-1; i++ ) {
    offsetCopyStack( Z, i*(stack.comSize+1)+1, i*stack.opSize, 
	   i*stack.cSize, i*stack.rSize, i*stack.cfSize, i*stack.rfSize );
    Z.stack.com[i*(stack.comSize+1)] = TIMES;
  }
  return Z;
}

Operator Operator::hc() const
//
// Hermitian conjugate of an operator.
{
  StackPtr stackPtr = {0,0,0,0,0,0};   // Initialize stack pointers.
  return dagger( stackPtr );           // Enter the recursive `dagger'.
}

Operator Operator::dagger( StackPtr& stackPtr ) const
{
  Command c=stack.com[stackPtr.com++];   // Take command-stack content and pop.
  switch (c) 
  {
  case OPERATOR: {
    #ifdef DEBUG_TRACE
    cout << "Operator::dagger for case OPERATOR entered." << endl;
    #endif
    Operator Z;
    Z.deallocate();
    Z.allocate(1,1);
    Z.stack.com[0] = OPERATOR_HC;
    Z.stack.op[0] = stack.op[stackPtr.op++];
    return Z;
  }
  case OPERATOR_HC: {
    #ifdef DEBUG_TRACE
    cout << "Operator::dagger for case OPERATOR_HC entered." << endl;
    #endif
    Operator Z;
    Z.deallocate();
    Z.allocate(1,1);
    Z.stack.com[0] = OPERATOR;
    Z.stack.op[0] = stack.op[stackPtr.op++];
    return Z;
  }
  case COMPLEX: {
    #ifdef DEBUG_TRACE
    cout << "Operator::dagger for case COMPLEX entered." << endl;
    #endif
    Operator Z;
    Z.deallocate();
    Z.allocate(1,0,1);
    Z.stack.com[0] = COMPLEX;
    Z.stack.c[0] = conj( stack.c[stackPtr.c++] );
    return Z;   // `Z' is not a well-formed operator, but the final result is.
  }
  case REAL: {
    #ifdef DEBUG_TRACE
    cout << "Operator::dagger for case REAL entered." << endl;
    #endif
    Operator Z;
    Z.deallocate();
    Z.allocate(1,0,0,1);
    Z.stack.com[0] = REAL;
    Z.stack.r[0] = stack.r[stackPtr.r++];
    return Z;
  }
  case IMAG: {
    #ifdef DEBUG_TRACE
    cout << "Operator::dagger for case IMAG entered." << endl;
    #endif
    Operator Z;
    Z.deallocate();
    Z.allocate(1);
    Z.stack.com[0] = M_IMAG;
    return Z;
  }
  case M_IMAG: {
    #ifdef DEBUG_TRACE
    cout << "Operator::dagger for case M_IMAG entered." << endl;
    #endif
    Operator Z;
    Z.deallocate();
    Z.allocate(1);
    Z.stack.com[0] = IMAG;
    return Z;
  }
  case CFUNC: {
    #ifdef DEBUG_TRACE
    cout << "Operator::dagger for case CFUNC entered." << endl;
    #endif
    Operator Z;
    Z.deallocate();
    Z.allocate(1,0,0,0,1);
    Z.stack.com[0] = CCFUNC;
    Z.stack.cf[0] = stack.cf[stackPtr.cf++];
    return Z;   // `Z' is not a well-formed operator, but the final result is.
  }
  case CCFUNC: {
    #ifdef DEBUG_TRACE
    cout << "Operator::dagger for case CCFUNC entered." << endl;
    #endif
    Operator Z;
    Z.deallocate();
    Z.allocate(1,0,0,0,1);
    Z.stack.com[0] = CFUNC;
    Z.stack.cf[0] = stack.cf[stackPtr.cf++];
    return Z;   // `Z' is not a well-formed operator, but the final result is.
  }
  case RFUNC: {
    #ifdef DEBUG_TRACE
    cout << "Operator::dagger for case RFUNC entered." << endl;
    #endif
    Operator Z;
    Z.deallocate();
    Z.allocate(1,0,0,0,0,1);
    Z.stack.com[0] = RFUNC;
    Z.stack.rf[0] = stack.rf[stackPtr.rf++];
    return Z;   // `Z' is not a well-formed operator, but the final result is.
  }
  case PLUS: {
    #ifdef DEBUG_TRACE
    cout << "Operator::dagger for case PLUS entered." << endl;
    #endif
    Operator Z = dagger( stackPtr );
    Operator X = dagger( stackPtr );
    return X += Z;
  }
  case MINUS: {
    #ifdef DEBUG_TRACE
    cout << "Operator::dagger for case MINUS entered." << endl;
    #endif
    Operator Z = dagger( stackPtr );
    Operator X = dagger( stackPtr );
    return X -= Z;
  }
  case TIMES: {
    #ifdef DEBUG_TRACE
    cout << "Operator::dagger for case TIMES entered." << endl;
    #endif
    Operator Z = dagger( stackPtr );
    Operator X = dagger( stackPtr );
    return Z * X;     // Notice the reversal of the order of terms.
  }
  case UNINITIALIZED:    
    error("Attempt to compute the conjugate of an uninitialized operator.");

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

//////////////////////////////////////////////////////////////////////
// The following algebraic operations on operators are all defined in
// terms of algebraic operations defined above.
//////////////////////////////////////////////////////////////////////

Operator operator-(const Operator& op)      // Unary `-'.
{
  Operator result = -1.0 * op;
  return result;
}

Operator operator+(const Operator& op)      // Unary `+'.
{
  Operator result = op;
  return result;
}

Operator& Operator::operator+=(const Operator& op)
{
  *this = *this + op;
  return *this;
}

Operator& Operator::operator-=(const Operator& op)
{
  *this = *this - op;
  return *this;
}

Operator& Operator::operator*=(const Complex& alpha)
{
  *this = *this * alpha;
  return *this;
}

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

Operator& Operator::operator*=(ImaginaryUnit im)
{
  *this = *this * im;
  return *this;
}

Operator operator*(const Complex& alpha, const Operator& op)
{
  return op * alpha;
}
  
Operator operator*(double a, const Operator& op)
{
  return op * a;
}
  
Operator operator*(ImaginaryUnit im, const Operator& op)
{
  return op * im;
}

Operator& Operator::operator*=(ComplexFunction f)
{
  *this = *this * f;
  return *this;
}

Operator& Operator::operator*=(RealFunction f)
{
  *this = *this * f;
  return *this;
}

Operator operator*(ComplexFunction f, const Operator& op)
{
  return op * f;
}
  
Operator operator*(RealFunction f, const Operator& op)
{
  return op * f;
}