Initiation au Calcul Scientifique Intensif¶
Rappels sur les matrices
Structure des matrices
Illustration du remplissage par la multiplication matricielle
import matplotlib.pylab as plt
import scipy.sparse as sps
import scipy.linalg as spl
import numpy as np
n=100
A = sps.rand(n,n, density=0.05)
show('creux : ', (len(sps.find(A)[0])/n^2*100).n(),'%')
A=A.todense()
plt.spy(A,precision=0.01, markersize=1)
plt.show();
B=np.dot(A,A)
show('creux : ', (len(sps.find(B)[0])/n^2*100).n(),'%')
plt.spy(B,precision=0.01, markersize=1)
plt.show();
C=np.dot(B,A)
show('creux : ', (len(sps.find(C)[0])/n^2*100).n(),'%')
plt.spy(C,precision=0.01, markersize=1)
plt.show();
D=np.dot(C,A)
show('creux : ', (len(sps.find(D)[0])/n^2*100).n(),'%')
plt.spy(D,precision=0.01, markersize=1)
plt.show();
Illustration de l'effet de la renumérotation
# Matrice aléatoire au format CSR
ran=sps.random(100,100,.05, format='csr')
A = sps.csr_matrix((np.ones(ran.data.shape).astype(np.uint8),ran.indices, ran.indptr))
plt.spy(A,precision=0.01, markersize=1)
plt.show();
# Largeur de bande
i,j=A.nonzero()
show(np.max(i-j))
# Tableau de renumérotation après algorithme Reverse Cuthill MacKee
perm = sps.csgraph.reverse_cuthill_mckee(A, False)
# Permutation
B=A[np.ix_(perm,perm)].A
plt.spy(B,precision=0.01, markersize=1)
plt.show();
# Largeur de bande
i,j=B.nonzero()
show(np.max(i-j))
Illustration de l'effet du conditionnement
A= matrix(4,4,[10,7,8,7,7,5,6,5,8,6,10,9,7,5,9,10])
b1=vector([32,23,33,31])
show('A=',A,'b=',b1)
sol1=A\b1
show(sol1)
b2=vector([32.1,22.9,33.1,30.9])
show('b=',b2)
sol2=A\b2
show(sol2.n())
err1=(b1-b2).norm()/b1.norm()
show(err1.n())
err2=(sol1-sol2).norm()/sol1.norm()
show(err2.n())
show(err2.n()/err1.n())
L'erreur commise sur la solution du système perturbée est 2460 fois plus importante que l'erreur commise sur le second membre !
Calculons le conditionnement de la matrice : A est symétrique et ses valeurs propres sont positives, son conditionnement en norme 2 est le rapport des la valeur propre maximale sur la valeur propre minimale.
show(A.eigenvalues())
condA=max(A.eigenvalues())/min(A.eigenvalues())
show(condA)
import numpy as np
A=np.eye(n,n);
for i in range(n):
A[i,n-1]=1
A[n-1,i]=1
plt.spy(A)
plt.show()
P, L, U = spl.lu(A)
plt.spy(L)
plt.show()
plt.spy(U)
plt.show()
import numpy as np
A=np.eye(n,n);
for i in range(n):
A[i,0]=1
A[0,i]=1
plt.spy(A)
plt.show()
P, L, U = spl.lu(A)
plt.spy(L)
plt.show()
plt.spy(U)
plt.show()
Illustration de l'importance du pivotage
M=matrix([[0,1],[1,1]])
y=matrix([[1],[2]])
show(M,'x=',y)
x=M\y
show('x=',x)
var('eps')
A=matrix([[eps,1],[1,1]])
a=matrix([[1],[2]])
show(A,'x1=',a)
B=matrix([[1,1],[eps,1]])
b=matrix([[2],[1]])
show(B,'x2=',b)
x1=(A\a); show('x1=',x1)
x2=(B\b); show('x2=',x2)
En terme de représentation machine, le premier résultat est (0,1), tandis que le deuxième est bien (1,1) !
1-1-3- Décomposition QR¶
Illustration du procédé de Gram-Schmidt modifié
Procédure de Gram-Schmidt pour une matrice 3x3.
Voir exercice 2 pour les algorithmes.
def GramSchmidt(A):
v1=vector(A.transpose()[0])
v2=vector(A.transpose()[1])
v3=vector(A.transpose()[2])
u1=v1;
n1=u1.norm();
q1=(u1/n1)
u2=v2 - u1.dot_product(v2)/n1^2*u1;
n2=u2.norm()
q2=(u2/n2)
u3=v3 - u1.dot_product(v3)/n1^2*u1 - u2.dot_product(v3)/n2^2*u2;
n3=u3.norm()
q3=(u3/n3)
r1=u1.dot_product(v3)/n1^2
r2=u2.dot_product(v3)/n2^2
# Base orthogonale obtenue : (q1,q2,q3)
return q1, q2, q3, r1, r2
var('epsilon')
assume(epsilon>0)
A=matrix([[1,1,1],[epsilon,0,0],[0,epsilon,0]]);
show('A=',A)
[q1,q2,q3,r1,r2]=GramSchmidt(A)
show(q1.simplify_full(),q2.simplify_full(),q3.simplify_full())
show('Verification de l orthogonalite')
show('(q1,q2)=',(q1.dot_product(q2)).simplify_full())
show('(q1,q3)=',(q1.dot_product(q3)).simplify_full())
show('(q2,q3)=',(q2.dot_product(q3)).simplify_full())
Regardons de plus près les valeurs suivantes
r1 = proj_{u1}(v3) = < u1,v3 >/< u1,u1 >
r2 = proj_{u2}(v3) = < u2,v3 >/< u2,u2 >
show('proj_{u1}(v3)')
show(r1)
show(limit(r1,epsilon=0))
show(r1(epsilon=1e-10))
show('proj_{u2}(v3)=')
show(r2)
show(r2.simplify_full())
show(limit(r2,epsilon=0))
show(r2(epsilon=1e-10))
Le terme de projection de v3 sur u2 est censé valoir 1/2 mais est numériquement nul !
Le vecteur proj_{u2}(v3) ne contribue donc pas au calcul du vecteur de base comme il le devrait.
epsilon=1e-10
A=matrix([[1,1,1],[epsilon,0,0],[0,epsilon,0]]);
[q1,q2,q3,R1,R2]=GramSchmidt(A)
show(q1,q2,q3)
show('Verification de l orthogonalite')
show('(q1,q2)=',q1.dot_product(q2))
show('(q1,q3)=',q1.dot_product(q3))
show('(q2,q3)=',q2.dot_product(q3))
La base otenue n'est pas orthogonale !
def GramSchmidt_Modif(A):
v1=vector(A.transpose()[0])
v2=vector(A.transpose()[1])
v3=vector(A.transpose()[2])
u1=v1;
n1=u1.norm();
q1=(u1/n1)
u2=v2 - u1.dot_product(v2)/n1^2*u1;
n2=u2.norm()
q2=(u2/n2)
u31=v3 - u1.dot_product(v3)/n1^2*u1;
u3= u31 - u2.dot_product(u31)/n2^2*u2;
n3=u3.norm()
q3=(u3/n3)
R1=u1.dot_product(v3)/n1^2
R2=u2.dot_product(u31)/n2^2
return q1, q2, q3, R1, R2
[q1,q2,q3,R1,R2]=GramSchmidt_Modif(A)
show(q1,q2,q3)
show('Verification de l orthogonalite')
show('(q1,q2)=',q1.dot_product(q2))
show('(q1,q3)=',q1.dot_product(q3))
show('(q2,q3)=',q2.dot_product(q3))
Avec l'algorithme modifié, la base obtenue est bien orthogonale.
show(r2.simplify_full())
show(R2)
En effet, dans ce cas, le terme de projection proj_{u2}(v3) est correctement approché.
1-2- Méthodes itératives classiques¶
Illustration de l'influence de omega sur la convergence de l'algorithme de la relaxation
A=matrix(5,[1., 0.1, 0.01, 0.001, 0.0001, 1., 1., 1., 1., 1., 1., 1.5, 2.25, 3.375, 5.0625, 1., 2., 4., 8., 16., 1., 3., 9., 27., 81.])
b=vector([1., 1.5, 2.25, 3.375, 5.0625])
show('A=',A)
show('b=',b)
A\b
Algorithme de la relaxation
def relax(A,b,x0,niter,omega):
n=(A.dimensions())[0];
D=matrix(RR,n)
E=matrix(RR,n)
F=matrix(RR,n)
for i in range(n):
D[i,i]=A[i,i]
for j in range(i):
E[i,j]=-A[i,j];
F[j,i]=-A[j,i];
sol=A\b;
u=x0;
res=[x0];
resn=[norm(u-sol)];
for k in range(niter):
u=(1/omega*D-E)\(((1-omega)/omega*D+F)*u+b);
res=res+[u];
resn=resn+[norm(u-sol)];
return resn;
Convergence logarithmique de la méthode de la relaxation pour différentes valeurs de omega
x0=vector(RR,[1,5,1,5,1])
n=50
resn05=relax(A,b,x0,n,0.5);
resn08=relax(A,b,x0,n,0.8);
resn12=relax(A,b,x0,n,1.2);
resn14=relax(A,b,x0,n,1.4);
resn16=relax(A,b,x0,n,1.6);
resn19=relax(A,b,x0,n,1.9);
p=point((1,log(resn05[0])),color='green')
p=point((1,log(resn08[0])),color='red')
p=point((1,log(resn12[0])),color='blue')
p=point((1,log(resn14[0])),color='cyan')
p=point((1,log(resn16[0])),color='magenta')
p=point((1,log(resn19[0])),color='black')
for i in range(n):
p=p+point((i,log(resn05[i])),color='green')
p=p+point((i,log(resn08[i])),color='red')
p=p+point((i,log(resn12[i])),color='blue')
p=p+point((i,log(resn14[i])),color='cyan')
p=p+point((i,log(resn16[i])),color='magenta')
p=p+point((i,log(resn19[i])),color='black')
p.show()
La meilleur valeur dans ce cas est donc omega=1.6
1-3- Méthodes du gradient¶
Illustration du gradient à pas constant
def f(x,y):
return x^2+2*y^2+x*y+x-y
def df(x,y):
return [2*x+y+1, x+4*y-1]
Le minimium est (-5/7,3/7), où le gradient est nul
import numpy as np
import matplotlib.pylab as plt
X=np.arange(-3,2,0.01)
Y=np.arange(-2,2.5,0.01)
[X,Y]=np.meshgrid(X,Y)
Z=f(X,Y)
plt.axis('equal')
plt.contour(X,Y,Z,[i*0.1 for i in srange(-100,40,5)])
plt.scatter(-5/7,3/7)
plt.show()
def grad(f,df,rho,X0,Y0,tol,N):
err=1; it=1;
X=X0; Y=Y0; Xs=[X]; Ys=[Y];
while ((err>tol) and (it<N)):
Xt = X - rho*df(X,Y)[0]
Yt = Y - rho*df(X,Y)[1]
err=abs(f(Xt,Yt)-f(X,Y))
X=Xt; Y=Yt;
Xs=Xs+[X]; Ys=Ys+[Y];
it=it+1;
if it==N:
show('Algorithm stopped before reaching convergence')
return Xs,Ys
[Xs,Ys]=grad(f,df,0.1,0,0,1E-10,100)
show([Xs[-1]+5/7,Ys[-1]-3/7])
Le résultat numérique est proche de la solution exacte
plt.axis('equal')
plt.contour(X,Y,Z,[i*0.1 for i in srange(-100,40,5)])
plt.scatter(Xs,Ys)
plt.show()
1-4- Méthodes de projection sur des espaces de Krylov¶
Illustration de la dégénérescence numérique de la base naturelle des espaces de Krylov
D=2*matrix([[1,0,0,0,0],[0,10,0,0,0],[0,0,100,0,0],[0,0,0,1000,0],[0,0,0,0,10000]])
P=matrix([[1,2,0,0,1],[0,1,0,-1,1],[2,0,1,0,1],[0,2,0,1,1],[1,-2,1,1,-1]])
A=P*D*P^(-1)
show('A=',A)
r0=vector([1,1,1,1,1])
show(r0)
show(A*r0)
show(A^2*r0)
show(A^3*r0)
show(A^4*r0)
def MethodeQR(A,n):
for i in range(n):
[Q,R]=A.QR();
A=R*Q;
return A;
A=matrix(QQbar,[[-17,36,27],[-33,34,33],[-84,102,94]])
show('A=',A)
show('Sp(A)=',A.eigenvalues())
An=MethodeQR(A,5);
show('An=',An)
3- Résolution de systèmes non linéaires¶
Illustration du bassin de convergence de l'algorithme de Newton
On s’intéresse à la recherche des solutions complexes de l’équation z^3 = 1. Pour cela nous allons utiliser la méthode de Newton appliquée à la fonction
F : (x,y)->(Re((x+i y)^3-1),Im((x+i y)^3-1)) qui s'annule exactement lorsque (x+i y)^3=1.
Calculons explicitement la fonction F et sa différentielle DF
var('x y')
F=vector([-1+x^3-3*x*y^2,3*x^2*y-y^3])
DF=matrix([[3*x^2-3*y^2,-6*x*y],[6*x*y,3*x^2-3*y^2]])
show('F=',F,', DF=',DF)
L'algorithme de Newton en 2D s'écrit de la manière suivante
def Newton(F,DF,x0,y0,eps):
z=vector([x0,y0]);
err=[norm(DF(x=z[0],y=z[1])\F(x=z[0],y=z[1]))];
n=0;
while ((norm(F(x=z[0],y=z[1]))>eps) and (n<100)):
z = z - (DF(x=z[0],y=z[1])\F(x=z[0],y=z[1]));
err=err+[norm(DF(x=z[0],y=z[1])\F(x=z[0],y=z[1]))];
n=n+1;
return n,z;
On travaille sur une grille couvrant le carré [−1.5, 1.5] × [−1.5, 1.5] de pas h = 3/n .
On va résoudre F(x, y) = 0 par la méthode de Newton initialisée successivement en chacun des points (x_i , y_j ) de cette grille.
Dans un tableau, on mémorisera le numéro k entre 1 et 3 de la racine exp(2ikπ/3) vers lequel l’algorithme aura convergé.
n=30;
T=matrix(n,n);
l=3;
h=l/n;
racines=[(cos(0),sin(0)),(cos(2*pi/3).n(),(sin(2*pi/3)).n()),(cos(4*pi/3).n(),(sin(4*pi/3)).n())];
tol=0.001;
for i in range(n):
xg=-l/2+h*(i-0.5);
for j in range(n):
yg=-l/2+h*(j-0.5);
[iter,z]=Newton(F,DF,xg,yg,tol);
for k in range(3):
if ((abs(racines[k][0]-z[0])<tol) and (abs(racines[k][1]-z[1])<tol)):
T[i,j]=k;
p=plot((0,0),color='red')
for i in range(n):
xg=-l/2+h*(i-0.5);
for j in range(n):
yg=-l/2+h*(j-0.5);
if (T[i,j]==0):
p=p+point((xg,yg),color='red')
if (T[i,j]==1):
p=p+point((xg,yg),color='green')
if (T[i,j]==2):
p=p+point((xg,yg),color='yellow')
show(p)
En raffinant, on obtiendrait une figure fractale comme celle-ci !
from IPython.display import Image
Image("Newton3.png")
Avec l'équation z^5-1=0, on obtiendrait la figure suivante
from IPython.display import Image
Image("Newton5.png")
Illustration de l'ordre de convergence des méthodes de résolution de systèmes non linéaires en 1D
def dichotomie(f,a,b,eps):
c=(a+b)/2.;
it=[c];
n=0;
while ((abs(f(c))>eps) and (n<100)):
n=n+1;
if f(a)*f(c)<0:
b=c;
else:
a=c;
c=(a+b)/2.
it=it+[c];
return n, c, it;
def secante(f,x0,x1,eps):
xt=x0;x=x1;
err=[(x-xt)/(f(x)-f(xt))];
n=0;
while ((abs(f(x))>eps) and (n<100)):
c = (x-xt)/(f(x)-f(xt));
xt = x;
x = x - c*f(x);
err=err+[abs(f(xt)*c)];
n=n+1;
return n,x,err;
def newton(f,df,x0,eps):
x=x0; err=[abs(f(x)/df(x))];
n=0;
while ((abs(f(x))>eps) and (n<100)):
x = x - f(x)/df(x);
err=err+[abs(f(x)/df(x))];
n=n+1;
return n,x,err;
On cherche la valeur approchée de sqrt(2) à 1E-15 près en utilisant la méthode de la dichotomie, de Newton et de la sécante.
var('x')
function('f')(x)
f(x)=x^2-2;
function('df')(x)
df(x)=2*x;
a=0; b=2;
x0=1.; x1=1.1;
eps=1E-15;
[n1,c1,err1]=dichotomie(f,a,b,eps)
show('Methode de la dichotomie')
show('nombre d iterations : ',n1)
show('Valeur obtenue : ',c1)
show('Erreur : ',abs(c1-sqrt(2.)))
[n2,c2,err2]=newton(f,df,x0,eps)
show('Methode de Newton')
show('nombre d iterations : ',n2)
show('Valeur obtenue : ',c2)
show('Erreur : ',abs(c2-sqrt(2.)))
[n3,c3,err3]=secante(f,x0,x1,eps)
show('Methode de la secante')
show('nombre d iterations : ',n3)
show('Valeur obtenue : ',c3)
show('Erreur : ',abs(c3-sqrt(2.)))
Dichotomie
On trace la courbe de convergence de l'algorithme : on sait que la suite d'approximation (x(n+1)) approchant la valeur alpha vérifie || x(n+1) - alpha || < C || x_n - alpha ||, c'est à dire que la méthode est d'ordre 1.
En traçant log(|| x_(n+1) - alpha ||) en fonction de log(|| x_n - alpha ||), on devrait obtenir une droite de pente 1
p=point((log(abs(err1[0]-sqrt(2))),log(abs(err1[1]-sqrt(2)))))
for i in range(n1):
p=p+point((log(abs(err1[i-1]-sqrt(2))),log(abs(err1[i]-sqrt(2)))))
p=p+plot(x,[-40,0],color="green")
show(p)
Newton
On trace la courbe de convergence de l'algorithme : on sait que la suite d'approximation (x(n+1)) approchant la valeur alpha vérifie || x(n+1) - x_n || < C || xn - x(n-1) ||^2, c'est à dire que la méthode est d'ordre 2.
En traçant log(|| x_(n+1) - x_n ||) en fonction de log(|| xn - x(n-1) ||), on devrait obtenir une droite de pente 2.
p=point((log(err2[0]),log(err2[1])))
for i in range(n2):
p=p+point((log(err2[i-1]),log(err2[i])))
p=p+plot(2*x,[-20,0],color="green")
show(p)
Sécante
On trace la courbe de convergence de l'algorithme : on sait que la suite d'approximation (x(n+1)) approchant la valeur alpha vérifie || x(n+1) - x_n || < C || xn - x(n-1) ||^alpha, où alpha=(1+sqrt(5))/2 est le nombre d'or, c'est à dire que la méthode est d'ordre alpha.
En traçant log(|| x_(n+1) - x_n ||) en fonction de log(|| xn - x(n-1) ||), on devrait obtenir une droite de pente alpha.
p=point((log(err3[0]),log(err3[1])))
for i in range(n3):
p=p+point((log(err3[i-1]),log(err3[i])))
p=p+plot((1+sqrt(5.))/2*x,[-10,0],color="green")
show(p)
Estimation du conditionnement¶
Illustration de certains estimateurs pour le calcul du conditionnement
import numpy as np
import scipy.linalg as spl
def cond_app(A):
P, L, U = spl.lu(A)
n=np.size(A,1);
p=np.zeros(n);
z=np.zeros(n);
for k in range(n):
zplus =(1-p[k])/U[k,k];
splus=abs(1-p[k]);
zminus=(-1-p[k])/U[k,k];
sminus=abs(-1-p[k]);
for i in srange(k+1,n):
splus =splus+abs(p[i]+U[k,i]*zplus);
sminus=sminus+abs(p[i]+U[k,i]*zminus);
if splus > sminus:
z[k]=zplus;
else:
z[k]=zminus;
for i in srange(k+1,n):
p[i]=p[i]+U[k,i]*z[k];
z=z.transpose();
x=np.linalg.solve(L.transpose(), z);
w=np.linalg.solve(L, x);
y=np.linalg.solve(U, w);
return(np.linalg.norm(A,1)*np.linalg.norm(y,1)/np.linalg.norm(x,1));
A=matrix(4,4,[10,7,8,7,7,5,6,5,8,6,10,9,7,5,9,10]);
show(A)
show(np.linalg.cond(A,1))
show(cond_app(A))
# Matrice Hilbert
n=10;
H=matrix(QQ,n,n);
for i in range(n):
for j in range(n):
H[i,j]=1/(i+j+1);
show(H)
show(np.linalg.cond(H,1))
show(cond_app(H))
A=matrix([[10,7,8,7],[7,5,6,5],[8,6,10,9],[7,5,9,10]])
D=matrix([[10,0,0,0],[0,5,0,0],[0,0,10,0],[0,0,0,10]])
E=-matrix([[0,0,0,0],[7,0,0,0],[8,6,0,0],[7,5,9,0]])
F=E.transpose()
show(A)
condA=max(A.eigenvalues())/min(A.eigenvalues())
show(condA)
M=D;
A1=M^(-1)*A;
condA1=max(A1.eigenvalues())/min(A1.eigenvalues())
show(condA1)
Comparons la valeur du conditionnement avec différents préconditionneurs
p=plot(condA,[0,2],color='green')
p=p+plot(condA1,[0,2],color='red')
Lomega=[0.2,0.5,0.8,1,1.2,1.5,1.8]
for i in range(7):
omega=Lomega[i]
M=omega/(2-omega)*(1/omega*D-E)*D^(-1)*(1/omega*D-F);
#show(max(M.eigenvalues())/min(M.eigenvalues()))
A2=M^(-1)*A;
condA2=max(A2.eigenvalues())/min(A2.eigenvalues())
show(condA2)
p=p+point((omega,condA2))
show(p)
Sur la figure, on a tracé
- en vert la valeur du conditionnement initial
- en rouge la valeur du conditionnement avec préconditionneur de Jacobi
- en bleu les valeurs du conditionnement avec SSOR pour plusieurs valeurs de omega
Illustration des effets des préconditionneurs classiques sur la matrice du Laplacien
n=20;
dx=1/(n-1);
A=2*np.diag(np.ones(n))-np.diag(np.ones(n-1),k=1)-np.diag(np.ones(n-1),k=-1)
plt.spy(A,precision=0.01, markersize=1)
plt.show()
D=np.diag(np.diag(A),0);
E=-np.tril(A,-1);
F=-np.triu(A,1);
# Conditionnement
condA=np.linalg.cond(A,1)
# Jacobi
M=D;
J=np.dot(np.linalg.inv(M),A);
condJ=np.linalg.cond(J)
p=plot(condJ/condA,[0,2],color='green')
# SSOR
for omega in [0.2,0.4,0.6,0.8,1,1.2,1.4,1.6,1.8]:
M=np.dot(omega/(2-omega)*(1/omega*D-E),np.linalg.inv(D))
M=np.dot(M,(1/omega*D-F));
A2=np.dot(np.linalg.inv(M),A);
cond2=np.linalg.cond(A2)
p=p+point((omega,cond2/condA))
# ILU
P, L, U = spl.lu(A)
for i in range(n):
for j in range(n):
if A[i,j]==0:
L[i,j]==0;
U[i,j]==0
LU=np.dot(np.linalg.inv(L),np.linalg.inv(U))
ILU=np.dot(np.linalg.inv(LU),np.linalg.inv(A))
condILU=np.linalg.cond(ILU)
p+=plot(condILU/condA,[0,2],color='red')
# MILU
for i in range(n):
s=sum(A[i,j] for j in range(n))
U[i,i]+=s
LU=np.dot(np.linalg.inv(L),np.linalg.inv(U))
MILU=np.dot(np.linalg.inv(LU),np.linalg.inv(A))
condMILU=np.linalg.cond(MILU)
p+=plot(condMILU/condA,[0,2],color='black')
show(p)
np.linalg.cond(A)
Sur la figure, on a tracé
- en vert le gain du conditionnement avec préconditionneur de Jacobi
- en bleu les gain du conditionnement avec SSOR pour plusieurs valeurs de omega
- en rouge le gain du conditionnement avec préconditionneur ILU
- en noir le gain du conditionnement avec préconditionneur MILU
5- Stockage creux¶
Mise en pratique
A=matrix([[7,0,0,1,3],[6,1,2,0,0],[0,2,4,0,1],[5,0,0,0,2],[0,0,8,0,0]])
x=matrix([1,2,3,4,5]); x=x.transpose()
show('A=',A,' x=', x)
valA=matrix([7,1,3,6,1,2,2,4,1,5,2,8]); valA=valA.transpose();
col_indA=matrix([0,1,4,5,1,2,3,2,3,5,1,5,3]); col_indA=col_indA.transpose();
row_ptrA=matrix([0,3,6,9,11,12]); row_ptrA=row_ptrA.transpose()
def pdmatvect(A,x):
n=(A.dimensions())[0];
y=matrix(n,1);
for i in range(n):
for j in range(n):
y[i]=y[i]+A[i,j]*x[j]
return y;
y=pdmatvect(A,x)
show('y_th=',A*x,'y=',y)
def pdmatvect_CSR(valA,col_indA,row_ptrA,x):
n=x.dimensions()[0];
y=matrix(n,1);
for i in range(n):
for j in range(row_ptrA[i][0],row_ptrA[i+1][0]):
show(i,j)
y[i]=y[i]+valA[j][0]*x[col_indA[j+1][0]-1]
return y;
y=pdmatvect_CSR(valA,col_indA,row_ptrA,x)
show('y_th=',A*x,'y=',y)
def pdmatvect(A,x):
n=(A.dimensions())[0];
y=matrix(n,1);
for i in range(n):
for j in range(n):
y[i]=y[i]+A[i,j]*x[j]
return y;
def pdmatvect_CSR(valA,col_indA,row_ptrA,x):
n=x.dimensions()[0];
y=matrix(n,1);
for i in range(n):
for j in range(row_ptrA[i][0],row_ptrA[i+1][0]):
y[i]=y[i]+valA[j][0]*x[col_indA[j+1][0]-1]
return y;