ich habe ein Problem mit pyamg.
Ich würde das gerne auf ein lineares Gleichungssystem mit sehr schlecht konditionierter Matrix anwenden. Ich habe das mal mit Jacobi und Gauss-Seidel ausprobiert und meine neue Matrix ist auch besser konditioniert, aber die Lösung des Verfahrens ist trotzdem noch sehr fehlerbehaftet.
Leider habe ich aber das Programm nicht richtig durchblickt und weiß daher nicht ob ich den Fehler vielleicht falsch berechnet habe.
Hier ist erstmal der Link zur Dokumentation:
http://pyamg.googlecode.com/svn/branche ... ation.html
Da ist als Beispiel folgender Code:
Code: Alles auswählen
>>> ## Use Jacobi as a Stand-Alone Solver
>>> from pyamg.relaxation import *
>>> from pyamg.gallery import poisson
>>> from pyamg.util.linalg import norm
>>> import numpy
>>> A = poisson((10,10), format='csr')
>>> x0 = numpy.zeros((A.shape[0],1))
>>> b = numpy.ones((A.shape[0],1))
>>> jacobi(A, x0, b, iterations=10, omega=1.0)
>>> print norm(b-A*x0)
5.83475132751
>>> #
>>> ## Use Jacobi as the Multigrid Smoother
>>> from pyamg import smoothed_aggregation_solver
>>> sa = smoothed_aggregation_solver(A, B=numpy.ones((A.shape[0],1)),
... coarse_solver='pinv2', max_coarse=50,
... presmoother=('jacobi', {'omega': 4.0/3.0, 'iterations' : 2}),
... postsmoother=('jacobi', {'omega': 4.0/3.0, 'iterations' : 2}))
>>> x0=numpy.zeros((A.shape[0],1))
>>> residuals=[]
>>> x = sa.solve(b, x0=x0, tol=1e-8, residuals=residuals)
Wenn ich das ganze jetzt als Vorkonditionierer mache, wie ich es hier versucht habe:
Code: Alles auswählen
from pyamg.relaxation import *
from pyamg.gallery import poisson
from pyamg.util.linalg import norm
from scipy.sparse.linalg import cg
from scipy.linalg import lu
import numpy
from scipy import *
from scipy.sparse import csr_matrix
A = numpy.matrix([[1,2,3],[2,3,1],[1,2,4]])
b = numpy.ones((A.shape[0],1))
from pyamg import smoothed_aggregation_solver
sa = smoothed_aggregation_solver(A, B=numpy.ones((A.shape[0],1)),
coarse_solver='pinv2', max_coarse=50,
presmoother=('jacobi', {'omega': 4.0/3.0, 'iterations' : 2}),
postsmoother=('jacobi', {'omega': 4.0/3.0, 'iterations' : 2}))
M = sa.aspreconditioner(cycle='V')
B=M.matmat(b)
a=M.matmat(A)
x=spsolve(a,B)
print "Fehler:", np.dot(a,x)-B