ページ

2013年12月1日日曜日

Mathmaticaを用いた数値対角化の計算時間

Mathematicaを用いた対角化の計算速度をplotしてみました.

対角化した行列のclassは
ランダムではなくスパーズではない一般のエルミート行列(赤ドット)と
ランダムではなくスパースではない一般のユニタリー行列(青ドット).
また,行列は200桁精度で近似評価されています.

厳密評価のまま計算するとキットモット遅い.

行列の構成方法は気にしないと言う事でご勘弁.

実線は計算時間を3次の多項式でfittingしたもの.
最後に描画とfittingのプログラムを置いときます.
(もちろん計算時間のデータファイルがないと動きませんが・・・
欲しいなら上げますが,欲しいか?)




import numpy as np
import matplotlib.pyplot as plt

def f(x,para):
    return para[0]*x**3 + para[1]*x**2 + para[2]*x + para[3]
    
data = np.loadtxt("time_prec200.dat").transpose()
fig = plt.figure(figsize=(8,4))
ax0 = fig.add_subplot(1,2,1)
ax1 = fig.add_subplot(1,2,2)
colors=["red", "blue"]
for i, ax in enumerate([ax0,ax1]):
    ax.plot(data[0],data[1]/3600,'o', label='Hermitian matrix',markersize=10, color=colors[0])
    ax.plot(data[0],data[2]/3600,'o', label='Unitary matrix',markersize=10,color=colors[1])
#    ax.plot(data[0],data[3]/3600,'og')
    for j in range(1,3):
        x = np.arange(100, 10000,1)    
        para = np.polyfit(data[0],data[j],3)
        y = f(x,para)
        ax.plot(x,y/3600,'-',linewidth=3,color=colors[j-1])
    if i==1:
        ax.loglog()
        ax.set_xlim(50,10000)
    else:
        ax.set_ylim(0,5)
        ax.set_xlim(0,1000)
        ax.set_ylabel("calculation time (hours)",fontsize=16)            
    le = ax.legend(loc=2,numpoints=1)
    le.draw_frame(False)
    ax.grid()
    ax.set_xlabel("matrix dimension",fontsize=16)
sup = plt.suptitle("diagonalization time with Mathematica (precision=200)",fontsize=15)
sup.set_position((0.5,1.0))
plt.tight_layout()
plt.savefig("calc_time.png")
plt.show()