Dqgmres

[dpo]

Introduce rotpos and jrotpos indices used to index into c and s. Previously, it was possible that kpos or jpos = mem + 1 but c and s have length mem. Matlab probably silently expanded c and s. One of my students noticed this while implementing standard DQGMRES for another project.

[dpo]

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)?!

[dpo]

I perform the following experiment:

>> n = 15; m = 8; A = rand(n,n); B = rand(m,n); C = eye(m); b = [rand(n,1) ; zeros(m,1)]; K = [A B' ; B -C]; G = eye(n); M = [G B' ; B -C];
>> opts
opts = 
      print: 1
    restart: 20
        mem: 30
>> [xx, flag, relres, iter, resvec] = gmres(K, b, 30, 1.0e-6, 30, M);  % Matlab's own GMRES with CP

>> [x, stats, flag] = reg_cpkrylov(@cpgmres, b, A, B, C, G, opts);
>> % the relative residuals match accurately, except the very last one
>> norm(x - xx) / norm(xx)
ans =
   2.1124e-12

>> [x, stats, flag] = reg_cpkrylov(@cpgmres, b, A, B, C, G, opts);
>> % the relative residuals match accurately, except the very last one
>> norm(x - xx) / norm(xx)
ans =
   2.1127e-12

[diserafi]

The fixes of indexing seem ok. I've not checked yet the inaccuracy in the very last iteration.

[diserafi]

I repeated your experiment several times with larger values of n and m, setting the restart parameter equal to n+m, so that the solution computed at the last iteration should be exact (but it is not because of the roundoff error). The difference between the residual norms computed by cpgmres and Matlab gmres seems to increase as the residual norm decreases (you clearly see it in the last iteration), although the solutions computed by the two gmres versions are very close (the relative error is close to the unit roundoff).

By looking into the gmres implementation provided by Matlab, I found that it uses Householder transformations in the generation of the Krylov basis, in order to reduce the numerical error. Quoting the abstract of H.F. Walker, Implementation of the GMRES Method Using Householder Transformations, SIAM J. Sci. Comp. 9 (1), 1988: "numerical experiments suggest that it [the implementation based on Householder transformations] is more stable, especially as the limits of residual reduction are reached. This could explain what we observed in the experiments. Furthermore, I suspect there are other improvements that can be applied to our codes to reduce the roundoff error.

However, I suggest you to proceed with the revision of the other codes I uploaded on GitHub and make the experiments we discussed about when we met in Bordeaux, so that we can finalize a version of the paper that can be (hopefully) submitted. We can (and certainly will) make further modifications later, taking also into account the comments we will receive from the reviewers. Please, let me know if you are still interested in this work.