Added BLAS/LAPACK support for matrix inverse inv()

src/Floquet/Floq2Op.cc

@@ -803,6 +803,18 @@ floq2_op exp(floq2_op &Op1)
   floq2_op EXpOp1(Op1.N1, Op1.N2, Op1.hs, Op1._omega1, Op1._omega2);
    EXpOp1.gen_op::operator= (exp((gen_op&) Op1));  // calculate exp(Op1)
   return EXpOp1;
+  }
+
+
+	// Input		Op1   : Floquet operator
+        // Return		Op    : exponential of Op1 using Pade
+	//				Op = exp(Op1)
+
+floq2_op expm(floq2_op &Op1)
+  {  
+  floq2_op EXpOp1(Op1.N1, Op1.N2, Op1.hs, Op1._omega1, Op1._omega2);
+   EXpOp1.gen_op::operator= (Op1.gen_op::expm());  // calculate exp(Op1)
+  return EXpOp1;
   }
 
 /*

src/Floquet/Floq2Op.h

@@ -243,13 +243,19 @@ MSVCDLL        void     operator /= (      double   d);
         // Output          gen_op    : Operator containing trace(s) over    
         //                             Photon space
 
-MSVCDLL friend floq2_op exp (floq2_op &Op1);
+MSVCDLL friend floq2_op exp(floq2_op &Op1);
 
         // Input		Op1  : Floquet operator
         // Return		Op   : exponential of Op1
 	//			       Op = exp(Op1)
         // Note			     : Computed in EBR of Op1
 
+MSVCDLL friend floq2_op expm(floq2_op &Op1);
+
+        // Input		Op1  : Floquet operator
+        // Return		Op   : exponential of Op1 using Pade
+	//			       Op = exp(Op1)
+
 /*  friend floq2_op prop (floq2_op &FLOQHAM , double &time);
 
 	// Input		Op1   : Floquet Hamiltonian

src/Floquet/FloqOp.cc

@@ -653,6 +653,32 @@ floq_op floq_op::exp()
   return ExpFOp;
   }
 
+        // Input		FOp	: Floquet operator
+        // Return		ExpFOp	: Exponential of FOp using Pade
+	//				   ExpFOp = exp(FOp)
+
+floq_op expm(const floq_op& FOp)
+  {
+  floq_op ExpFOp(FOp.N, FOp.hs, FOp.Om);	// Allocate return operator
+  ExpFOp.gen_op::operator= (FOp.gen_op::expm());	// Set exponential array
+  return ExpFOp;				// Return exponentiated FOp
+  }
+
+
+floq_op floq_op::expm()
+
+        // Input                FlOp    : Floquet operator (*this)
+        // Return               ExpFlOp : Exponential of Op1 using Pade
+        //                                   Op = exp(Op1)
+	// sosi				: Be nice to set this up as in
+	//				  class gen_op.....
+
+  {
+  floq_op ExpFOp(N, hs, Om);	// Allocate return operator
+  ExpFOp=*this;
+  ExpFOp.gen_op::operator= (ExpFOp.gen_op::expm());	// Set exponential array
+  return ExpFOp;
+  }
 
 floq_op prop(floq_op& FH, double time)
 

src/Floquet/FloqOp.h

@@ -268,6 +268,18 @@ MSVCDLL floq_op exp();
         // Return               ExpFlOp : Exponential of Op1
         //                                   Op = exp(Op1)
 
+MSVCDLL friend floq_op expm(const floq_op& Op1);
+
+        // Input		Op1  : Floquet operator
+        // Return		Op   : exponential of Op1 using Pade
+	//			       Op = exp(Op1)
+
+MSVCDLL floq_op expm();
+
+        // Input                FlOp    : Floquet operator (*this)
+        // Return               ExpFlOp : Exponential of Op1 using Pade
+        //                                   Op = exp(Op1)
+
 MSVCDLL friend floq_op prop(floq_op& FLOQHAM, double time);
 
 	// Input		Op1   : Floquet Hamiltonian

src/HSLib/GenOp.cc

@@ -1165,13 +1165,14 @@ complex gen_op::proj(const gen_op &Op2, int normf) const
   }
 
 
+	// Input		Op   : General operator (this)
+        // Return		int  : Op Liouville space dimension
+
 int gen_op::dim()    const { return (!WBR)?0:WBR->RepMx.cols(); }
 int gen_op::HS()     const { return (!WBR)?0:WBR->RepMx.cols(); }
 int gen_op::LS()     const { return (!WBR)?0:WBR->RepBs.dim_LS(); }
 int gen_op::dim_LS() const { return (!WBR)?0:WBR->RepBs.dim_LS(); }
 
-	// Input		Op   : General operator (this)
-        // Return		int  : Op Liouville space dimension
 
 
 	// Input		Op1   : General operator (this)
@@ -1228,6 +1229,52 @@ gen_op gen_op::exp(const complex& t, double cutoff) const
   return ExpOp;
   }
 
+	// Input		Op   : General operator (this)
+        // Return		int  : Op Liouville space dimension
+
+
+	// Input		Op1   : General operator (this)
+        // Return		Op    : exponential of Op1
+	//				Op = exp(Op1)
+        // Note			      : Computed in same base as Op1
+        // Note			      : Returns I matrix for Null Op
+
+gen_op gen_op::expm() const
+  {
+  int hs = dim();				// Operator Hilbert space
+  if(!WBR)					// If we have no reps then
+    {						// return an identity matrix
+    if(!hs) GenOpfatality(3, "exp"); 		// or Op exponential error
+    return gen_op(matrix(hs,hs,i_matrix_type)); // if there is no dimension
+    }
+  gen_op ExpOp(*this); 				// Copy Op in diagonal form
+  matrix mx1 = this->get_mx();
+  ExpOp.put_mx(mx1.expm());
+  return ExpOp;
+  }     
+
+
+
+        // Input                Op    : Operator (this)
+        //                      t     : Exponential factor
+        //                      cutoff: Exponential factor roundoff
+        // Return               ExpOp : Exponential of Op
+        //                              ExpOp = exp(t*Op)
+        // Note                       : Exponential output in same base as Op
+        // Note                       : Value of t is considered 0 if
+        //                              it's magnituded is less than cutoff
+
+gen_op gen_op::expm(const complex& t, double cutoff) const
+  {
+  int hs = dim();				// Operator Hilbert space
+  if(!WBR && !hs) GenOpfatality(3,"exp"); 	// Op exponential error
+  if(!WBR || norm(t) < fabs(cutoff))
+    return gen_op(matrix(hs, hs, i_matrix_type));
+  gen_op ExpOp(*this); 
+  matrix mx1 = this->get_mx()*t;
+  ExpOp.put_mx(mx1.expm());
+  return ExpOp;
+  }
 
 
         // Input                Op      : Operator (this)
@@ -1682,9 +1729,14 @@ void gen_op::status(int pf) const
       case i_matrix_type:
       case d_matrix_type:
       case n_matrix_type:
-        if((REP->RepMx).test_hermitian()) std::cout << ", Hermitian"; 
-        else                              std::cout << ", Non-Hermitian"; break;
-      default:                            std::cout << ", Non-Hermitian"; break;
+        if((REP->RepMx).test_hermitian())
+	  std::cout << ", Hermitian"; 
+        else                              
+	  std::cout << ", Non-Hermitian";
+	break;
+      default:                            
+	std::cout << ", Non-Hermitian";
+	break;
       }
 
     std::cout << "\n\t     - basis  = " << (REP->RepBs).refs() << " refs";
@@ -1702,9 +1754,14 @@ void gen_op::status(int pf) const
       case i_matrix_type:
       case d_matrix_type:
       case n_matrix_type:
-        if((REP->RepMx).test_hermitian()) std::cout << ", Hermitian";
-        else                              std::cout << ", Non-Hermitian"; break;
-      default:                            std::cout << ", Non-Hermitian"; break;
+        if((REP->RepMx).test_hermitian())
+	  std::cout << ", Hermitian";
+        else                             
+	  std::cout << ", Non-Hermitian";
+	break;
+      default:                           
+	std::cout << ", Non-Hermitian";
+	break;
       }
     REP++;
     item++;
@@ -1904,8 +1961,7 @@ void gen_op::SetLimits(int limit) const
   std::cout << "\n\tLimiting No. Of ";
   if(OpName.length()) std::cout << OpName << " ";
   std::cout << "Operator Representations";
-  std::cout << "\n\tCurrent Representations: " 
-            << size();
+  std::cout << "\n\tCurrent Representations: " << size();
   std::cout << ", Maximum Allowed: " << limit; 
 # endif
   gen_op* OpTmp = const_cast<gen_op*>(this);	// Cast away const (mutable!)

src/HSLib/GenOp.h

@@ -429,6 +429,24 @@ MSVCDLL gen_op exp(const complex& t, double cutoff=1.e-12) const;
         // Note                       : Value of t is considered 0 if
         //                              it's magnituded is less than cutoff
 
+MSVCDLL gen_op expm() const;
+
+	// Input		Op1  : General operator (this)
+        // Return		Op   : exponential of Op1
+	//			       Op = exp(Op1) using pade
+        // Note			     : Computed in the same base as Op1
+
+
+MSVCDLL gen_op expm(const complex& t, double cutoff=1.e-12) const;
+
+        // Input                Op    : Operator (this)
+        //                      t     : Exponential factor
+        //                      cutoff: Exponential factor roundoff
+        // Return               ExpOp : Exponential of Op
+        //                              ExpOp = exp(t*Op) using Pade
+        // Note                       : Exponential output in same base as Op
+        // Note                       : Value of t is considered 0 if
+        //                              it's magnituded is less than cutoff
 
 MSVCDLL gen_op Pow(int power) const;
 

src/LSLib/SuperOp.cc

@@ -1327,6 +1327,124 @@ super_op exp(const super_op& LOp1, const complex& t)
   return LOp;
   }
 
+// ----------------------------------------------------------------------------
+//                 Exponential Superoperators using Pade Method
+// ----------------------------------------------------------------------------
+
+
+	// Input		LOp   : Superoperator (this)
+        // Return		ExpLOp: Exponential of LOp
+	//				ExpLOp = exp(LOp) (Pade method)
+        // Note			      : Exponential output in same base as LOp
+
+super_op super_op::expm() const
+  {
+  if(!HSp)				// Check for NULL LOp
+    {
+    LOperror(5, "exp", 1);		// Error during LOp exp function
+    LOpfatal(7);			// Accessing Null LOp
+    }
+  super_op ExpLOp(*this);		// Copy LOp in its current base
+  matrix mx1 = this->mx;
+
+  ExpLOp.mx = mx1.expm();
+
+  return ExpLOp;
+  }
+
+
+super_op super_op::expm(const complex& t, double cutoff) const
+
+        // Input                LOp   : Superoperator (this)
+        //                      t     : Exponential factor
+        //                      cutoff: Exponential factor roundoff
+        // Return               ExpLOp: Exponential of LOp
+        //                              ExpLOp = exp(t*LOp)
+        // Note                       : Exponential output in EBR of LOp
+        // Note                       : L0p's EBR is generated herein
+        // Note                       : Value of t is considered 0 if
+        //                              it's magnituded is less than cutoff
+
+  {
+  if(!HSp)				// Check for NULL LOp
+    {
+    LOperror(5, "exp", 1);		// Error during LOp exp function
+    LOpfatal(7);			// Accessing NULL Lop
+    }
+  super_op ExpLOp(*this);		// Copy LOp in its current base
+  if(norm(t) < fabs(cutoff))
+    {
+    matrix mx(LSp, LSp, i_matrix_type);	// ExpLOp is the identity superop
+    ExpLOp.mx = mx;
+    }
+  else
+    {
+    matrix mx1 = this->mx*t;
+    ExpLOp.mx = mx1.expm();
+    }
+  return ExpLOp;
+  }
+
+
+
+
+	// Input		LOp1  : Superoperator
+        // Return		LOp   : Exponential of LOp1
+	//				LOp = exp(LOp1)
+        // Note			      : Computed in EBR of LOp1
+        // Note			      : Superoperator output in EBR of LOp1
+
+super_op expm(const super_op& LOp1)
+  {
+  super_op LOp = LOp1;
+  if(LOp1.HSp)				// Check for NULL LOp
+    {
+    matrix m1 = LOp.get_mx();
+    LOp.put_mx(m1.expm());
+    return LOp;
+    }
+  else
+    {
+    LOp1.LOperror(5, "exp", 1);	// Error during LOp exp function
+    LOp1.LOpfatal(7);		// Accessing Null LOp
+    }
+  return LOp;
+  }
+
+
+super_op expm(const super_op& LOp1, const complex& t)
+//return LOp(matrix(LOp1.LSp, LOp1.LSp, d_matrix_type));
+
+	// Input		LOp1  : Superoperator
+	//			t     : A time
+        // Return		LOp   : Exponential of LOp1
+	//				LOp = exp(LOp1*t)
+        // Note			      : Computed in EBR of LOp1
+
+  {
+  super_op LOp = LOp1;
+  if(LOp.HSp) 				// Check for NULL LOp
+    {
+    if(t != complex0)
+      {
+      matrix m1 = LOp.get_mx()*t;
+      LOp.put_mx(m1.expm());
+      return LOp;
+      }
+    else
+      {
+      matrix mx(LOp1.LSp, LOp1.LSp, i_matrix_type);	// LOp is the identity superoperator
+      LOp.mx = mx;
+      }
+    }
+  else
+    {
+    LOp1.LOperror(5, "exp", 1);	// Error during LOp exp function
+    LOp1.LOpfatal(7);		// Accessing Null LOp
+    }
+  return LOp;
+  }
+
 
 // ----------------------------------------------------------------------------
 //                      Superoperators Taken To A Power

src/LSLib/SuperOp.h

@@ -398,6 +398,41 @@ MSVCDLL friend super_op exp(const super_op& LOp1, const complex& t);
 	//				LOp = exp(LOp1*t)
         // Note			      : Computed in EBR of LOp1
 
+MSVCDLL super_op expm() const;
+
+        // Input                LOp   : Superoperator (this)
+        // Return               ExpLOp: Exponential of LOp
+        //                              ExpLOp = exp(LOp) (Pade method)
+        // Note                       : Exponential output in same base as LOp
+
+
+MSVCDLL super_op expm(const complex& t, double cutoff=1.e-12) const;
+
+        // Input                LOp   : Superoperator (this)
+	//			t     : Exponential factor
+	//			cutoff: Exponential factor roundoff
+        // Return               ExpLOp: Exponential of LOp
+        //                              ExpLOp = exp(t*LOp) (Pade method)
+        // Note                       : Exponential output in same base as LOp
+	// Note			      : Value of t is considered 0 if
+	//				it's magnituded is less than cutoff
+
+
+MSVCDLL friend super_op expm(const super_op& LOp1);
+
+	// Input		LOp1  : Superoperator
+        // Return		LOp   : Exponential of LOp1
+	//				LOp = exp(LOp1) (Pade method)
+        // Note			      : Computed in same base as LOp1
+
+
+MSVCDLL friend super_op expm(const super_op& LOp1, const complex& t);
+
+	// Input		LOp1  : Superoperator
+        // Return		LOp   : Exponential of LOp1
+	//				LOp = exp(LOp1*t) (Pade method)
+        // Note			      : Computed in same base as LOp1
+
 
 MSVCDLL friend super_op pow(const super_op& LOp, int power);