Optimize using in-place math; small cleanups

src/expmv_fun.jl

@@ -1,25 +1,20 @@
 export expmv
 
-function expmv(t, A, b; M = [], prec = "double", shift = false, full_term = false, prnt = false)
+function expmv(t, A, b; M::Array{Float64, 2} = Array(Float64, 0, 0), prec = "double", shift = false, full_term = false, prnt = false)
                # bal = false, 
 
     #EXPMV   Matrix exponential times vector or matrix.
-    #   [F,S,M,MV,MVD] = EXPMV(t,A,B,[],PREC) computes EXPM(t*A)*B without
+    #   [F,S,M] = EXPMV(t,A,B,[],PREC) computes EXPM(t*A)*B without
     #   explicitly forming EXPM(t*A). PREC is the required accuracy, 'double',
     #   'single' or 'half', and defaults to CLASS(A).
     #
-    #   A total of MV products with A or A^* are used, of which MVD are
-    #   for norm estimation.
-    #
     #   The full syntax is
     #
-    #     [f,s,m,mv,mvd,unA] = expmv(t,A,b,M,prec,shift,bal,full_term,prnt).
-    #
-    #   unA = 1 if the alpha_p were used instead of norm(A).
+    #     f = expmv(t,A,b,M,prec,shift,bal,full_term,prnt).
     #
     #   If repeated invocation of EXPMV is required for several values of t
     #   or B, it is recommended to provide M as an external parameter as
-    #   M = SELECT_TAYLOR_DEGREE(A,m_max,p_max,prec,shift,bal,true).
+    #   M = SELECT_TAYLOR_DEGREE(A,b_columns,m_max,p_max,prec,shift,bal,true).
     #   This also allows choosing different m_max and p_max.
 
     #   Reference: A. H. Al-Mohy and N. J. Higham, Computing the action of
@@ -41,31 +36,27 @@ function expmv(t, A, b; M = [], prec = "double", shift = false, full_term = fals
     n = size(A,1)
 
     if shift
-        mu = trace(A)/n 
-        #mu = full(mu) # Much slower without the full!
+        mu = trace(A)/n
         A = A-mu*speye(n)
     end
 
     if isempty(M)
-        tt = 1
-        (M,mvd,alpha,unA) = select_taylor_degree(t*A,b)
-        mv = mvd
+        tt = 1.0
+        M = select_taylor_degree(t * A, size(b, 2))
     else
         tt = t
-        mv = 0
-        mvd = 0
     end
 
-    tol =   
-      if prec == "double"
-          2.0^(-53)
+    tol =
+      if prec == "half"
+          2.0^(-10)
       elseif prec == "single"
           2.0^(-24)
-      elseif prec == "half"   
-          2.0^(-10)
+      else
+          2.0^(-53)
       end
 
-    s = 1;
+    s = 1.0;
 
     if t == 0
         m = 0;
@@ -83,15 +74,13 @@ function expmv(t, A, b; M = [], prec = "double", shift = false, full_term = fals
             cost,m = findmin(C);  # when C is one column. Happens if p_max = 2.
         end
         if cost == Inf
-            cost = 0 
+            cost = 0.0
         end
-        s = max(cost/m,1);
+        s = max(cost/m,1.0);
     end
 
-    eta = 1;
-
     if shift 
-        eta = exp(t*mu/s)
+        eta = exp(t * mu / s)
     end
 
     f = b;
@@ -104,7 +93,6 @@ function expmv(t, A, b; M = [], prec = "double", shift = false, full_term = fals
         c1 = norm(b,Inf);
         for k = 1:m
             b = (t/(s*k))*(A*b);
-            mv = mv + 1;
             f =  f + b;
             c2 = norm(b,Inf);
             if !full_term
@@ -118,7 +106,9 @@ function expmv(t, A, b; M = [], prec = "double", shift = false, full_term = fals
             end
 
         end
-        f = eta*f
+        if shift
+            f = eta*f
+        end
         b = f
     end
 
@@ -130,6 +120,5 @@ function expmv(t, A, b; M = [], prec = "double", shift = false, full_term = fals
     #    f = D*f
     #end
 
-    #return (f,s,m,mv,mvd,unA)
     return f
 end

src/normAm.jl

@@ -56,15 +56,13 @@ function normAm(A,m)
             e = A'*e;
         end
         c = norm(e,Inf);
-        mv = m;
     else
         #(c,_,_,it) = normest1(@afun_power,t);
         #mv = it[2]*t*m;
         c = norm(A^m,1);
-        mv = 6*t*m;
     end
 
-    return (c,mv)
+    return c
 
 end
 

src/select_taylor_degree.jl

@@ -1,8 +1,5 @@
-#function  [M,mv,alpha,unA] = ...
-#           select_taylor_degree(A,b,m_max,p_max,prec,shift,bal,force_estm)
-
 function select_taylor_degree(A,
-                              b;
+                              b_columns;
                               m_max = 55,
                               p_max = 8,
                               prec = "double",
@@ -12,10 +9,9 @@ function select_taylor_degree(A,
 
 
     #SELECT_TAYLOR_DEGREE   Select degree of Taylor approximation.
-    #   [M,MV,alpha,unA] = SELECT_TAYLOR_DEGREE(A,m_max,p_max) forms a matrix M
+    #   [M,alpha,unA] = SELECT_TAYLOR_DEGREE(A,m_max,p_max) forms a matrix M
     #   for use in determining the truncated Taylor series degree in EXPMV
     #   and EXPMV_TSPAN, based on parameters m_max and p_max.
-    #   MV is the number of matrix-vector products with A or A^* computed.
 
     #   Reference: A. H. Al-Mohy and N. J. Higham, Computing the action of
     #   the matrix exponential, with an application to exponential
@@ -48,24 +44,19 @@ function select_taylor_degree(A,
         A = A-mu*speye(n);
     end
 
-    mv = 0;
-
     if !force_estm
         normA = norm(A,1)
     end
 
-    if !force_estm && normA <= 4*theta[m_max]*p_max*(p_max + 3)/(m_max*size(b,2));
-        unA = 1;
+    if !force_estm && normA <= 4*theta[m_max]*p_max*(p_max + 3)/(m_max*b_columns);
         c = normA;
         alpha = c*ones(p_max-1,1);
     else
-        unA = 0;
         eta = zeros(p_max,1); 
         alpha = zeros(p_max-1,1);
         for p = 1:p_max
-            (c,k) = normAm(A,p+1);
+            c = normAm(A,p+1);
             c = c^(1/(p+1));
-            mv = mv + k;
             eta[p] = c;
         end
         for p = 1:p_max-1
@@ -81,6 +72,5 @@ function select_taylor_degree(A,
         end
     end
 
-    return (M,mv,alpha,unA)
-
-end
+    return M
+end