added expm

miscmaths.cc

@@ -1408,6 +1408,73 @@ ReturnMatrix exp(const Matrix& mat)
   return res;
 }
 
+  // optimised code for calculating matrix exponential
+ReturnMatrix expm(const Matrix& mat){
+  float nmat = sum(mat).Maximum();
+  int nc=mat.Ncols(),nr=mat.Nrows();
+  Matrix res(nr,nc);
+  IdentityMatrix id(nr);
+  Matrix U(nr,nc),V(nr,nc);
+
+  if(nmat <= 1.495585217958292e-002){ // m=3
+    Matrix mat2(nr,nc);
+    mat2=mat*mat;
+    U = mat*(mat2+60.0*id);
+    V = 12.0*mat2+120.0*id;
+    res = (-U+V).i()*(U+V);
+  }
+  else if(nmat <= 2.539398330063230e-001){ // m=5
+    Matrix mat2(nr,nc),mat4(nr,nc);
+    mat2=mat*mat;mat4=mat2*mat2;
+    U = mat*(mat4+420.0*mat2+15120.0*id);
+    V = 30.0*mat4+3360.0*mat2+30240.0*id;
+    res = (-U+V).i()*(U+V);
+  }
+  else if(nmat <= 9.504178996162932e-001){ // m=7
+    Matrix mat2(nr,nc),mat4(nr,nc),mat6(nr,nc);
+    mat2=mat*mat;mat4=mat2*mat2,mat6=mat4*mat2;
+    U = mat*(mat6+1512.0*mat4+277200.0*mat2+8648640.0*id);
+    V = 56.0*mat6+25200.0*mat4+1995840.0*mat2+17297280.0*id;
+    res = (-U+V).i()*(U+V);
+  }
+  else if(nmat <= 2.097847961257068e+000){
+    Matrix mat2(nr,nc),mat4(nr,nc),mat6(nr,nc),mat8(nr,nc);
+    mat2=mat*mat;mat4=mat2*mat2,mat6=mat4*mat2,mat8=mat6*mat2;
+    U = mat*(mat8+3960.0*mat6+2162160.0*mat4+302702400.0*mat2+8821612800.0*id);
+    V = 90.0*mat8+110880.0*mat6+30270240.0*mat4+2075673600.0*mat2+17643225600.0*id;
+    res = (-U+V).i()*(U+V);
+  }
+  else if(nmat <= 5.371920351148152e+000){
+    Matrix mat2(nr,nc),mat4(nr,nc),mat6(nr,nc);
+    mat2=mat*mat;mat4=mat2*mat2,mat6=mat4*mat2;
+    U = mat*(mat6*(mat6+16380.0*mat4+40840800.0*mat2)+
+	     +33522128640.0*mat6+10559470521600.0*mat4+1187353796428800.0*mat2+32382376266240000.0*id);
+    V = mat6*(182.0*mat6+960960.0*mat4+1323241920.0*mat2)
+      + 670442572800.0*mat6+129060195264000.0*mat4+7771770303897600.0*mat2+64764752532480000.0*id;
+    res = (-U+V).i()*(U+V);
+  }
+  else{
+    double t;int s;
+    t = frexp(nmat/5.371920351148152,&s);
+    if(t==0.5) s--;
+    t = std::pow(2.0,s);
+    res = (mat/t);
+    Matrix mat2(nr,nc),mat4(nr,nc),mat6(nr,nc);
+    mat2=res*res;mat4=mat2*mat2,mat6=mat4*mat2;
+    U = res*(mat6*(mat6+16380*mat4+40840800*mat2)+
+	     +33522128640.0*mat6+10559470521600.0*mat4+1187353796428800.0*mat2+32382376266240000.0*id);
+    V = mat6*(182.0*mat6+960960.0*mat4+1323241920.0*mat2)
+      + 670442572800.0*mat6+129060195264000.0*mat4+7771770303897600.0*mat2+64764752532480000.0*id;
+    res = (-U+V).i()*(U+V);
+    for(int i=1;i<=s;i++)
+      res = res*res;
+  }
+
+  res.Release();
+  return res;
+}
+
+
 ReturnMatrix tanh(const Matrix& mat)
 {
   Matrix res = mat;

miscmaths.h

@@ -217,6 +217,7 @@ namespace MISCMATHS {
   ReturnMatrix sqrtm(const Matrix& mat);
   ReturnMatrix log(const Matrix& mat);
   ReturnMatrix exp(const Matrix& mat);
+  ReturnMatrix expm(const Matrix& mat);
   ReturnMatrix tanh(const Matrix& mat);
   ReturnMatrix pow(const Matrix& mat, const double exp);
   ReturnMatrix sum(const Matrix& mat, const int dim = 1);