Skip to content

Dqgmres #1

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
dpo wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Dqgmres #1

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
2411fae
fix line endings
dpo Sep 14, 2018
203bf8d
fix indexing into c and s
dpo Sep 14, 2018
201c27a
fix indexing
dpo Sep 14, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
38 changes: 20 additions & 18 deletions kernels/cpdqgmres.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -138,18 +138,18 @@
c = zeros(mem, 1); % Givens cosines.
s = zeros(mem, 1); % Givens sines.
H = zeros(itmax, mem+2); % Upper Hessenberg form of A, with mem upper
% diagonals. Its mem+2 diagonals are stored
% diagonals. Its mem+2 diagonals are stored
% as column vectors, according to the tranformation
% (j,k) --> (j,2+k-j). ALL THIS MEMORY IS
% NOT NEEDED. TODO: REDUCE MEMORY FOR H.

% Set up zero vectors
zeron = zeros(n,1);
zerom = zeros(m,1);

% Initialize some vectors
x = zeron;
y = zerom; % yk = y0 - qk, y0 = 0 ==> yk = -qk
y = zerom; % yk = y0 - qk, y0 = 0 ==> yk = -qk
u = b; % u_0 = b - A * x_0 = b
t = zerom; % t_0 = C * q_0 = 0

Expand Down Expand Up @@ -179,18 +179,19 @@
end

% Main loop.

% The following stopping criterion compensates for the lag in the
% residual, but usually increases the number of iterations.
% residual, but usually increases the number of iterations.
% while sqrt(max(1, k-mem+1)) * residNorm > stopTol && k < itmax

while residNorm > stopTol && k < itmax % less accurate, but acceptable

k = k + 1;

% Set position in circular stack where (k+1)-st Krylov vector should go.
kpos = mod(k-1, mem+1) + 1; % Position corresponding to k in the circular stack.
kp1pos = mod(k, mem+1) + 1; % Position corresponding to k+1 in the circular stack.
rotpos = mod(k-1, mem) + 1; % Position of the current rotation parameters

% Compute next Krylov vectors from the modified Gram-Schmidt process.
% Only orthogonalize against the most recent min(k,mem) vectors.
Expand All @@ -201,40 +202,41 @@
Q(:,kp1pos) = Q(:,kpos) - w(n+1:n+m,1);
for j = max(1, k-mem+1) : k
jpos = mod(j-1, mem+1) + 1;
kk = 2+k-j;
kk = 2+k-j;
Copy link
Copy Markdown
Member Author

@dpo dpo Sep 14, 2018

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There might be something else going on here or at line 211 below. At the very first pass through the main loop, k=1, j=1 in this loop, so that kk = 2. We write to H(1,2) and on line 211, we write to H(2,1). We never write anything in H(1,1)?!

H(j,kk) = dot(V(:,jpos), u) + dot(Q(:,jpos), t);
V(:,kp1pos) = V(:,kp1pos) - H(j,kk) * V(:,jpos);
Q(:,kp1pos) = Q(:,kp1pos) - H(j,kk) * Q(:,jpos);
end
% kk = k-(k+1)+2 = 1
H(k+1,1) = sqrt(dot(u, V(:,kp1pos)) + dot(t, Q(:,kp1pos)));
H(k,1) = sqrt(dot(u, V(:,kp1pos)) + dot(t, Q(:,kp1pos)));

if H(k+1,1) ~= 0 % Lucky breakdown if = 0.
V(:,kp1pos) = V(:,kp1pos) / H(k+1,1);
Q(:,kp1pos) = Q(:,kp1pos) / H(k+1,1);
if H(k,1) ~= 0 % Lucky breakdown if = 0.
V(:,kp1pos) = V(:,kp1pos) / H(k,1);
Q(:,kp1pos) = Q(:,kp1pos) / H(k,1);
end

% Apply previous (symmetric) Givens rotations.
for j = max(1,k-mem) : k-1
jpos = mod(j-1, mem+1) + 1;
jp1pos = mod(j, mem+1) + 1;
jrotpos = mod(j-1, mem) + 1;
kk = k-j+1; % kk = 2+k-(j+1)
kk1 = kk+1; % kk1 = 2+k-j
Hjk = c(jpos) * H(j,kk1) + s(jpos) * H(j+1,kk);
H(j+1,kk) = s(jpos) * H(j,kk1) - c(jpos) * H(j+1,kk);
Hjk = c(jrotpos) * H(j,kk1) + s(jrotpos) * H(j+1,kk);
H(j+1,kk) = s(jrotpos) * H(j,kk1) - c(jrotpos) * H(j+1,kk);
H(j,kk1) = Hjk;
end

% Compute and apply current (symmetric) Givens rotation:
% [ck sk] [H(k,k) ] = [*]
% [sk -ck] [H(k+1,k)] [0].
% Indices for H:
% (k+1,k) --> (k+1,2+k-(k+1)) = (k+1,1)
% (k,k) --> (k,2+k-k) = (k,2)
[c(kpos), s(kpos), H(k,2)] = SymGivens(H(k,2), H(k+1,1));
H(k+1,1) = 0;
g(kp1pos) = s(kpos) * g(kpos);
g(kpos) = c(kpos) * g(kpos);
[c(rotpos), s(rotpos), H(k,2)] = SymGivens(H(k,2), H(k,1));
H(k,1) = 0;
g(kp1pos) = s(rotpos) * g(kpos);
g(kpos) = c(rotpos) * g(kpos);

% Update directions PV and PQ, and solution [x; y].
PV(:,kpos) = V(:,kpos);
Expand Down