# -*- coding: utf-8 -*-

import numpy as np
from . import qr

from logging import getLogger, DEBUG
logger = getLogger(__name__)
logger.setLevel(DEBUG)


class Arnoldi(object):

    def __init__(self, A, b, eps=1e-6, initialize=True, dot=np.dot):
        self.A = A
        self.dot = dot
        self.ortho = qr.MGS(eps=eps, dot=dot)
        self.ortho(b)
        self.eps = eps
        self.coefs = []
        if initialize:
            self.calc()

    def __iter__(self):
        return self.ortho.__iter__()

    def __getitem__(self, i):
        return self.ortho[i]

    def basis(self):
        return np.stack(self).T

    def calc(self):
        while True:
            Av = self.A * self.ortho[-1]
            coef = self.ortho(Av)
            logger.debug("Residual of Arnoldi iteration = {}".format(coef[-1]))
            self.coefs.append(coef)
            if coef[-1] < self.eps:
                return

    def projected_matrix(self):
        N = len(self.coefs)
        H = np.zeros((N, N))
        for i, c in enumerate(self.coefs):
            n = min(N, len(c))
            H[:n, i] = c[:n]
        return H


def arnoldi(A, b, **kwds):
    O = Arnoldi(A, b, **kwds)
    return O.projected_matrix(), O.basis()


def solve_Hessenberg(H, b):
    N = len(H)
    g = np.zeros((N, 1))
    g[0, 0] = b
    Hg = np.concatenate((H, g), axis=1)
    for i in range(N):
        Hg[i, i+1:] /= Hg[i, i]
        Hg[i, i] = 1
        if i == N-1:
            break
        Hg[i+1, i:] -= Hg[i+1, i] * Hg[i, i:]
    for i in reversed(range(1, N)):
        Hg[:i, N] -= Hg[i, N] * Hg[:i, i]
        Hg[:i, i] = 0
    return Hg[:, N]


def gmres(A, b, eps=1e-6, dot=np.dot):
    H, V = arnoldi(A, b, eps=eps, dot=dot)
    g = solve_Hessenberg(H, np.sqrt(dot(b, b)))
    return dot(V, g)

@@ -37,16 +37,54 @@
 y as np%0A
+from itertools import count as icount%0A
 from . i
@@ -188,30 +188,317 @@
 )%0A%0A%0A
-class Arnoldi(object):
+def norm(v, dot=np.dot):%0A    return np.sqrt(dot(v, v))%0A%0A%0Aclass Arnoldi(object):%0A    %22%22%22 Construct Krylov subspace (Arnoldi process)%0A%0A    Attributes%0A    -----------%0A    residual : float%0A        The residual of Arnoldi Process%0A    matrix_norm : float%0A        Approximated (projected) matrix norm of %60A%60%0A%0A    %22%22%22
 %0A%0A  
@@ -537,25 +537,8 @@
 e-6,
- initialize=True,
  dot
@@ -708,35 +708,8 @@
  %5B%5D%0A
-        if initialize:%0A    
     
@@ -717,16 +717,17 @@
    self.
+_
 calc()%0A%0A
@@ -856,21 +856,21 @@
     def 
-basis
+_calc
 (self):%0A
@@ -881,147 +881,351 @@
     
-return np.stack(self).T%0A%0A    def calc(self):%0A        while True:%0A            Av = self.A * self.ortho%5B-1%5D%0A            coef = self.ortho(
+%22%22%22 Main process of Arnoldi process %22%22%22%0A        self.residual = 1.0%0A        self.matrix_norm = 0.0%0A        for c in icount():%0A            v = self.ortho%5B-1%5D%0A            Av = self.A * v%0A            norm_Av = norm(Av, dot=self.dot)%0A            self.matrix_norm = max(self.matrix_norm, norm_Av)%0A            logger.debug(%22%7CAv%7C=%7B%7D%22.format(norm_
 Av)
+)
 %0A   
@@ -1251,59 +1251,237 @@
 ug(%22
-Residual of Arnoldi iteration = %7B%7D%22.format(coef%5B-1%5D
+(v, Av)/%7CAv%7C=%7B%7D%22.format(self.dot(v, Av) / norm_Av))%0A            coef = self.ortho(Av)%0A            self.residual *= coef%5B-1%5D%0A            logger.info(%22Arnoldi: Count=%7B%7D, Residual=%7B%7D%22%0A                        .format(c, self.residual
 ))%0A 
@@ -1530,24 +1530,29 @@
      if 
-coef%5B-1%5D
+self.residual
  %3C self.
@@ -1572,22 +1572,147 @@
         
-return
+logger.info(%22Matrix Norm=%7B%7D%22.format(self.matrix_norm))%0A                return%0A%0A    def basis(self):%0A        return np.stack(self).T
 %0A%0A    de
@@ -2540,24 +2540,22 @@
 H, n
-p.sqrt(dot(b, b)
+orm(b, dot=dot
 ))%0A 
@@ -2569,13 +2569,25 @@
 rn dot(V
+%5B:, :len(g)%5D
 , g)%0A

Revise coding style