function shifted_block_gcrodr_examples
    load('conf5_0-4x4-10.mat');
    A = Problem.A; clear Problem;
    kc = 0.2061;
    A = -A + (1/kc + .001)*speye(size(A)); %scale and shift so that all positive shifts produce physically realistic problems according to matrix documentaiton
    n = size(A,1);
    B = rand(n,5);
    shifts = [0 .0001 .0002 .01 .02];
    opts.isComplex = true;
    
    s = RandStream('mt19937ar','Seed',0);
    RandStream.setGlobalStream(s);

    testShortTimeInterval = 1;
    testMediumTimeInterval = 3;
    testLongTimeInterval = 6;
    
    %% PREAMBLE
    fprintf('Thank you for downloading our shifted block GCRO-DR code!\n\n');
    pause(testMediumTimeInterval);
    fprintf('The following is a series of example runs of our code\nto demonstrate how one can use our codes.\n\n');
    pause(testMediumTimeInterval);
    fprintf('Please cite this code and the associated paper\nif you use them in your work.\n\n');
    pause(testMediumTimeInterval);
    fprintf('All examples use the conf5_0-4x4-10 matrix (scaled and shifted) from\nthe QCD collection in the Univ. Florida Sparse Matrix library.\nFive shifts are used (except in the last experiment) and right-hand sides are\ngenerated randomly.\n\n');
    pause(testMediumTimeInterval);
    fprintf('>> shifts = [0 .0001 .0002 .01 .02];\n\n');
    pause(testLongTimeInterval);clc
    
    %% EXAMPLE 1
    fprintf('Example 1: One option given.\nMethod at a minimum requires only A, shifts, and right-hand sides.\nAll other parameters have defaults\n\n');pause(testMediumTimeInterval);
    fprintf('>> opts.isComplex=true;%%This option set because QCD matrices are complex\n\n');pause(testShortTimeInterval);
    fprintf('>> [X,out]=rgmres_bl_sh(A,B,shifts,opts)\n\n');pause(testShortTimeInterval);
    opts.isOutU = true;
    [~,out]=rgmres_bl_sh(A,B,shifts,opts);
    semilogy(out.resvecFro/out.resvecFro(1)); set(gca,'fontsize',18); xlabel('Number of iterations','interpreter','latex'); ylabel('Relative frobenius residual norm','interpreter','latex');
    pause(testLongTimeInterval);clc; close all;
    
    %% EXAMPLE 2
    fprintf('Example 2: Rerun first experiment using recycled subspace from first run to accelerate\nconvergence.\nAll other parameters have defaults\n\n');pause(testMediumTimeInterval);
    fprintf('>> opts.isComplex=true;%%One option set because QCD matrices are complex\n\n');pause(testShortTimeInterval);
    fprintf('>> opts.isOutU=true;%%Output subspace to recycle\n\n');pause(testShortTimeInterval);
    fprintf('>> [X,out]=rgmres_bl_sh(A,B,shifts,opts)\n\n');pause(testShortTimeInterval);
    fprintf('>> opts.U = out.U;\n\n');pause(testShortTimeInterval);
    fprintf('>> [X,out]=rgmres_bl_sh(A,B,shifts,opts)\n\n');pause(testShortTimeInterval);
    opts.U = out.U;
    [~,out2]=rgmres_bl_sh(A,B,shifts,opts);
    semilogy(out.resvecFro/out.resvecFro(1),'k'); hold on; semilogy(out2.resvecFro/out2.resvecFro(1),'r'); set(gca,'fontsize',18); xlabel('Number of iterations','interpreter','latex'); ylabel('Relative frobenius residual norm','interpreter','latex');
    legend('First run, no initial recycling','Second run, recycling from first run')
    pause(testLongTimeInterval);clc; close all;
    
    %% EXAMPLE 3
    fprintf('Example 3: Options\nAll other options can be set in an options structure.\n');
    fprintf('See documentation for full set of options.\n');
    fprintf('Here we set recycle space size and cycle length.\n\n');pause(testMediumTimeInterval);
    fprintf('>> opts.m = 100; %%cycle length\n');pause(testShortTimeInterval);
    fprintf('>> opts.k = 20; %%recycled subspace dimension\n');pause(testShortTimeInterval);
    fprintf('>> opts.isComplex=true;%%QCD matrices are complex\n');pause(testShortTimeInterval);
    fprintf('>> [X,out]=rgmres_bl_sh(A,B,shifts,opts)\n\n');pause(testShortTimeInterval);
    opts.m = 100;
    opts.k = 20;
    [~,out]=rgmres_bl_sh(A,B,shifts,opts);
    semilogy(out.resvecFro/out.resvecFro(1)); set(gca,'fontsize',18); xlabel('Number of iterations','interpreter','latex'); ylabel('Relative Frobenius residual norm','interpreter','latex');
    pause(testLongTimeInterval);clc; close all; clear opts;
    
    %% EXAMPLE 4
    fprintf('Example 4: Frobenius and 2-norm residuals available.\n');
    fprintf('These results come from Experiment 2.\n');
    fprintf('Located in the out structure.\nSee documenation for full output list.\n\n');pause(testMediumTimeInterval);
    fprintf('>> opts.m = 100; %%cycle length\n');pause(testShortTimeInterval);
    fprintf('>> opts.k = 20; %%recycled subspace dimension\n');pause(testShortTimeInterval);
    fprintf('>> opts.isComplex=true;%%QCD matrices are complex\n');pause(testShortTimeInterval);
    fprintf('>> [X,out]=rgmres_bl_sh(A,B,shifts,opts)\n\n');pause(testShortTimeInterval);
    pause(testMediumTimeInterval);
    fprintf('>> semilogy(out.resvecFro/out.resvecFro(1), ''k''); hold on;\n>>semilogy(out.resvec2norm(:,1)/out.resvec2norm(1,1), ''r-'');\n>>semilogy(out.resvec2norm(:,2)/out.resvec2norm(1,2), ''r-.'');\n>>semilogy(out.resvec2norm(:,3)/out.resvec2norm(1,3), ''r--'');\n>>semilogy(out.resvec2norm(:,4)/out.resvec2norm(1,4), ''k-.'');\n>>semilogy(out.resvec2norm(:,5)/out.resvec2norm(1,5), ''k--'');\n\n');pause(testMediumTimeInterval);
    semilogy(out.resvecFro/out.resvecFro(1), 'k'); hold on; semilogy(out.resvec2norm(:,1)/out.resvec2norm(1,1), 'r-'); semilogy(out.resvec2norm(:,2)/out.resvec2norm(1,2), 'r-.'); semilogy(out.resvec2norm(:,3)/out.resvec2norm(1,3), 'r-.'); semilogy(out.resvec2norm(:,4)/out.resvec2norm(1,4), 'k--'); semilogy(out.resvec2norm(:,5)/out.resvec2norm(1,5), 'k-.');  hold off;   set(gca,'fontsize',18); xlabel('Number of iterations','interpreter','latex'); ylabel('Relative residual norm','interpreter','latex');
    pause(testLongTimeInterval);clc; close all; clear opts;
    
    %% EXAMPLE 5
    fprintf('Example 5: Recycling with Ritz versus harmonic Ritz vectors.\n');
    fprintf('We compare performance for the two avaible recycling routines.\n''ritz'' default computed in last example.\n\n');pause(testMediumTimeInterval);
    fprintf('>> opts.m = 100; %%cycle length\n');pause(testShortTimeInterval);
    fprintf('>> opts.k = 20; %%recycled subspace dimension\n');pause(testShortTimeInterval);
    fprintf('>> opts.isComplex=true;%%QCD matrices are complex\n');pause(testShortTimeInterval);
    fprintf('>> opts.recycleMethod = ''harmritz''; \n');pause(testShortTimeInterval);
    fprintf('>> [X,out]=rgmres_bl_sh(A,B,shifts,opts)\n\n');pause(testShortTimeInterval);
    opts.recycleMethod = 'harmritz';
    opts.m = 100;
    opts.k = 20;
    [~,out2]=rgmres_bl_sh(A,B,shifts,opts);
    semilogy(out.resvecFro, 'k'); hold on; semilogy(out2.resvecFro, 'r'); hold off; set(gca,'fontsize',18); xlabel('Number of iterations','interpreter','latex'); ylabel('Relative residual norm','interpreter','latex');
    legend('Ritz recycling','harmonic Ritz recycling');
    pause(testLongTimeInterval);clc; close all; clear opts;
    
    %% EXAMPLE 6
    fprintf('Example 6: Continuous dependence of the right-hand side on the shift.\n');
    fprintf('If right-hand side depends sufficiently continuously on the shift\n(in some intervale of shifts), solutions for shifts in this interval span a ''small'' subspace.\nWe use this to solve 100 shifted systems by solving the first 10 and using those\nsolutions as a recycling subspace.\n\n');pause(testMediumTimeInterval);
    fprintf('>> shifts = rand(1, 100)+1;%%100 random positive shifts\n');pause(testShortTimeInterval);
    fprintf('>> B = sin(kron(shifts,(1:n)''*pi/n));%%columns of B depend continuously on the shifts\n');pause(testShortTimeInterval);
    fprintf('>> opts.m = 100;\n');pause(testShortTimeInterval);
    fprintf('>> opts.k = 20;\n\n');pause(testShortTimeInterval);
    fprintf('>> opts.isComplex=true;%%QCD matrices are complex;\n\n');pause(testShortTimeInterval);
    fprintf('>> [X(:,1:10),out1{1}]=rgmres_bl_sh(A,B(:,1:10),shifts(1:10),opts);%%solve for first ten shifts\n\n');pause(testShortTimeInterval);
    shifts = rand(1, 100)+1;
    B = sin(kron(shifts,(1:n)'*pi/n));
    X = zeros(n,100);
    opts.m = 100;
    opts.k = 20;
    [X(:,1:10),out1{1}]=rgmres_bl_sh(A,B(:,1:10),shifts(1:10),opts);
    pause(testMediumTimeInterval);
    fprintf('>> opts.U = X(:,1:10);%%Explicitely include X(:,1:10) as augmenting subspace\n');pause(testShortTimeInterval);
    fprintf('>> for i=2:10\n        [X(:,(i-1)*10+1:i*10),out1{i}]=rgmres_bl_sh(A,B(:,(i-1)*10+1:i*10),shifts((i-1)*10+1:i*10),opts);%%solve for first ten shifts\nend\n\n');pause(testShortTimeInterval);
    opts.U = X(:,1:10);
    for i=2:10
        [X(:,(i-1)*10+1:i*10),out1{i}]=rgmres_bl_sh(A,B(:,(i-1)*10+1:i*10),shifts((i-1)*10+1:i*10),opts);
    end
    semilogy(out1{1}.resvec2norm,'r--');  hold on;
    for i=2:10
        semilogy(size(out1{1}.resvec2norm,1) + (1:size(out1{i}.resvec2norm,1)),out1{i}.resvec2norm,'k--'); 
    end
    hold off;
    set(gca,'fontsize',18); xlabel('Number of iterations','interpreter','latex'); ylabel('Relative residual norm','interpreter','latex');
    title(sprintf('\\textbf{red}: first ten convergence curves;\n\\textbf{black}: other 90 convergence curves with augmentation'), 'interpreter', 'latex');
    pause(testLongTimeInterval);clc; close all; clear opts;
    
    %% CITATION
    fprintf('Kirk M. Soodhalter\n"Block Krylov Subspace Recycling for Shifted\nSystems with Unrelated Right-Hand Sides"\nSISC 38(1), A302-A324, 2016\n\n')
end