|
|
|
|
@ -16,7 +16,6 @@ satisfy the constraints for the 'conjugate gradient' algorithm
|
|
|
|
|
|
|
|
|
|
from typing import Dict, Any
|
|
|
|
|
import time
|
|
|
|
|
import logging
|
|
|
|
|
|
|
|
|
|
import numpy
|
|
|
|
|
from numpy.linalg import norm
|
|
|
|
|
@ -28,11 +27,6 @@ import fdfd_tools.solvers
|
|
|
|
|
from . import ops
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
__author__ = 'Jan Petykiewicz'
|
|
|
|
|
|
|
|
|
|
logger = logging.getLogger(__name__)
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
class CSRMatrix(object):
|
|
|
|
|
"""
|
|
|
|
|
Matrix stored in Compressed Sparse Row format, in GPU RAM.
|
|
|
|
|
@ -55,6 +49,7 @@ def cg(A: 'scipy.sparse.csr_matrix',
|
|
|
|
|
err_threshold: float = 1e-6,
|
|
|
|
|
context: pyopencl.Context = None,
|
|
|
|
|
queue: pyopencl.CommandQueue = None,
|
|
|
|
|
verbose: bool = False,
|
|
|
|
|
) -> numpy.ndarray:
|
|
|
|
|
"""
|
|
|
|
|
General conjugate-gradient solver for sparse matrices, where A @ x = b.
|
|
|
|
|
@ -65,6 +60,7 @@ def cg(A: 'scipy.sparse.csr_matrix',
|
|
|
|
|
:param err_threshold: Error threshold for successful solve, relative to norm(b)
|
|
|
|
|
:param context: PyOpenCL context. Will be created if not given.
|
|
|
|
|
:param queue: PyOpenCL command queue. Will be created if not given.
|
|
|
|
|
:param verbose: Whether to print statistics to screen.
|
|
|
|
|
:return: Solution vector x; returned even if solve doesn't converge.
|
|
|
|
|
"""
|
|
|
|
|
|
|
|
|
|
@ -106,11 +102,13 @@ def cg(A: 'scipy.sparse.csr_matrix',
|
|
|
|
|
|
|
|
|
|
_, err2 = rhoerr_step(r, [])
|
|
|
|
|
b_norm = numpy.sqrt(err2)
|
|
|
|
|
logging.debug('b_norm check: ', b_norm)
|
|
|
|
|
if verbose:
|
|
|
|
|
print('b_norm check: ', b_norm)
|
|
|
|
|
|
|
|
|
|
success = False
|
|
|
|
|
for k in range(max_iters):
|
|
|
|
|
logging.debug('[{:06d}] rho {:.4} alpha {:4.4}'.format(k, rho, alpha))
|
|
|
|
|
if verbose:
|
|
|
|
|
print('[{:06d}] rho {:.4} alpha {:4.4}'.format(k, rho, alpha), end=' ')
|
|
|
|
|
|
|
|
|
|
rho_prev = rho
|
|
|
|
|
e = xr_step(x, p, r, v, alpha, [])
|
|
|
|
|
@ -118,7 +116,8 @@ def cg(A: 'scipy.sparse.csr_matrix',
|
|
|
|
|
|
|
|
|
|
errs += [numpy.sqrt(err2) / b_norm]
|
|
|
|
|
|
|
|
|
|
logging.debug('err {}'.format(errs[-1]))
|
|
|
|
|
if verbose:
|
|
|
|
|
print('err', errs[-1])
|
|
|
|
|
|
|
|
|
|
if errs[-1] < err_threshold:
|
|
|
|
|
success = True
|
|
|
|
|
@ -129,7 +128,7 @@ def cg(A: 'scipy.sparse.csr_matrix',
|
|
|
|
|
alpha = rho / dot(p, v, e)
|
|
|
|
|
|
|
|
|
|
if verbose and k % 1000 == 0:
|
|
|
|
|
logging.info('iteration {}'.format(k))
|
|
|
|
|
print(k)
|
|
|
|
|
|
|
|
|
|
'''
|
|
|
|
|
Done solving
|
|
|
|
|
@ -138,16 +137,17 @@ def cg(A: 'scipy.sparse.csr_matrix',
|
|
|
|
|
|
|
|
|
|
x = x.get()
|
|
|
|
|
|
|
|
|
|
if success:
|
|
|
|
|
logging.info('Solve success')
|
|
|
|
|
else:
|
|
|
|
|
logging.warning('Solve failure')
|
|
|
|
|
logging.info('{} iterations in {} sec: {} iterations/sec \
|
|
|
|
|
'.format(k, time_elapsed, k / time_elapsed))
|
|
|
|
|
logging.debug('final error {}'.format(errs[-1]))
|
|
|
|
|
logging.debug('overhead {} sec'.format(start_time2 - start_time))
|
|
|
|
|
if verbose:
|
|
|
|
|
if success:
|
|
|
|
|
print('Success', end='')
|
|
|
|
|
else:
|
|
|
|
|
print('Failure', end=', ')
|
|
|
|
|
print(', {} iterations in {} sec: {} iterations/sec \
|
|
|
|
|
'.format(k, time_elapsed, k / time_elapsed))
|
|
|
|
|
print('final error', errs[-1])
|
|
|
|
|
print('overhead {} sec'.format(start_time2 - start_time))
|
|
|
|
|
|
|
|
|
|
logging.info('Final residual: {}'.format(norm(A @ x - b) / norm(b)))
|
|
|
|
|
print('Final residual:', norm(A @ x - b) / norm(b))
|
|
|
|
|
return x
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
jan