Update arnoldiqr.m

arnoldiqr.m

@@ -1,26 +1,45 @@
-
-function [Q, R, Qn] = arnoldiqr(diag, E, q1, m) % initialization
-    d = length(diag);
-    tau = size(E,1);
+function [Q, R, Qn] = arnoldiqr(D, E, q1, m, a, precond) % initialization
 
+    d = length(D);
+    tau = size(E,1);
     n = d + tau;
 
     if m <= 0
         display("Error. m must be > 0");
-    return;
+        return;
     elseif m > n
         m = n;
         display("You don't need more than n iterations. m has been set to n");
     end
+
     q1 = q1/norm(q1);
     Qn = zeros(n, m+1); 
     Qn(:, 1) = q1;
     Hn1 = zeros(m,1); %diagonal of H_n
     Hn2 = zeros(m,1); %subdiagonal of H_n
     Q = zeros(m+1, m+1);
     R = zeros(m+1, m);
+
+    if precond
+        D3 = D.^3; 
+    end
+
     for k = 1:m
-        z = [ (diag.*Qn(1:d, k)) + (E'*Qn(d+1:n, k)); E*Qn(1:d, k)]; %we want to avoid full computation since the matrices are sparse
+
+        qk1 = Qn(1:d, k);
+        qk2 = Qn(d+1:n, k);
+
+        if precond %generating the new vector of the Krylov Subspace (precond or not)
+            z1 = (D3.*qk1) + E'*(E*(D.*qk1)) + D.*(E'*(E*qk1)) -D.*(E'*(E*((E'*qk2)./D)));
+            z2 = -E*((E'*(E*(D.*qk1)))./D);
+            z = [z1; z2];
+
+            z = [z(1:d) + a*(D.*(E'*qk2)); z(d+1:end)+ a*(E*(D.*qk1))];
+        else
+            z = [ (D.*qk1) + (E'*qk2); E*qk1]; 
+        end
+
+        %COMPUTING NEW COMPONENTS OF Hn
         if k==1
             Hn1(1) = Qn(:,1)'*z;
             z = z - Hn1(1)*Qn(:, 1);
@@ -31,11 +50,15 @@
             z = z - Hn1(k)*Qn(:, k);
             Hn2(k) = norm(z);
         end
-        if abs(Hn2(k)) < 1e-10 %breakdown%
-            display("BREAKDOWN HAPPENED!");
+
+        if abs(Hn2(k)) < 1e-10 % check for breakdown
+            disp("BREAKDOWN HAPPENED!");
             return;
         end
-        Qn(:, k+1) = z/Hn2(k);
+
+        Qn(:, k+1) = z/Hn2(k); %new element of the orthonormal basis of Krylov Subspace
+
+        %QR DECOMPOSITION OF Hn (1st AND OTHER ITERATIONS)
         if k == 1
             x = [Hn1(1); Hn2(1)];
             y = norm(x)*[1; 0];
@@ -54,13 +77,14 @@
             v = x - y;
             v_norm = norm(v)^2;
             R(1:k, k) = [c(1:k-1); x_norm];
+
             %Building Q tilde%
-            Q(k+1,k+1)=1; %build Q' signed%
+            Q(k+1,k+1)=1; %build \Bar{Q'}
             u = Q(1:k-1, k);
-            a = Q(k, k);
-            p = [u zeros(k-1, 1); a 0; 0 1]; %last two rows of Q signed transponed
-            gamma = [v(1)^2*u v(1)*v(2)*u; v(1)^2*a v(1)*v(2)*a; v(1)*v(2) v(2)^2];
+            a2 = Q(k, k);
+            p = [u zeros(k-1, 1); a2 0; 0 1]; %last two rows of Q signed transponed
+            gamma = [v(1)^2*u v(1)*v(2)*u; v(1)^2*a2 v(1)*v(2)*a2; v(1)*v(2) v(2)^2];
             Q(1:k+1, k:k+1) = p-(2/(v_norm))*gamma;
         end
     end
-
+end