forked from skambuilds/IMAT
-
Notifications
You must be signed in to change notification settings - Fork 0
/
SelfPreConiugGradient.m
45 lines (44 loc) · 1.29 KB
/
SelfPreConiugGradient.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
% Metodo del gradiente coniugato precondizionato
% x - soluzione calcolata
% k - numero di iterazioni eseguite
% resvec - vettore che contiene il residuo ad ogni iterazione
function [x,k,resvec] = SelfPreConiugGradient(A,b,tau,maxn,Rt,R,x)
% Determino la dimensione della matrice
n=size(A,1);
% Inizializzo il vettore x0
x0=100*ones(n,1);
% Inizializzazione residuo
r = b - A*x;
% Risoluzine di un sistema lineare
y=Rt\r;
z=R\y;
% Inizializzazione direzione iniziale
p = z;
% Inizializzazione contatore iterazioni
k = 0;
% Avvio ciclo dell'algoritmo con condizione di controllo, posso scegliere
% tra:
% - controllo del residuo: norm(r)>tau*norm(b)
% - condizione di Cauchy: norm(x-x0)>tau*norm(x)
% Imposto comunque un numero di iterazioni massimo
% Preallocazione risorse per il vettore residuo
resvec=zeros(maxn,1);
while(norm(x-x0)>tau*norm(x)) && (k<maxn)
x0=x;
k=k+1;
% Memorizzo la norma del residuo nel vettore
resvec(k)=norm(r);
% Ottimizzo calcolando solo una volta il prodotto matrice vettore:
s = A*p;
delta = p'*s;
% Calcolo del passo:
alpha = (p'*r)/delta;
% Calcolo nuovo soluzione
x=x0+alpha*p;
% Aggiornamento residuo
r=r-alpha*s;
y=R'\r;
z=R\y;
beta = s'*z/delta;
p = z-beta*p;
end