"""
Quadratic eigenvalue problem solvers.
"""
from __future__ import absolute_import
import time

import numpy as nm
import scipy.sparse as sps

from sfepy.base.base import output, get_default
from sfepy.linalg.utils import max_diff_csr
from sfepy.solvers.solvers import QuadraticEVPSolver

def standard_call(call):
    """
    Decorator handling argument preparation and timing for quadratic
    eigensolvers.
    """
    def _standard_call(self, mtx_m, mtx_d, mtx_k, n_eigs=None,
                       eigenvectors=None, status=None, conf=None, **kwargs):
        timer = Timer(start=True)

        conf = get_default(conf, self.conf)
        mtx_m = get_default(mtx_m, self.mtx_m)
        mtx_d = get_default(mtx_d, self.mtx_d)
        mtx_k = get_default(mtx_k, self.mtx_k)
        n_eigs = get_default(n_eigs, self.n_eigs)
        eigenvectors = get_default(eigenvectors, self.eigenvectors)
        status = get_default(status, self.status)

        result = call(self, mtx_m, mtx_d, mtx_k,
                      n_eigs, eigenvectors, status, conf,
                      **kwargs)

        elapsed = timer.stop()
        if status is not None:
            status['time'] = elapsed

        return result

    return _standard_call

class LQuadraticEVPSolver(QuadraticEVPSolver):
    """
    Quadratic eigenvalue problem solver based on the problem linearization.

    (w^2 M + w D + K) x = 0.
    """
    name = 'eig.qevp'

    _parameters = [
        ('method', "{'companion', 'cholesky'}", 'companion', False,
         'The linearization method.'),
        ('solver', 'dict', {'kind': 'eig.scipy', 'method': 'eig'}, False,
         """The configuration of an eigenvalue solver for
            the linearized problem (A - w B) x = 0."""),
        ('mode', "{'normal', 'inverted'}", 'normal', False,
         'Solve either A - w B (normal), or B - 1/w A (inverted).'),
        ('debug', 'bool', False, False,
         'If True, print debugging information.'),
    ]

    @standard_call
    def __call__(self, mtx_m, mtx_d, mtx_k, n_eigs=None,
                 eigenvectors=None, status=None, conf=None):
        if conf.debug:
            ssym = status['matrix_info'] = {}
            ssym['|M - M^T|'] = max_diff_csr(mtx_m, mtx_m.T)
            ssym['|D - D^T|'] = max_diff_csr(mtx_d, mtx_d.T)
            ssym['|K - K^T|'] = 
max_diff_csr(mtx_k, mtx_k.T)