Residencia mínima general (GMRES) con precondicionador ILU

Estoy tratando de implementar el preacondicionador ILU en este código GMRES que escribí (para resolver el sistema lineal Ax = b. Estoy tratando con una matriz tridigonal SPD de dimensión 25×25. Como puede ver, estoy calculando el precondicionador con método spilu. El código se ejecuta sin errores, pero la solución es claramente errónea ya que, al final del código, estoy imprimiendo la norma de b y la norma del producto A * x. No son casi las igual .. El código Ejecutar bien sin acondicionador previo y converger con 13 iteraciones para la misma matriz.

Este es el código que he seguido.

import numpy as np import scipy as sp import matplotlib.pyplot as plt 'Size controller' matrixSize =25 'Building a tri-diagonal matrix' def Atridiag(val_0, val_sup, val_inf, mSize): cen = np.ones((1, mSize))*val_0 sup = np.ones((1, mSize-1))*val_sup inf = np.ones((1, mSize-1))*val_inf diag_cen = np.diagflat(cen, 0) diag_sup = np.diagflat(sup, 1) diag_inf = np.diagflat(inf, -1) return diag_cen + diag_sup + diag_inf A = Atridiag(2, -1, -1, matrixSize) A = sp.sparse.csc_matrix (A) 'Plot matrix sparsity' plt.clf() plt.spy(A, marker ='.', markersize=2) plt.show() 'random b and x0 vectors' b = np.matrix(np.ones((matrixSize, 1))) x = np.matrix(np.ones((matrixSize, 1))) 'Incomplete LU' M = sp.sparse.linalg.dsolve.spilu(A) M1 = lambda x: M.solve(x) M2=sp.sparse.linalg.LinearOperator((matrixSize,matrixSize),M1) 'Initial Data' nmax_iter = 30 rstart = 2 tol = 1e-7 e = np.zeros((nmax_iter + 1, 1)) rr = 1 'Starting GMRES' for rs in range (0, rstart+1): 'first check on residual' if rr < tol : break else: r0 = (b - A.dot(x)) betha = np.linalg.norm(r0) e[0] = betha H = np.zeros((nmax_iter + 1, nmax_iter)) V = np.zeros((matrixSize, nmax_iter+1)) V[:, 0:1] = r0/betha for k in range (1, nmax_iter+1): 'Appling the Preconditioner' t = A.dot(V[:, k-1]) V[:, k] = M2.matvec(t) 'Ortogonalizzazione GS' for j in range (k): H[j, k-1] = np.dot(V[:, k].T, V[:, j]) V[:, k] = V[:, k] - (np.dot(H[j, k-1], V[:, j])) H[k, k-1] = np.linalg.norm(V[:, k]) V[:, k] = V[:, k] / H[k, k-1] 'QR Decomposition' n=k Q = np.zeros((n+1, n)) R = np.zeros((n, n)) R[0, 0] = np.linalg.norm(H[0:n+2, 0]) Q[:, 0] = H[0:n+1, 0] / R[0,0] for j in range (0, n+1): t = H[0:n+1, j-1] for i in range (0, j-1): R[i, j-1] = np.dot(Q[:, i], t) t = t - np.dot(R[i, j-1], Q[:, i]) R[j-1, j-1] = np.linalg.norm(t) Q[:, j-1] = t / R[j-1, j-1] g = np.dot(QT, e[0:k+1]) Z = np.dot(np.linalg.inv(R), g) Res = e[0:n] - np.dot(H[0:n, 0:n], Z[0:n]) rr = np.linalg.norm(Res) 'second check on residual' if rr < tol: break 'Updating the solution' x = x + np.dot(V[:, 0:k], Z) print(sp.linalg.norm(b)) print(sp.linalg.norm(np.dot(A.todense(),x))) 

Realmente espero que alguien pueda resolverlo !!

Tal vez sea demasiado tarde, pero para futuras referencias:

Olvidaste multiplicar por el acondicionador al actualizar x:

 x = x + M2.dot(np.dot(V[:, 0:k], Z) # M2.matvec() works the same 

Ver aqui

Con esa corrección, el algoritmo converge en 1 iteración.


Otros comentarios:

  • Puede hacerlo directamente: M2 = sp.sparse.linalg.LinearOperator((matrixSize,matrixSize),M.solve)
  • Al final, para comparar Ax y b , es mejor imprimir la diferencia (residual) porque obtendrá un resultado mucho más preciso: print(sp.linalg.norm(b - np.dot(A.todense(),x)))