#!/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
from naginterfaces.base import utils as b_utils
[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=b_utils.EngineIntNumPyType)
icola = np.empty(3*nnz, dtype=b_utils.EngineIntNumPyType)
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
)
