1

I am trying to implement the algorithm of GMRES with right-preconditioner P for solving the linear system Ax = b enter image description here. The code is running without error; however, it pops into unprecise result for me because the error I have is very large. For the GMRES method (without preconditioning matrix - remove P in the algorithm), the error I get is around 1e^{-12} and it converges with the same matrix.

import numpy as np
from scipy import sparse
import matplotlib.pyplot as plt
from scipy.linalg import norm as norm
import scipy.sparse as sp
from scipy.sparse import diags

"""The program is to split the matrix into D-diagonal; L: strictly lower matrix; U strictly upper matrix
    satisfying: A = D - L - U  """
def splitMat(A):
    n,m = A.shape
    if (n == m):
        diagval = np.diag(A)
        D = diags(diagval,0).toarray()
        L = (-1)*np.tril(A,-1)
        U = (-1)*np.triu(A,1)
    else: 
        print("A needs to be a square matrix")
    return (L,D,U)

"""Preconditioned Matrix for symmetric successive over-relaxation (SSOR): """
def P_SSOR(A,w): 
    ## Split up matrix A: 
    L,D,U = splitMat(A)
    Comp1 = (D - w*U)
    Comp2 = (D - w*L)
    Comp1inv = np.linalg.inv(Comp1)
    Comp2inv = np.linalg.inv(Comp2)
    P = w*(2-w)*np.matmul(Comp1inv, np.matmul(D,Comp2inv))
    return  P

"""GMRES_SSOR using right preconditioning P: 
A - matrix of linear system Ax = b
x0 - initial guess
tol - tolerance
maxit - maximum iteration """

def myGMRES_SSOR(A,x0, b, tol, maxit):
    matrixSize = A.shape[0] 
    e = np.zeros((maxit+1,1))
    rr = 1
    rstart = 2
    X = x0
    w = 1.9 ## in ssor
    P = P_SSOR(A,w) ### preconditioned matrix 
    ### Starting the GMRES #### 
    for rs in range(0,rstart+1):
        ### first check the residual: 
        if rr<tol: 
            break 
        else: 
            r0 = (b-A.dot(x0))
            rho = norm(r0)
            e[0] = rho 
            H = np.zeros((maxit+1,maxit))
            Qcol = np.zeros((matrixSize, maxit+1))
            Qcol[:,0:1] = r0/rho 
        for k in range(1, maxit+1):
            ### Arnodi procedure ##
            Qcol[:,k] =np.matmul(np.matmul(A,P), Qcol[:,k-1])  ### This step applies P here: 
            for j in range(0,k):
                H[j,k-1] = np.dot(np.transpose(Qcol[:,k]),Qcol[:,j])
                Qcol[:,k] = Qcol[:,k] - (np.dot(H[j,k-1], Qcol[:,j]))
            
            H[k,k-1] =norm(Qcol[:,k])
            Qcol[:,k] = Qcol[:,k]/H[k,k-1]

            ###  QR decomposition step ### 
            n = k 
            Q = np.zeros((n+1, n))
            R = np.zeros((n, n))
            R[0, 0] = 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] = norm(t)
                Q[:, j-1] = t / R[j-1, j-1]

            g = np.dot(np.transpose(Q), e[0:k+1])
            Y = np.dot(np.linalg.inv(R), g) 

            Res= e[0:n] - np.dot(H[0:n, 0:n], Y[0:n])
            rr = norm(Res) 

            #### second check on the residual ###
            if rr < tol: 
                break 

        #### Updating the solution with the preconditioned matrix ####
        X = X  + np.matmul(np.matmul(P,Qcol[:, 0:k]), Y) ### This steps applies P here: 
    return X 

######
A = np.random.rand(100,100)
x = np.random.rand(100,1)
b = np.matmul(A,x) 
x0 = np.zeros((100,1))
maxit = 100
tol = 0.00001
x = myGMRES_SSOR(A,x0,b,tol,maxit) 
res = b - np.matmul(A,x) 
print(norm(res))
print("Solution with gmres\n", np.matmul(A,x))
print("---------------------------------------")
print("b matrix:", b)

I hope anyone could help me figure out this!!!

Dong Le
  • 179
  • 6
  • Additional info: The right-preconditioning P here is computed using symmetric successive over-relaxation with some w between 0 and 2. The first 2 functions I have is to compute this P preconditioner. In my GMRES SSOR I use w = 1.9 – Dong Le Nov 07 '20 at 21:50
  • Try replacing any `np.linalg.inv` by according `solve`s for a start. – Peter Meisrimel Nov 08 '20 at 11:03
  • I have replaced with `scipy.linalg.solve`; however, the result still remains. I believe the problem is not coming from this. – Dong Le Nov 08 '20 at 16:02

1 Answers1

1

I'm not sure where you got you "Symmetric_successive_over-relaxation" SSOR code from, but it appears to be wrong. You also seem to be assuming that A is symmetric matrix, but in your random test case it is not.

Following SSOR's Wikipedia entry, I replaced your P_SSOR function with

def P_SSOR(A,w): 
    L,D,U = splitMat(A)
    P = 2/(2-w) * (1/w*D+L)*np.linalg.inv(D)*(1/w*D+L).T
    return P

and your test matrix with

A = np.random.rand(100,100)
A = A + A.T

and your code works up to a 12 digit residual error.

Martín Maas
  • 108
  • 8