% Diary from Matlab demos in CS 290H, Monday, April 28, 2008


clc
A = grid5(10);
n = length(A)
n =
   100
spy(A),shg
b = rand(n,1);
help pcg
 PCG   Preconditioned Conjugate Gradients Method.
    X = PCG(A,B) attempts to solve the system of linear equations A*X=B for
    X. The N-by-N coefficient matrix A must be symmetric and positive
    definite and the right hand side column vector B must have length N.
 
    X = PCG(AFUN,B) accepts a function handle AFUN instead of the matrix A.
    AFUN(X) accepts a vector input X and returns the matrix-vector product
    A*X. In all of the following syntaxes, you can replace A by AFUN.
 
    X = PCG(A,B,TOL) specifies the tolerance of the method. If TOL is []
    then PCG uses the default, 1e-6.
 
    X = PCG(A,B,TOL,MAXIT) specifies the maximum number of iterations. If
    MAXIT is [] then PCG uses the default, min(N,20).
 
    X = PCG(A,B,TOL,MAXIT,M) and X = PCG(A,B,TOL,MAXIT,M1,M2) use symmetric
    positive definite preconditioner M or M=M1*M2 and effectively solve the
    system inv(M)*A*X = inv(M)*B for X. If M is [] then a preconditioner
    is not applied. M may be a function handle MFUN returning M\X.
 
    X = PCG(A,B,TOL,MAXIT,M1,M2,X0) specifies the initial guess. If X0 is
    [] then PCG uses the default, an all zero vector.
 
    [X,FLAG] = PCG(A,B,...) also returns a convergence FLAG:
     0 PCG converged to the desired tolerance TOL within MAXIT iterations
     1 PCG iterated MAXIT times but did not converge.
     2 preconditioner M was ill-conditioned.
     3 PCG stagnated (two consecutive iterates were the same).
     4 one of the scalar quantities calculated during PCG became too
       small or too large to continue computing.
 
    [X,FLAG,RELRES] = PCG(A,B,...) also returns the relative residual
    NORM(B-A*X)/NORM(B). If FLAG is 0, then RELRES <= TOL.
 
    [X,FLAG,RELRES,ITER] = PCG(A,B,...) also returns the iteration number
    at which X was computed: 0 <= ITER <= MAXIT.
 
    [X,FLAG,RELRES,ITER,RESVEC] = PCG(A,B,...) also returns a vector of the
    residual norms at each iteration including NORM(B-A*X0).
 
    Example:
       n1 = 21; A = gallery('moler',n1);  b1 = A*ones(n1,1);
       tol = 1e-6;  maxit = 15;  M = diag([10:-1:1 1 1:10]);
       [x1,flag1,rr1,iter1,rv1] = pcg(A,b1,tol,maxit,M);
    Or use this parameterized matrix-vector product function:
       afun = @(x,n)gallery('moler',n)*x;
       n2 = 21; b2 = afun(ones(n2,1),n2);
       [x2,flag2,rr2,iter2,rv2] = pcg(@(x)afun(x,n2),b2,tol,maxit,M);
 
    Class support for inputs A,B,M1,M2,X0 and the output of AFUN:
       float: double
 
    See also <a href="matlab:help bicg">bicg</a>, <a href="matlab:help bicgstab">bicgstab</a>, <a href="matlab:help cgs">cgs</a>, <a href="matlab:help gmres">gmres</a>, <a href="matlab:help lsqr">lsqr</a>, <a href="matlab:help minres">minres</a>, <a href="matlab:help qmr">qmr</a>, <a href="matlab:help symmlq">symmlq</a>,
    <a href="matlab:help cholinc">cholinc</a>, <a href="matlab:help function_handle">function_handle</a>.


    Reference page in Help browser
       <a href="matlab:doc pcg">doc pcg</a>



[x,flag,relres,iter,resvec] = pcg(A,b,0,n);



size(x)
ans =
   100     1
norm(b - A*x)
ans =
  3.3433e-014
size(resvec)
ans =
    46     1
resvec
resvec = (omitted)
   
semilogy(resvec / resvec(1),'k'),shg
[x,flag,relres,iter,resvec] = pcg(A,b,1e-8,n);
semilogy(resvec / resvec(1),'k'),shg
A = grid5(100);
b = rand(length(A),1);
[x,flag,relres,iter,resvec] = pcg(A,b,1e-8,n);
hold on
semilogy(resvec / resvec(1),'r'),shg
n = length(A)
n =
       10000
[x,flag,relres,iter,resvec] = pcg(A,b,1e-8,n);
semilogy(resvec / resvec(1),'r'),shg
A = grid5(20);
n = length(A)
n =
   400
b = rand(length(A),1);
[x,flag,relres,iter,resvec] = pcg(A,b,1e-8,n);
semilogy(resvec / resvec(1),'g'),shg

A = BCSSTK08;  % This matrix comes from from CGmats.mat
spy(A),shg
clf
spy(A),shg
n = length(A)
n =
        1074
b = rand(length(A),1);
[x,flag,relres,iter,resvec] = pcg(A,b,1e-8,n);
semilogy(resvec / resvec(1),'g'),shg
[x,flag,relres,iter,resvec] = pcg(A,b,1e-8,10*n);
semilogy(resvec / resvec(1),'k'),shg
clc

s = (1:100 )/ 100;
s
s =
  Columns 1 through 9 
    0.0100    0.0200    0.0300    0.0400    0.0500    0.0600    0.0700    0.0800    0.0900
  Columns 10 through 18 
    0.1000    0.1100    0.1200    0.1300    0.1400    0.1500    0.1600    0.1700    0.1800
  Columns 19 through 27 
    0.1900    0.2000    0.2100    0.2200    0.2300    0.2400    0.2500    0.2600    0.2700
  Columns 28 through 36 
    0.2800    0.2900    0.3000    0.3100    0.3200    0.3300    0.3400    0.3500    0.3600
  Columns 37 through 45 
    0.3700    0.3800    0.3900    0.4000    0.4100    0.4200    0.4300    0.4400    0.4500
  Columns 46 through 54 
    0.4600    0.4700    0.4800    0.4900    0.5000    0.5100    0.5200    0.5300    0.5400
  Columns 55 through 63 
    0.5500    0.5600    0.5700    0.5800    0.5900    0.6000    0.6100    0.6200    0.6300
  Columns 64 through 72 
    0.6400    0.6500    0.6600    0.6700    0.6800    0.6900    0.7000    0.7100    0.7200
  Columns 73 through 81 
    0.7300    0.7400    0.7500    0.7600    0.7700    0.7800    0.7900    0.8000    0.8100
  Columns 82 through 90 
    0.8200    0.8300    0.8400    0.8500    0.8600    0.8700    0.8800    0.8900    0.9000
  Columns 91 through 99 
    0.9100    0.9200    0.9300    0.9400    0.9500    0.9600    0.9700    0.9800    0.9900
  Column 100 
    1.0000
A = randspd(s);
clc
sort(eig(A))
ans =
    0.0100
    0.0200
    0.0300
    0.0400
    0.0500
    0.0600
    0.0700
    0.0800
    0.0900
    0.1000
    0.1100
    0.1200
    0.1300
    0.1400
    0.1500
    0.1600
    0.1700
    0.1800
    0.1900
    0.2000
    0.2100
    0.2200
    0.2300
    0.2400
    0.2500
    0.2600
    0.2700
    0.2800
    0.2900
    0.3000
    0.3100
    0.3200
    0.3300
    0.3400
    0.3500
    0.3600
    0.3700
    0.3800
    0.3900
    0.4000
    0.4100
    0.4200
    0.4300
    0.4400
    0.4500
    0.4600
    0.4700
    0.4800
    0.4900
    0.5000
    0.5100
    0.5200
    0.5300
    0.5400
    0.5500
    0.5600
    0.5700
    0.5800
    0.5900
    0.6000
    0.6100
    0.6200
    0.6300
    0.6400
    0.6500
    0.6600
    0.6700
    0.6800
    0.6900
    0.7000
    0.7100
    0.7200
    0.7300
    0.7400
    0.7500
    0.7600
    0.7700
    0.7800
    0.7900
    0.8000
    0.8100
    0.8200
    0.8300
    0.8400
    0.8500
    0.8600
    0.8700
    0.8800
    0.8900
    0.9000
    0.9100
    0.9200
    0.9300
    0.9400
    0.9500
    0.9600
    0.9700
    0.9800
    0.9900
    1.0000

edit randspd
size(A)
ans =
   100   100
nnz(A)
ans =
       10000
n = length(A)
n =
   100
b = rand(length(A),1);
[x,flag,relres,iter,resvec] = pcg(A,b,1e-8,n);
close all
semilogy(resvec / resvec(1),'k'),shg
s = (s+1)/2;
min(s)
ans =
    0.5050
max(s)
ans =
     1
A = randspd(s);
[x,flag,relres,iter,resvec] = pcg(A,b,1e-8,n);
hold on
semilogy(resvec / resvec(1),'r'),shg
figure
hist(s,100)
clf
hist(s',100)
s = [.1 ones(1,98)/2 1];

hist(s',100)
A = randspd(s);
[x,flag,relres,iter,resvec] = pcg(A,b,1e-8,n);
figure(1)
semilogy(resvec / resvec(1),'g'),shg
x = [.1  .5+.01*randn(1,98) 1];
s = x;
size(s)
ans =
     1   100
figure(2)
hist(s',100)
A = randspd(s);
[x,flag,relres,iter,resvec] = pcg(A,b,1e-8,n);
figure(1)
semilogy(resvec / resvec(1),'k'),shg
x = [.1  .5+.2*randn(1,98) 1];
s=x;
figure(2)
hist(s',100)
A = randspd(s);
[x,flag,relres,iter,resvec] = pcg(A,b,1e-8,n);
figure(1)
semilogy(resvec / resvec(1),'k'),shg