Source code for naginterfaces.library.examples.sparse.real_gen_basic_solver_ex

#!/usr/bin/env python3
"``naginterfaces.library.sparse.real_gen_basic_solver`` Python Example."

# NAG Copyright 2021.

# pylint: disable=invalid-name,too-many-locals

import numpy as np

from naginterfaces.library import sparse

[docs]def main(): """ Example for :func:`naginterfaces.library.sparse.real_gen_basic_solver`. This example solves an 8×8 nonsymmetric system of simultaneous linear equations using the bi-conjugate gradient stabilized method of order l=1, where the coefficient matrix A has a random sparsity pattern. An incomplete LU preconditioner is used. >>> main() naginterfaces.library.sparse.real_gen_basic_solver Python Example Results. Monitoring at iteration no. 1 residual no rm: 1.4059e+02 <BLANKLINE> Monitoring at iteration no. 2 residual no rm: 3.2742e+01 <BLANKLINE> Final results Number of iterations for convergence: 3 Residual norm: ... Right-hand side of termination criterion: 5.5900e-06 1-norm of matrix A: 1.1000e+01 <BLANKLINE> Solution vector Residual vector 1.0000e+00 ... 2.0000e+00 ... 3.0000e+00 ... 4.0000e+00 ... 5.0000e+00 ... 6.0000e+00 ... 7.0000e+00 ... 8.0000e+00 ... """ print( 'naginterfaces.library.sparse.real_gen_basic_solver ' 'Python Example Results.' ) # Define the matrix A in Coordinate Storage format n = 8 nnz = 24 irowa = np.empty(3*nnz, dtype=int) icola = np.empty(3*nnz, dtype=int) a = np.empty(3*nnz) irowa[:nnz] = [ 1, 1, 1, 2, 2, 2, 3, 3, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 7, 7, 8, 8, 8, ] icola[:nnz] = [ 1, 4, 8, 1, 2, 5, 3, 6, 1, 3, 4, 7, 2, 5, 7, 1, 3, 6, 3, 5, 7, 2, 6, 8, ] a[:nnz] = [ 2., -1., 1., 4., -3., 2., -7., 2., 3., -4., 5., 5., -1., 8., -3., -6., 5., 2., -5., -1., 6., -1., 2., 3., ] # Define the system's right hand side b = np.array([6., 8., -9., 46., 17., 21., 22., 34.]) # Compute incomplete LU factorization as a preconditioner lu_f = sparse.real_gen_precon_ilu(nnz, a, irowa, icola, n=n) # Initialize the solver comm = sparse.real_gen_basic_setup(n, method='BICGSTAB', precon='P', norm='1', weight='N', iterm=1, m=1, tol=1.0e-08, maxitn=20, anorm=0.0, sigmax=0.0, monit=1) # Note that the arrays b and x are overwritten # in the reverse-communication loop below. # On final exit, x will contain the solution and b the # residual vector x = np.zeros(n) irevcm = 0 while irevcm != 4: irevcm = sparse.real_gen_basic_solver(irevcm, x, b, comm) if irevcm in (-1, 1): # Compute b = A^t x resp. b = A x b = sparse.real_gen_matvec( 'T' if irevcm == -1 else 'N', lu_f.a[:nnz], lu_f.irow[:nnz], lu_f.icol[:nnz], x, check='N', ) elif irevcm == 2: # Compute preconditioning equation C b = x b = sparse.real_gen_precon_ilu_solve( 'N', lu_f.a, lu_f.irow, lu_f.icol, lu_f.ipivp, lu_f.ipivq, lu_f.istr, lu_f.idiag, x, check='N', ) elif irevcm == 3: # Monitoring step, print some statistics stats = sparse.real_gen_basic_diag(comm) print(" Monitoring at iteration no. %4d" % stats.itn) print(" residual no rm: %14.4e" % stats.stplhs) print() # Print statistics and the solution before leaving the program stats = sparse.real_gen_basic_diag(comm) print(' Final results') print(' Number of iterations for convergence: %4d' % stats.itn) print(' Residual norm: %14.4e' % stats.stplhs) print(' Right-hand side of termination criterion:%14.4e' % stats.stprhs) print(' 1-norm of matrix A: %14.4e' % stats.anorm) print() print(' Solution vector Residual vector') for i in range(n): print('%16.4e %16.4e' % (x[i], b[i]))
if __name__ == '__main__': import doctest import sys main() sys.exit( doctest.testmod( None, verbose=True, report=False, optionflags=doctest.ELLIPSIS | doctest.REPORT_NDIFF, ).failed )