#!/usr/bin/env python
 # -*- coding: utf-8 -*-

 import numpy, random

 ## Hilfsprozeduren

 # Gibt eine zufällige hermitesche n*m-Matrix zurück.
 def random_hermite_matrix(n,m):
     def ran(): return random.random() * 10 - 5
     xs = numpy.array([ ran() for i in range(0,n*m) ]).reshape(n,m)
     return (xs + xs.transpose().conjugate()) / 2

 # Konstruiert aus dem Hauptdiagonalenvektor alphas und dem unteren
 # Nebendiagonalvektor betas die zugehörige hermitische n*n Tridiagonalmatrix.
 # Als obere Nebendiagonale wird die komplexe Konjugation der übergebenen
 # unteren Nebendiagonale verwendet.
 # Es wird in jedem Fall eine Matrix mit komplexen Einträgen zurückgeliefert,
 # auch, wenn die Eingabevektoren reell waren.
 # len(alphas) = n, len(betas) = n - 1
 def tridiag(alphas, betas):
     if len(alphas) != len(betas) + 1:
	 raise "..."
     n = len(alphas)

     matrix = numpy.zeros((n,n), dtype=complex)

     for i in range(0,n):
	 matrix[i,i] = alphas[i]
	 if i < n - 1:
	     matrix[i,i+1] = complex(betas[i]).conjugate()
	     matrix[i+1,i] = complex(betas[i])
     return matrix

 # Wie numpy.linalg.eigvalsh, nur werden die Eigenwerte als reelle float-Objekte
 # zurückgegeben, statt unnötigerweise (und Sortierung verhindernde) komplexe
 # Zahlen.
 def eigvalsh(matrix):
     return map(lambda x: x.real, numpy.linalg.eigvalsh(matrix))

 # Schreibt einige Informationen über das Spektrum auf die Standardausgabe.
 # matrix muss eine hermitesche Matrix sein.
 def dump_info(matrix, comment):
     eigenvalues = eigvalsh(matrix)
     eigenvalues.sort()

     print "[%s] Dim: %4i, kleinster: %7.2f, zw.-kl.. %7.2f, groesster: %7.2f" % \
	 (comment, matrix.shape[0], eigenvalues[0], eigenvalues[min(1,len(eigenvalues)-1)], eigenvalues[-1])


 ## Hauptprogramm

 # Dimension der zufällig erzeugten Testmatrix
 N       = 100

 # Länge, ab der wir einen Vektor als Nullvektor ansehen und das Verfahren
 # abbrechen
 epsilon = 0.01


 # Eine Zufallsmatrix als Eingabe
 H = random_hermite_matrix(N,N)
 dump_info(H, "orig")

 # Zufälliger normierter Startvektor
 v0  = numpy.array([ random.random() + 1 for i in range(0,N) ])
 v0 /= numpy.linalg.norm(v0)

 # Beginn des Lanczos-Algorithmus.
 # Initialisierung der Variablen:
 # vs:     Liste der erzeugten Orthonormalbasis des Krylowraums
 #         span { v0, H v0, H^2 v0, H^3 v0, ... }
 # alphas: Hauptdiagonalelemente in der Darstellung von H bzgl. vs
 # betas:  Untere = obere Nebendiagonalelemente (rein reell)
 vs     = [v0]
 alphas = []
 betas  = [0]

 # Iterationsschleife
 for j in range(0,N):
     w  = numpy.dot(H, vs[j])
     if j > 0:
	 w -= betas[j] * vs[j-1]
     alphas.append(numpy.dot(vs[j], w))
     w -= alphas[j] * vs[j]

     dump_info(tridiag(alphas, betas[1:len(betas)]), "%4d" % j)

     # Mehr als N Basisvektoren kann es nicht geben, wird sind fertig.
     if j == N - 1: break

     # Haben wir es geschafft, einen vom Nullvektor echt verschiedenen
     # Basisvektor zu finden?
     beta = numpy.linalg.norm(w)
     if beta < epsilon:
	 print "Norm < kleiner epsilon: %.2f, j=%d" % (beta,j)
	 break
     else:
	 betas.append(beta)
	 vs.append(w / beta)


 # Fertig; nun noch Code zur Kontrolle.
 # Aufgrund sich verstärkender Rundungsfehler werden die Kontrollnormen nicht
 # genau null sein. Insbesondere erhält man bessere Näherungen für die kleinsten
 # und größten Eigenwerte, wenn man nicht alle Iterationen durchführt, sondern
 # vorzeitig abbricht.
 H_           = tridiag(alphas,betas[1:len(betas)])
 dump_info(H_, "ende")

 # Matrix der Basisvektoren (ist eine möglicherweise rechteckige unitäre Matrix)
 trafo_matrix = numpy.concatenate(vs).reshape(len(vs),N).transpose()
 # Darstellung von H bzgl. der gefundenen Orthonormalbasis
 H__          = numpy.dot(trafo_matrix.transpose().conjugate(), numpy.dot(H, trafo_matrix))

 # Kontrollen:
 # 1) trafo_matrix muss unitär sein
 print "||U* U - I|| = %.2f" % \
     numpy.linalg.norm(\
	 numpy.dot(trafo_matrix.transpose().conjugate(), trafo_matrix) - \
	 numpy.eye(len(vs)))
 print "||U U* - I|| = %.2f" % \
     numpy.linalg.norm(\
	 numpy.dot(trafo_matrix, trafo_matrix.transpose().conjugate()) - \
	 numpy.eye(N))

 # 2) H_ vom Lanczos-Algorithmus müsste gleich H__ sein
 print "||H_ (Lanczos) - H__ (direkt transformiert)|| = %.2f" % \
     numpy.linalg.norm(H_ - H__)
 print "||H|| = %.2f, ||H_|| = %.2f, ||H__|| = %.2f" % \
     (numpy.linalg.norm(H), numpy.linalg.norm(H_), numpy.linalg.norm(H__))