以周道積分解三維光子晶體之廣義特徵值問題

全文

(1)國立臺灣師範大學數學系碩士班碩士論文. 指導教授：黃聰明博士. 論文題目. 以周道積分解三維光子晶體之廣義特徵值問題. 研究生：陳俞廷. 中華民國 103 年 7 月.

(2) 誌謝光陰似箭，歲月如梭，轉眼間三年已過。感謝母親毫無保留的付出，令我於金錢方面不虞匱乏。感謝姐姐與姐夫不時的與身在台北的我聚餐聊天，以度過異鄉的七年。感謝黃聰明教授對我緩慢的學習進度仍耐心教導，以至完成此篇論文。感謝王偉仲教授提供多方意見，以及良好的設備環境。感謝夏佑槐同學，於程式方面的指導。感謝連城與黃鴻霖同學於基本科目上的指導。並且感謝所有共同研討的台大，師大，中央的教授與同學。. i.

(3) 摘要在本文中，主要測試周道積分法 (Contour integral)[1] 的效能，測試並歸納出一些較好的使用原則。我們以此方法針對三維光子晶體的馬克斯威爾方程 (Maxwell equation) 取時諧波 (time harmonic) 後經過 K. Yee 的離散化程序 [3] 得到的廣義特徵值問題作效能的測試，觀察各參數間的交互作用，並且將其與 E. Polizzi 的 FEAST[10, 11, 12] 結合而提出了混合型的演算法。最後對這混合型的演算法與傳統的 Lanczos 作效能上的比較。文中的三維光子晶體馬克斯威爾方程包含兩種情形分別為簡易立方晶格 (SC,Simple Cube) 和面心立方晶格 (FCC,Face Centered cube)，並分為是否除去零空間 (null space free) 的兩種情況，於理論部分僅提供較簡易的簡易立方晶格的介紹，而我們的數值結果則著重於應用較廣的面心立方晶格。關鍵字廣義特徵值問題，週道積分，三維光子晶體，平行運算 Abstract In this papper, we consider to solve general eigenvalue problem for three dimensional photonic crystal by contour integral, and focus on the solver’s efficacy. At first, we take the time harmonic for three dimensional photonic crystal’s Maxwell equation , and Discrete by Yee’s scheme,then test the parameter for the sovler. We explain the implications of parameter for CIRR,and compare it with FEAST.After all, We propose a hybrid solver MLCIRR, it Combine CIRR and FEAST. Keyword General eigenvalue problem, contour integral, three dimensional photonic crystal, parallel computing. ii.

(4) Contents 1 導論. 1. 2 周道積分應用至瑞雷 -瑞茲投影特徵值解法. 2. 2.1. CIRR 的核心定理 . . . . . . . . . . . . . . . . . . . . . . . . . . .. 2. 2.2. FEAST 演算法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 4. 2.3. CIRR 各項參數的意涵與作用 . . . . . . . . . . . . . . . . . . . . .. 6. 2.3.1. M . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 7. 2.3.2. L. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 8. 2.3.3. N . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 9. 2.3.4. 積分路徑相關參數 . . . . . . . . . . . . . . . . . . . . . . .. 9. 2.4. 混合型的 CIRR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10. 3 馬克斯威爾方程與建構的矩陣. 12. 3.1. 從馬克斯威爾方程到離散出的矩陣 . . . . . . . . . . . . . . . . . . 14. 3.2. 以快速傅立葉變換得到的預處理矩陣 . . . . . . . . . . . . . . . . . 16. 3.3. 轉化為沒有零空間的問題 . . . . . . . . . . . . . . . . . . . . . . . 19. 4 周道積分法與三為光子晶體的馬克斯威爾方程. 22. 4.1. Multi-level 的引入 . . . . . . . . . . . . . . . . . . . . . . . . . . . 23. 4.2. 線性系統與平移值之關聯 . . . . . . . . . . . . . . . . . . . . . . . 24. 5 數值結果. 25. 5.1. 線性系統迭代代價關係 . . . . . . . . . . . . . . . . . . . . . . . . . 25. 5.2. 移去零空間的必要性 . . . . . . . . . . . . . . . . . . . . . . . . . . 27. 5.3. 參數 M 與 FEAST 的探討 . . . . . . . . . . . . . . . . . . . . . . . 27. 5.4. 參數 r 與 M 的交互影響 . . . . . . . . . . . . . . . . . . . . . . . . 29. 5.5. N 與 L 的比較 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30. 5.6. MLCIRR 與 Lancozs 的比較 . . . . . . . . . . . . . . . . . . . . . . 31. 6 結論. 32. Reference. 34. iii.

(5) 1. 導論光子晶體是一種週期性的介電質結構，由於形狀和介電質材料的不同，各區. 位產生了對應的頻率以及因頻率不同所產生的光子能帶，進而產生所謂的光子能帶隙。光子能隙的問題，在現今不論是研究或者實作上，都已經被廣泛的討論與應用。而如果要應用光子能帶隙的話，就必須了解光子結構能帶，為此，我們必需解出光子晶體取時諧變波的馬克斯威爾方程。在時諧變波中的馬克斯威爾方程中，可以選擇消去磁場或電場其中之一，進而將馬克斯威爾方程離散化後，轉換成一個廣義特徵值問題，於本文中所採取的方式為消去磁場的方式。然而，要解一個大型廣義特徵值問題對於不同的每個問題，其合適的方法都不盡相同，雖然能從矩陣書 (matrix pencil) 的形勢推測出可能較佳的解法，但真實情況仍受內部元素以及特徵值的分佈的影響，而並沒有一個標準的答案。二維的光子晶體問題，由於將其馬克斯威爾方程離散之後鎖產生的矩陣書還不至過於龐大，所以當解廣義特徵值問題中需解線性系統時，使用傳統的線性系統解法 (如:Jacobi, Gauss Seidel, SSOR, Generalized Conjugate Gradient method(GCG) 或者 LU 分解) 尚可勉強進行運算。然而當問題轉變為三維的時候，這些較古老的解法往往會因記憶體的限制，或者是收斂速度過慢而無法正常運作，不過黃聰明教授等人於 [6] 中比較了不同的預處理矩陣對於 Jacobi-Davidson 和 Lanczos 效果。而後於 [5] 中分析了簡易立方晶格的矩陣結構，並找出如何快以速傅立葉變換 (FFT) 為預處理矩陣 (preconditioner) 的方法。雖然這是僅限定在簡易晶格下的形式，但其解線性系統的速度得到了大幅度的躍進。又因三維光子晶體的廣義特徵值問題中含有大量的零根，所以後來又發展出了所謂的 Null Space Free 的轉換 [5]，將所有的零根移出。之後，同一團隊又完成了面心立方晶格的分解 [3, 4, 5]，而使這個以 FFT 為預處理矩陣的高速解法得以套用。由於近代多核心平行運算的發展，T. Sakurai 和 H. Sugiura 所提出的周道積分 (Contour integral) 為基礎發展出的特徵值演算法 [9, 13, 14, 15]，普遍稱為 SS 法。而這個方法於 2007 年中，又被提出一種加以應用的方式，即為將周道積分與瑞雷. 瑞茲 (Rayleigh-Ritz) 程序作結合，稱為 Contour Integral with Rayleigh-Ritz type (CIRR)[16]。而其後於 2008 年，T. Sakurai 的團隊又將其初始向量推廣為可區塊型 (block) 的版本 [13, 14]，並且於 2012 年由 E. Polizzi 談到了此方法迭代的可能性 [10, 12]。2013 年後 CIRR 快速的發展，除了 E. Polizzi 給出了自己的程式碼外，著名的涵式庫 Slpec 也將其納入了解廣義特徵值的方 1.

(6) 法之一。本文主要探討將三維光子晶體的廣義特徵值問題以 CIRR 運算時的情形，主要目的是觀察 CIRR 的各項參數間的交互作用，並歸納出如何恰當的使用這個符合現代多線程概念的演算法。CIRR 有兩個重要的特性，一是同時計算了大量並且不同平移值 (Shift) 的線性系統，二是特徵簇的存在對其會有非常不良的影響。因此黃聰明教授等人的貢獻，對於使用 CIRR 是非常重要的，除了縮短了大量解線性系統的時間外，移除大量的零根將令 CIRR 更容易圈選範圍，將不需擔心遺漏了極小的非零解和接近零根時受到特徵簇的影響。將會於第二章中介紹 CIRR 與其特性，並於第三章中介紹黃聰明教授等人如何重新建構矩陣，並以快速傅立葉變換為預處理矩陣，快速解出線性系統問題，再來於第四章統整 CIRR 與三維光子晶體的ㄧ些關聯和現象，接著於第五章呈現用以佐證我們推論的數值結果，並於第六章作出結論。. 2. 周道積分應用至瑞雷 -瑞茲投影特徵值解法在本節中將介紹將周道積分應用至瑞雷 -瑞茲投影特徵值解法 (CIRR)，這是. 一個以周道積分法 (SS method)[9, 13, 14, 15] 為基礎，並將其生成的空間，作為瑞雷 -瑞茲 (Rayleigh-Ritzl) 的投影空間 [16] 的演算法，是用於求一個目標區間的高精度特徵值與特徵向量的演算法。在此列出此方法的核心定理與概念，以及其收斂性。並歸納出一些我們測試之後對於此演算法的推論. 2.1 CIRR 的核心定理讓 (λi , xi ), 1 ≤ i ≤ n 代表一對矩陣書 (A, B) 的特徵值與特徵向量。假設複數平面上 m 個相異的特徵值 λ1 , ..., λm 在以 γ ∈ C 為中心，半徑為 ρ > 0 的周道積分路徑 Γ 中. 若有一個非零向量 v ∈ Cn ，令 ∫ 1 sk = z k (zB − A)−1 Bvdz, k = 0, 1, ..., m − 1 2πi Γ. (1). 並且令 S = [s0 , ..., sm − 1] ∈ Cn×m ，若此非零向量 v 由 x1 , ..., xn 所構成 v=. n ∑ i=1. 會有以下定理 [16] 2. αi xi .. (2).

(7) 定理 2.1. 若 λ1 , ..., λm 是不同的值，並且 αj ̸= 0 對於所有的 1 ≤ j ≤ m，則 span{s0 , ..., sm−1 } = span{x1 , ..., xm }.. 因此，若定義 Q = [q1 , ..., qm ] ∈ Cn×m 是一個從 S = [s0 , ..., sm−1 ]. 提取出的正交基底空間，那麼由 Theorem 2.1[9] 就可以經由使用瑞雷 -瑞茲程序由 (QT AQ, QT BQ) 中提取出所要求的範圍內之特徵對 (λi , xi ), 1 ≤ i ≤ m。我們使用 N 個等分點的梯形逼近 (1) 的積分得到下式: N −1 1 ∑ ωj − γ k+1 sˆk = ( ) (ωj B − A)−1 Bv, N j=0 ρ. 其中 ωj = γ + ρ exp(2πi. j+ 21 N. k = 0, 1, ..., m − 1. (3). ) 且 N 是一個正整數。在這個計算中，必須耗費大. 量的計算成本在於解 (ωj B − A)yj = Bv, j = 0, 1, ..., N − 1 的線性系統，而非常幸運的，可以由 [3, 4, 5] 中快速的解決此問題。. 接著關於此方法的收線性質，有以下兩個定理 [13, 14] 定理 2.2. 若 η ∈ R ，且 |η| = ̸ 1。若 k ∈ Z 且 (1 ≤ k < N ) ，則有下列關係式: N −1 k+1 1 ∑ θj ηk = . N j=0 θj − η 1 + ηN. 定理 2.3. 令 ηi = (λi − γ)/ρ, 1 ≤ i ≤ n, 則 sˆk =. n ∑ αi i=1. ηik˙ xi . ρ 1 + ηiN. 由定理 2.3[1, 2] 可們可以得到由 sˆk 中所提取的特徵向量與整個空間基底的 B-orthogonal 的誤差估計如下 [13, 14]。定理 2.4. 若 x1 , ..., xn 是一組 B-orthogonal 的正交基底，i.e., ⟨xi , xj ⟩B = δij . 則，   αi η k + O(|η |N +k ) (1 ≤ i ≤ m) i ρ i ⟨xi , xj ⟩B =  O(|η |−N +k ) (m + 1 ≤ i ≤ n) i. 3.

(8) 從此定理中，可以看出由 sˆk 提取出以 λi 為特徵值的特徵向量，如果 λi 是在積分路徑 Γ 外的話，將會以 |ηi | = |λi − ρ|/ρ 的指數速度退化。而區塊型的向量推導也已經於 [13, 14] 中證明，因此可以使用一組而並非單一的初始向量來作運算。我們可以將其寫成以下形式。令 V = [v1 , ..., vL ] ∈ Rn×L ，其中 v1 , ..., vL 是線性獨立的，而 L 則為區塊的寬度。則可以重新改寫 (3) 成為. N −1 1 ∑ ωj − γ k+1 ˆ ( ) Yj , Sk = N j=0 ρ. k = 0, 1, ..., M − 1. 其中 (ωj B − A)Yj = BV,. j = 0, 1, ..., N − 1. 顯然得，Sˆ = [Sˆ0 , ..., SˆM −1 ] 能容納的特徵對數成為了 m = M × L，為了避免多餘的計算，必需將已經退化的向量使用奇異值分離 (SVD)(singular value decomposition)。令 U T ΣW = Sˆ 為 Sˆ 的 SVD，且令 Σ = diag(σ1 , ..., σM ) ，而 m ˆ 代表. σj σ1. ˜ ∈ Rn×mˆ ≥ δ, 1 ≤ j ≤ m 的個數。則最後作為瑞雷 -瑞茲的提取空間 Q. ˆ 1 : m) 將由 S(:, ˆ 所構成. 下列即為瑞雷 -瑞茲版的周道積分法的演算法。 Algorithm 1 CIRR method 1: 輸入:A, B ∈ Rn×n , V ∈ Rn×L , N, M, γ, ρ, r, δ ˆ 1 , xˆ1 ), (λ ˆ 2 , xˆ2 )...(λ ˆ m , xˆm ) 2: 輸出:(λ 2πi(j+ 21 ) ))1, ≤ N −1. 3: 令 ωj = γ + ρ(exp( 4: 解 Yi = (ωi B − A) 5: 計算 Sˆk =. 1 N. ∑N −1 j=0. i≤N. BV, 1 ≤ i ≤ N. (. ωj −γ k+1 ) Yj , ρ. k = 0, 1, ..., M − 1. ˆ1 , ..., SˆM −1 ]，僅選取 |σj |/|σ1 | > δ 6: 使用 SVD 將已經退化的值分離 U T ΣW = [S for 1 ≤ j ≤ m ˆ ˆ 建構正交基 Q ˜ 7: 由 S ˜=Q ˜ T AQ, ˜ B ˜=Q ˜T BQ ˜ 8: 使用瑞雷 -瑞茲程序投影 A ˆ i , wî ) ˜ B), ˜ 1≤i≤m 9: 計算投影矩陣書 (A, ˆ 之特徵對 (λ ˜ wî , 1 ≤ i ≤ m 10: 還原並計算特徵向量 x î = Q ˆ. 2.2 FEAST 演算法 FEAST[10, 12] 是由 E. Polizzi 最早發表於 2009 年，並於 2012 年寫好了以 C 為基礎的套件 (package)，而其中與 CIRR 最大的差異在於，FEAST 開創子 4.

(9) 空間並不需要用 M 來創造虛擬的空間，意即是說，FEAST 中的 M = 1。如此一來，所有的瑞雷 -瑞茲空間的開拓，皆需倚靠初始向量來張成，這將造成一個很直觀的想法被引進。如果靠大量的初始向量來張成整個瑞雷 -瑞茲的投影空間，而很自然的在精準度不足的情況下，將這些算出來的特徵向量重新丟回程序中成為初始向量進行迭代，也因此使得整個進行的過程中將不再使用 SVD(才不會造成初始向量數的變化)。FEAST 的演算法如下:. Algorithm 2 FEAST 1: 輸入:A, B ∈ Rn×n , V ∈ Rn×L , N, γ, ρ ˆ 1 , xˆ1 ), (λ ˆ 2 , xˆ2 )...(λ ˆ L , xˆL ) 2: 輸出:(λ 2πi(j+ 21 ) ))1, ≤ N −1. 3: 令 ωj = γ + ρ(exp( 4: 解 Yi = (ωi B − A). ˆ= 5: 計算 S. 1 N. ∑N −1 j=0. i≤N. BV, 1 ≤ i ≤ N. (. ωj −γ k+1 ) Yj ρ. ˜ = SˆT AS, ˆ B ˜ = SˆT B Sˆ 6: 使用瑞雷 -瑞茲程序投影 A ˆ i , wî ) ˜ B), ˜ 1 ≤ i ≤ L 之特徵對 (λ 7: 計算投影矩陣書 (A, 8: 還原並計算特徵向量 x î = Sˆwî , 1 ≤ i ≤ L 9: 若尚未收斂至要求精度，則 V = BX，X = [ˆ x1 , ..., xˆL ] 重新回到第 4 步. 為了方便對照，在此在逼近積分時依然使用了梯型法則，而並沒有依照 E. Polizzi 原始論文所使用的高斯求積法 (gaussian legendre quadrature)，除此之外是不變的。由上一個小節中，可以明確的了解到，整個演算法的計算成本將花費在 N × L 的線性系統計算上。因此，若在 N 跟 CIRR 取相等，且要開拓出相同大小的投影空間，則 FEAST 將會足足多花費 M 倍的計算成本，這是一個非常沉重代價。所以一般而言，FEAST 所取的分割數，會相對少非常多，以 E. Polizzi 自身的範例中而言，通常是取 8 或 16 分割。經由我們實際測試發現，若以相同的計算成本運行 (即 N × L 相等，而 CIRR 的開拓子空間以 M 補足)，則 FEAST 僅計算一次時的精度，通常遠低於 CIRR，而通常需要至第 3 次左右的迭代才會相同的精度，藉由定理 2.4 也能佐證這個測試結果。. FEAST 的優點在於，比起 CIRR 的依靠 M 所擬似出的向量空間，FEAST 往往能有著更好的縝密性。所謂的縝密性，是指若在一個較大的範圍，且內部有較多的特徵對的情況下，CIRR 往往會漏掉內部較偏中心部分的特徵對 (eigenpair)，相較之下 FEAST 則較沒有這一個問題。而 CIRR 會產生這種問 5.

(10) 題，其主要因為定理 2.4 和 SVD 的程序所導致，若初始向量本身所含的該特徵向量過少則會消失。CIRR 使用 SVD 程序的原因，是因為 CIRR 並非一個迭代的演算法，若其中出現了精度非常差的數組特徵對，也無法判別其真偽。所以在同樣無法判別的情況下，雖然其相較於解大量線性系統時的時間是相對少的，但 CIRR 通常依然選擇將其捨去，用以減少計算瑞雷瑞茲投影時的計算成本。而 CIRR 無法迭代的原因非常的顯而易見，因為的初始向量並沒有那麼多的空間去容納所有的特徵向量。. 2.3 CIRR 各項參數的意涵與作用在這個小節中，將逐一的討論所有可以變動的參數，並借由我們數值結果的觀察，而有了一些推論。在此稍做回顧，CIRR 解一個廣義特徵值問題需要由使用者決定的參數，其中與解線性系統和特徵對較有關的是: • N 路徑分割數 • M 虛擬的投影空間倍數 • V 初始的向量值 (又或者是 L 初始的向量寬度) 而另外關於積分路徑的有： • γ 中心點 • ρ 半徑長 • r 長短軸比先從與線性系統較有關聯的參數進行分析。回顧初始向量寬度的參數 L 這個參數在於優化迭代中扮演著重要的角色，而且與瑞雷 -瑞茲程序所開拓的投影空間息息相關，並且也與計算線性系統的代價有直接的正相關。尚且還需選擇決定積分逼近準確度的路徑上的等分割數 N ，此參數決定了不僅僅是積分路徑的逼近程度，同時亦為計算線性系統得成本主要依據之一，因為所需解的線性系統總數為 N × L，並且於定理 2.4 可以看出其對於收斂與 M 的交互影響。而最後一個參數 M ，由於 M × L 為所開拓出之投影空間最大的維度，所以此參數決定了初步預計此空間可能涵有的特徵值數，除此之外，由於其影響的地方出現在 (3) 中的冪次上，所以當該值過大時，由於機器的位元數限制，會對整個計算造成不良的影響，導致前半段值或者視後半段值無效，而於定理 2.4 中可以看到其 6.

(11) 對收斂性的影響亦非常重大，其中與 N 有著交互的影響。因此 N, M, L 三個參數彼此之間息息相關，根據問題的特性以及選擇區塊的不同，更有著不盡相同的答案，但希望經過這邊的數據比較，可以對以後使用 CIRR 有一個明確的思考方向。. 2.3.1 M 首先先對 M 著手，發現顯然的，在 M × L 不夠大的情況下，即使將 N 擴大使其積分值逼近準確，且解了足夠量的線性系統，卻不足以找到所有的特徵值。. • M × L 應大於目標區域內特徵對數 2 倍以上，若存在特徵簇則需要更多。理論上，M × L 僅需剛好大於該區域特徵對數量即可，但實際上在 FEAST 中，此項的要求是 1.5 倍，而我們實際測試，在不迭代的情況下，需要 2 倍以上，如果該區域內的特徵對眾多，或者相鄰區域有著特徵簇，則需要更多。 • 以 M 擴展投影空間不會增加解線性系統的代價，但其相較於以 L 擴展投影空間，相對容易遺漏特徵對。對於開拓投影子空間來說，增加 M 的效果明顯比增加 L 還要低落。原因在於，即使假定了足夠多的線性系統解的組合，但卻沒有足夠的向量來組出不同的向量，意即是說將組出線性相依的系統，進而導致其實組出的向量數並沒有想要的多。 • 以 M 增大投影空間可以同時增加精度，但過量則會造成反蝕以及特徵對的遺漏。根據定理 2.4 得知，M 同時對於精確度有顯著的影響，然而這個值並不能過量，否則會反而侵蝕到精準度，甚至因為位元數的限制，使得一部份的向量組合自動倍拋棄故。 • 建議 M ≤ 8 且 N ≥ 2M ，若減半線性系統運算則需 N ≥ 4M 。在 T. Sakurai 的論文中 [8, 9] 提到 M 的參數條件最好是 N ≥ 2M ，而若要使用漢米頓矩陣的性質則為 N ≥ 4M 。因此若由 2 的冪次方選擇一個量值，則 8 通常是一個不錯的選擇。 7.

(12) 2.3.2 L 當決定了 M 之後，可以接續對與其息息相關的 L 做探討。 • L 的增長同時倍增開拓投影空間，組成的向量數，以及線性系統個數。 • L 為最後才考慮增加的參數雖然 L 同樣的可以開拓投影子空間，且其不像 M 有著侵蝕精準度的副作用，並且同時因為有相當於 N 增加時，增加向量數的功用，故其也能一定程度上的提升精準度。然而一切都如此美好的 L 卻是整個系統中最昂貴的參數，以開拓投影子空間維目的來說，M 並不會增加所需解的線性系統數量，可以說是一個幾乎免費的參數，但如果是由 L 來增加的話，則會毫不留情的使需解的線性系統量倍增。而從增加向量的角度而言，如果花一樣的代價去增加 N 而不是 L 的話，可以同時得到更精確的積分值，以及收斂性上更強的退化效果。整體上而言，L 的增加雖然全無壞處，但卻是最不想增加的數量。 • FEAST 因迭代的緣故，L 必須大於目標區域內的特徵對數 2 倍以上以 E. Polizzi 所提到的迭代觀點而言，L 則是一個非常重要的角色，其原因在於當把整個 CIRR 當作一步的迭代式時，那就相當於在使用轉向反冪法 (Shift inverse power method)，而這會出現一個重大的關鍵點，可以回憶定理 2.1[13]，其假設初始向量空間包含著積分路徑內特徵對的特徵向量，意即是說，當此初始向量純化到不包含其他向量時，將永遠只能收到這一個特徵對，即使區域內有著其他的特徵對。因此，如果要用整個 CIRR 來迭代，為了避免特徵對的消失，的 L 必須自身就比區域內的特徵向量還要多。 • 當區間內特徵值過多，應切割區間而非增加 L 當我們面臨 L 過大時只有兩個選擇，其一是讓選擇區域成為非常小的區域借此來消除 N 分割數過少的問題，但這同時也削弱了退化的速度，並且能同時解出的特徵向量就會非常少，雖然 CIRR 本身就是一個對特定區塊求解更精準的特徵對的解法。其二是保持分割數的數量，但這顯然會使得所需計算的線性系統量變得非常龐大。因此，絕大多數的情況下，除非在此區間的特徵對多到需要用 L 來增加投影子空間，但如果事實如此的話， 8.

(13) 我們會建議將整個積分路徑縮小，因為密集的特徵簇除了會造成精準度的下降外，甚至會使得 CIRR 完全無法運作。 • 除非確定搜尋範圍內特徵對數非常稀少，否則不建議低於 4。然而也不建議這個值低於 4 原因如同上面所述，若所找得向量本身不包涵這區域內特徵對的特徵向量的話，則無論如何都將找不到特徵向量，因此為了避免這幾發生的機率，還是妥當的讓此值不低於 4 是較好的，但若是對於小區域的單點針對性收斂的情況來說，以 L = 1 並將其他運算成本全部移到 N 上會是ㄧ個不錯的選擇。 2.3.3 N • N 所需的大小與積分路徑大小有關。 N 本身雖然說與選定積分路徑無關，但其本身卻又受到選定路徑的影響。如果以有限的分割逼近無限的積分，切割越密是越好的，但是如果本身積分路徑周長就非常的小，那即使分割數並不多，但所切割的間距依然是相對小的，進而逼近積分值。是故 N 的選定，會與所選取的積分路徑有關聯性，範圍越大的積分路徑，所需的 N 也就相對越大，這也是為什麼 CIRR 只能是針對小型區域求更精確的特徵對的算法的最主要原因。 • 在 M × L 足夠大的情況下，還有餘裕的計算成本應耗費在 N 而非 L。回到與演算法較相關的層面，N 本身的作用，除了讓積分值更準確外，其也影響到了定理 2.4 中的退化速度，並使得有更多的向量得以作線性組合。因此若不考慮計算線性系統的代價問題，N 跟 L 相同，是個越大越好的參數，並且其擁有的效果比 L 還要卓越。然而，即使所有的線性系統都可以平行計算，但我們並不會擁有無限的核心，因此應用時，應該視的硬體設備選取 N × L 的適當值，並依照要解決的問題決定 L 的大小後，將可接受的運算成本用來增加 N 。同時因為 N > 2M 的原則，所以如果是如前段 M 所建議，則 N 的最小建議值則為 16。 2.3.4 積分路徑相關參數 γ 與 ρ 的選定取決於我們有興趣的區間，與問題本身是否有特定的需求有關，而 r 本身大幅度的影響了積分路徑的大小，相對的影響到對於積分路徑切割逼近的準度。 9.

(14) • r 的大小將影響線性系統的穩定性，收斂精度，以及迭代時的收斂速度由於針對的問題是一個漢米頓矩陣 (Hermitian matrix)，所以特徵值皆分布於實數軸上，因此在要求精度的情況下，能將此值取近乎於 0，使得積分路徑最小化，並且所解之線性系統之平移值貼近所想要的特徵值。然而當平移值接近特徵值時，往往容易破壞矩陣的穩定性，進而造成解線性系統之迭代次數遽增的現象，故若將 CIRR 定為為尋找特徵值分布的初步判定方式，則可以讓此值稍為增大，用以遠離特徵值所在。必須注意的是，這其中存在著因為該平移值與內部所求之特徵值相距太遠，而無法找到特徵對的可能性。 • 區域的選取與實作目標結合因本文中所測試的問題為三維光子晶體的馬克斯威爾方程的光子能隙問題 (將於下一章節介紹)，為了與物理現象和實作結合，我們並不對高能階的光子能隙感到興趣，因此就得到了明確的有興趣的區域，即是低能階的區段，更棒的是由於其 5.1 小節中所介紹的特質，可以很快得找到第一次該有的停損點，進而明確的知道該選定的區域，所以在以三維光子晶體的馬克斯威爾方程的光子能隙問題為目標時，就已經解決了關於選定積分路徑的問題。. 2.4. 混合型的 CIRR. 從上一小節中，了解到了 CIRR 與 FEAST 兩個演算法的優缺點，在此我們提出一個折衷的辦法來改善 CIRR 對大型區域內的特徵向量有所遺漏的問題，卻又不花費有如 FEAST 般龐大的計算代價，由於我們並未搜尋到相關的論文，所以在此暫且稱之為 MLCIRR(Muti-Level CIRR)。由上小節得知，CIRR 由於初始向量空間無法容納下所有的特徵向量，所以無法迭代，而 FEAST 又因為龐大的計算代價和精度問題而讓人詬病。顯然得，若想要讓整個演算法得迭代得以成立，我們必須將架構偏向 FEAST，但我們可以於一些地方做出改進。. • 將 CIRR 中的 SVD 步驟移除 CIRR 的問題在於 M 會使得內部得特徵向量經過 SVD 而消逝，而又因無法迭代而不能取消 SVD。因此，若今日的架構可以迭代，則我們可以去除 10.

(15) SVD 的步驟，進而使得縝密性上升。如此一來就直接改善了 FEAST 第一次的代價和精度，即使從這步開始進行 FEAST，都已經是有所改善。 • 隨著每次的迭代改變積分路徑由於需要迭代的緣故，不得不把架構偏向 FEAST，但這不代表不能再度加快其收斂。回憶積分路徑在是一個大的橢圓曲線，而我們所求的特徵對又由於矩陣是漢米頓矩陣的緣故，所有的特徵值都在於實數軸上。若將大橢圓在進行迭代時分化成小橢圓，而依然包覆著整個實數軸，且這個舉動，比起直接再用大橢圓進行迭代，並不會增加任何的計算量，甚至會有減少，因為僅將所對應該區域的的特徵向量成為該小橢圓中的初始向量，而不在所求區域內之向量，自然就不再被進行計算。而注意到由於線性系統的計算量是 N ∗ L 所以比起 FEAST 的迭代，計算量顯然是下降的 (因為 N 可以固定，而 L 只會減少)。更棒的事情是，比起 FEAST 的大橢圓用非常少的 N 去切割，如今小橢圓一樣用 N 去切割，則相對於原本的積分值逼近程度，是有非常大幅度的提升的，因此迭代收斂速度也會遠快於原始的 FEAST，於我們的實際測試結果而言，最多 2 次就可以直接達到解線性系統的精準度，而若以原始的 FEAST 則需要約 6 次左右的迭代。 • 視殘差決定迭代所需包覆的區域關於如何選擇小橢圓，我們視第一輪 CIRR 的結果，來決定該如何選擇。如果第一次的 CIRR，有部分的特徵對殘差過高 (> δ)，則我們建議使用將整條實數軸覆蓋的方式，然後以每兩個特徵值為端點，並設置最小的距離以避免特徵簇或重根問題。假設 l 為分割的小橢圓數，而該區域內有 m ˆ l 個 CIRR 算出的特徵值，則給予此小橢圓 m ˆ l − 1 個隨機向量，並且以 2M ≤ N 的原則決定 M 。 • 殘差過高時，以隨機向量為初始值由於殘差過高的原因，此處的特徵向量本身通常並沒有太好的修正性值，而且由於不能精準的知道其到底是在兩個相鄰小橢圓中的哪一個，導致我們無法得知將其至入何者，為了避免這些煩惱，故我們採取的是隨機向量。以數值結果來說，即使並沒有使用到特徵向量來修正，卻因為積分路徑的大幅度縮小，以及內部的特徵對稀疏的原因，於此步驟就可以得到解線性系統精度的殘差，而鮮少需要用到第三步，亦即是若殘差足夠小的情況下所使用的第二種劃分法。 11.

(16) • 殘差足夠小時，以極小的區域包覆單根，並導入特徵向量為初始向量如果已有足夠小的殘差 (ε)，則可以直接以每個特徵值為圓心，用非常小的數值為半徑 (ρ = 10ε)，分別將每個特徵值包覆，以求直接達到計算目標精度 (ϵ)，若有重根或特徵簇現象，則置入相對應之特徵向量，且不再重新計算一次該區域。該注意的是，如果此廣義特徵值問題中，有重根，亦或是非常鄰近之特徵簇，那就會使得我們的半徑幾乎區近於 0，雖然並不會使整個程序無法運行，卻會造成大量的計算浪費。因此，我們需要再附加兩個條件來避開此類問題。第一個條件是，我們設定一個最小的半徑 rmin ，若兩個特徵值 λj 與 λj+1 的距離小於 rmin 則 λj+1 = λj+2 以此類推直至 |λj − λj+l | > rm in。如此一來就解決了過度浪費計算量的問題。然而若問題本身擁有特徵簇，則無論始任何的特徵問題解法都會對其感到棘手，而感謝區塊形式的解法被推導完成，我們可以藉由增加初始向量的數量 L，來達到開拓足夠的投影空間來容納其數量，雖然其也同時造成了數倍的計算量，但其是無法被避免的。以有限的核心數來說，比起 FEAST 對全部的特徵向量再做一次迭代，我們此處可以僅對不足的區域做迭代，而若核心數足夠的情況下，亦可同時改善所有的區域，故多了非常大的彈性，然而此處可能牽扯到技術層面上如何分配核心的問題，但由於我們的數值測試是用 8 核心進行測試，故只要 8 ≤ N 就已經可以達到全分配，所以並無處理此問題，但應該並非無法解決的課題。此演算法最大的變動，除了將第一輪的 FEAST 以 CIRR 取代而節省計算成本並提升精準度外，在於之後的迭代將積分路徑切為小圈使得精準度大幅提升，並其代價與原始 FEAST 的迭代相較不增反減。然而，重新劃分積分路徑會使得整個演算法與傳統的線性系統解法 LU 分解無法共用，因為若不改變積分路徑和化分數，則整個 FEAST 中都僅需解同一組線性系統問題，而當在迭代時一併改變積分路徑的話，則無法反覆的利用於第一輪解出的 LU 分解之結果，但在矩陣問題都非常龐大的情況下，通常並不會使用 LU 分解來解線性系統問題，因此這並不能算上非常大的缺憾。. 3. 馬克斯威爾方程與建構的矩陣在本節中，將介紹如何生成並離散三維光子晶體的馬克斯威爾方程，並如何. 倚靠快速傅立葉 (FFT) 變換為預條件 (preconditioner) 快速解出線性系統 12.

(17) Algorithm 3 MLCIRR method 1: 輸入:A, B ∈ Rn×n , V ∈ Rn×L , N, N2 , M, M2 , γ, ρ, ρmin , r, δ ˆ 1 , xˆ1 ), (λ ˆ 2 , xˆ2 )...(λ ˆ m , xˆm ) 2: 輸出:(λ 2πi(j+ 21 ) ))1, ≤ N −1. 3: 令 ωj = γ + ρ(exp( 4: 解 Yi = (ωi B − A) 5: 計算 Sˆk =. 1 N. ∑N −1 j=0. i≤N. BV, 1 ≤ i ≤ N. (. ωj −γ k+1 ) Yj , ρ. k = 0, 1, ..., M − 1. ˆ 建構正交基 Q ˜ 6: 由 S ˜=Q ˜ T AQ, ˜ B ˜=Q ˜T BQ ˜ 7: 使用瑞雷 -瑞茲程序投影 A ˆ i , wî ) ˜ B), ˜ 1≤i≤m 8: 計算投影矩陣書 (A, ˆ 之特徵對 (λ ˜ wî , 1 ≤ i ≤ m 9: 還原並計算特徵向量 x î = Q ˆ 10: 計算殘差 εi 11: if ∃εi > δ then 12:. i = 1, j = 1, k = 1. 13:. while i ≤ m − 1 do ˜ k −λ ˜ i+1 λ 2. 14:. k = i, ρj =. 15:. while ρi < ρmin 且 i < m − 1 do. 16:. i=i+1. 17:. ρj =. 18:. end while. 19:. γj =. 20:. end while. 21:. 以 γj , ρj 為中心和半徑，N2 , M2 , Lj 以 CIRR 解各區域 GEP。. ˜ k −λ ˜ i+1 λ 2. ˜ k +λ ˜ i+1 λ ,i 2. = i + 1, Lj = i − k. 22: end if 23: while εi > ϵ do 24:. 令 γ = λi , ρi = 10εi , Vi = xi ，以 FEAST 解該區域 GEP. 25: end while. 13.

(18) Figure 1: 左為簡易立方晶格之單一原求與圓柱垂直圓柱示意圖，右為面心立晶格之單一圓球與橢圓柱示意圖. 3.1. 從馬克斯威爾方程到離散出的矩陣. 光子晶體的馬克斯威爾方程可以寫成.    ▽ × H = iωE,      ▽ × E = −iωu H 0   ▽ · εE = 0,      ▽ · H = 0, H 和 E 分別代表磁場和電場，ω 代表頻率，u0 代表磁場常數，ε 為材質的介電系數。令 λ = u0 ω 2 並將磁場 H 以電場 E 的關係式帶入，則我們可以得到一個特徵值問題如下   ▽ × ▽ × E = λεE,  ▽ · εE = 0, 取時諧波後又由於布洛赫定理 [19]，可得 E(x + aℓ ) = ei2πk·aℓ E(x) 其中 l = 1, 2, 3 分別為三個維度方向的序號，而 aℓ 為光子晶體的晶格移動的方向向量，因此我們所要解的三維光子晶格的兩種標準型態即為當 al 的夾角 14.

(19) 為. π 2. 時的簡易晶格 (SC,Simple Cube)，與夾腳為. π 3. 的面心立方晶格 (FCC,Face. Centered Cubic)。顯然的 SC 的三個移動向量剛好分別為各方向的單位向量，而 FCC 的三個分量則分別為. [ √ ]T [ √ ]T a a 3 a 1 2 1 1 a1 = √ [1, 0, 0]T , a2 = √ , , 0 , a3 = √ , √ , 3 2 2 2 2 2 2 2 3. 而由於.    ▽ × E ==  . ∂E3 ∂y. −. ∂E2 ∂z. ∂E1 ∂z. −. ∂E3 ∂x. ∂E2 ∂x. −. ∂E1 ∂y. . .     =  . 0. − ∂∂z. ∂ ∂z. 0. − ∂∂y. ∂ ∂x. . . ∂ ∂y. E  1     − ∂∂x   E2    0 E3. (4). 因此，借由 Yee’s Scheme [23] 將 (4) 中的 ▽ × ▽× 離散為 C ∗ C(6) 且 E = e 則若使 A = C ∗ C 則可以將馬克斯威爾方程重新改寫成一個廣義特徵值問題， Ae = λBe 其中 B 為由 ε 所離散出的正定對角矩陣，若定 δ1 , δ2 , δ3 為分割後的長度，而 n1 , n2 , n3 為分割數，則可以將 C 的結構進行以下的分析   0 −C3 C2     3 ×3 C =  C3 0 −C1  ∈ C n n ,   −C2 C1 0 C1 = In2 n3 ⊗ K, C2 = In3 ⊗ K2 , C3 = K3 ,    1   K1 = δx       1  K2 =  δy       1  K3 =  δz   . 1. . −1 .. .. . . −1. −1 1. ei2πk·a1 In1. 1. −In1 .. .. ... . .. −In1. In1 −In1 In1. ei2πk·a2 J2 In1 n2.     ∈ Cn1 ×n1   . In1 n2 .. .. ...     ∈ Cn1 n2 ×n1 n2   . .. −In1 n2. In1 n2 −In1 n2 In1 n2. ei2πk·a3 J3 15.      ∈ Cn1 n2 n3 ×n1 n2 n3   . (5).

(20) 其中若是簡易立方晶格則 J2 = In1 , J3 = In1 n2 ，若是面心立方晶格則     0 e−i2πk·a1 I n21 0 e−i2πk·a2 I n32 ⊗ In1  , J3 =   J2 =  n n I 21 0 I 32 ⊗ J2 0 至此即為我們用來測試 CIRR 對於三維光子晶體的馬克斯威爾方程所離散出的廣義特徵值問題的測試對象。. 3.2. 以快速傅立葉變換得到的預處理矩陣. 上一小節中，我們已經將三維光子晶體的馬克斯威爾方程離散為一個廣義特徵值問題，其中 A = C ∗ C 為一個漢米頓矩陣 (Hermitian)，B 是一個正定矩陣 (positive difine)，因此可以使用 CIRR 對其求解。然而回憶起 CIRR 最大量的計算代價，在於解大量 Yi = (ωi B − A)−1 BV 的線性系統問題。雖然 LU 分解在於 FEAST 時可以達到重複利用的效果，但由於矩陣大小過於龐大，實際上連開始運算都會造成記憶體不足的問題，更加不用談論效率。而較傳統的迭代法諸如 Jacobi，GS，SOR，SSOR 如老牛拖車，而 MINRES，GMRES 亦差強人意。由於 (ωi B − A) 並非對稱正定，故無法使用 CG，但我們可以退而求其次選擇廣義的 CG(GCG)，其中最廣為人知的是 CGS 和 BiCGstab。而要加速 GCG 的方法，是找一個良好的預處理矩陣 M，於 GCG 中插入解 M 的線性系統的程序，其中 M 與 A 越相似則收斂越快，但我們要求 M 必須是一個能快速求解的線性系統，通常這兩件事情是難以同時發生的。然而，黃聰明教授等人研究出對於此廣義特徵值問題，一個能快速解出 (A − τ I)y = d 的方法，顯然其與我們所須的 Yi = (ωi B − A)−1 BV 相去不遠，因為 B 僅是一個正定對角線矩陣。是故當與其極為相似的 (A − τ I) 能被快速求解時，那也離有效率的解出 Yi = (ωi B − A)−1 BV 不遠了。我們在這個章節僅呈現簡易立方晶格的定理結果，和其如何串連在一起進而快速解出線性系統問題，其所有內容皆引自 [5]，而關於面心立方晶格的定理及推導請參照 [3, 4]。首先令 . Ka,m.    =   . 1. . −1 .. .. . . −1. 1 −1 1. ei2πk·a. 16.     ∈ Cm×m   .

(21) 則 SC 中 K1 =. 1 1 1 ka1 ,n1 , K2 = Ka2 ,n2 ⊗ In1 , K3 = Ka3 ,n3 ⊗ In2 ⊗ In1 δx δy δz. 定理 3.1. Ka,m 的特徵對是 (eθm,i +θa,m − 1, Da,m um, i)，其中 θm,i + θa,m =. ı2πi ı2πk · a + m m. Da,m = diag(1, eθa,m , · · · , e(m−1)θa,m ) [ ]⊤ um,i = 1 eθm,i · · · e(m−1)θm,i 其中 i = 0, · · · , m − 1 若令 Um = [um,0. ···. um,m−1 ] ∈ Cm×m. Λa,m = diag(eθm,0 +θa,m − 1. · · · eθm,m−1 +θa,m − 1),. 則 Ka,m (Da,m Um ) = (Da,m Um )Λa,m. 藉由定理 3.1 我們可以得到 K2 ((Da2 ,n2 Un2 ) ⊗ (Da1 ,n1 Un1 )) =δy−1 ((Da2 ,n2 Un2 ) ⊗ (Da1 ,n1 Un1 ))(Λa2 ,n2 ⊗ In1 ). (6). 並且 K3 ((Da3 ,n3 Un3 ) ⊗ (Da2 ,n2 Un2 ) ⊗ (Da1 ,n1 Un1 )) =δz−1 ((Da3 ,n3 Un3 ) ⊗ (Da2 ,n2 Un2 ) ⊗ (Da1 ,n1 Un1 ))(Λa3 ,n3 ⊗ In2 ⊗ In1 ). (7). 收集 K1 , K2 , K3 的特徵向量後，我們定義一個么正 (unitary) 矩陣 1 T = √ ((Da3 ,n3 Un3 ) ⊗ (Da2 ,n2 Un2 ) ⊗ (Da1 ,n1 Un1 )) n 1 = √ (Da3 ,n3 ⊗ Da2 ,n2 ⊗ Da1 ,n1 )(Un3 ⊗ Un2 ⊗ Un1 ) n. (8). 定理 3.2. [4] 簡易立方晶格中的 C1 , C2 , C3 可以由么正矩陣 T 對角化如下: C1 T = δx−1 T (In3 ⊗ In2 ⊗ Λa1 ,n1 ) ≡ T Λx. (9). C2 T = δy−1 T (In3 ⊗ Λa2 ,n2 ⊗ In1 ) ≡ T Λx. (10). C3 T = δz−1 T (Λa3 ,n3 ⊗ In2 ⊗ In1 ) ≡ T Λx. (11). 17.

(22) [ 若我們定義 G = C1⊤. C2⊤. C3⊤. ]⊤. ，則藉由 (6) 可得. A = I3 ⊗ (G∗ G) − GG∗. (12). 將其帶入 (A − τ I)y = d 整理後可得 {I3 ⊗ (G∗ G) − τ I} y = d − τ −1 GG∗ d.. (13). 於是有以下結果定理 3.3. [5] (A − τ I)y = d 等價於   Λ   x   (I3 ⊗ Λq − τ I)˜ y = I − τ −1  Λy   Λz. .  [   Λ∗x Λ∗y Λ∗z . ]   (I3 ⊗ T )∗ d ≡ d˜ . (14). y = (I3 ⊗ T )˜ y,. (15). 其中 Λq = Λ∗x Λx + Λ∗y Λy + Λ∗z Λz. (16). 將原本的 (A − τ I)y = d 藉由定理 3.3 轉換為 (17) 的好處在於 Λq , Λx , Λy , Λz 全都是對角線矩陣，主要運算花費的成本將會在計算 T ∗ p 和 T q 上，其中 p, q 為所需要計算的向量。回憶 T 的定義 (8) 其前半段的運算亦是不費力的對角線矩陣，而當寫開 Um , m = 1, 2, 3  1 1 ···    1 eθm ,1 ··· Um =   .. ..  . .  1 e(m−1)θm ,1 · · ·.  1 e. θm ,m−1. .. ..       . e(m−1)θm ,m−1. 可以發現其實為一個快速傅立葉變換的矩陣，而 Un3 ⊗ Un2 ⊗ Un1 實為三維的快速傅立葉變換矩陣，所以無論是 T ∗ p 或 T q 都能非常迅速的用固有的涵式庫算出結果，至此，得到了一個與 (A − σB) 相近且容易運算的預處理矩陣。. 18.

(23) 3.3. 轉化為沒有零空間的問題. 如同第二章所說，當使用 CIRR 時，無論是內部有特徵簇，或者積分路徑外之周圍有特徵簇，都會使 CIRR 之運作，受到不良的影響。所面對的三維光子晶體問題中，其矩陣 A 擁有 1/3 的特徵值為零根，然而所求又是低能階區塊的特徵值，顯然得，這些零根當直接將其包覆於積分路經中時，將會佔去所有的投影空間。而當為了避開零點，又因為不能漏去趨近於零的特徵對，而取一個極小值的平移時，這切零根，依舊會顯著的侵蝕我們的投影空間。是故，為了避免這些惱人的現象，若能分離出這些零根，而直接解一個滿秩 (Full rank) 的問題，將對我們有極大的幫助，因此，我們在這小節簡略的介紹，如何將原本 3n 擁有 n 的零空間的問題，轉化為 2n 的滿秩問題。. 定理 3.4. 若 . 2Λz − 3Λy   Λp =  3Λx − Λz  Λy − 2Λx.  } {  1  ⊤  , B = k = (k1 , k2 , k3 ) ̸= 0 | 0 ≤ k1 , k2 , k3 ≤ 2a . 當 k ∈ B，Λp 是一個行滿秩 (column full rank) 矩陣。引理 3.1. [ C C1⊤. C2⊤. C3⊤. [ A C1⊤. C2⊤. C3⊤. ]⊤. =0. 意即 ]⊤. =0. (17). 定理 3.5. 令 k ∈ B，定義. .  T Λx.   Q0 =  T Λy  T Λz.  1  −2  Λq ∈ C3n×n . 則 Q0 是 A 的零空間的正交基底。. 19. (18).

(24) 當我們有了 A 的零空間的政交基底後，現在設法推導出 A 的值域的基底，定義一個行滿秩矩陣 T1 = [T ⊤ , 2T ⊤ , 3T ⊤ ]⊤ 並且對 Q0 投影, 則有下列式子: .  T.   1   (I − Q0 Q∗0 )T1 =  2T  − Q0 (Λ∗q )− 2 (Λ∗x + 2Λ∗y + 3Λ∗z )   3T   ∗ ∗ T ((Λy − 2Λx )Λy + (Λz − 3Λx )Λz )     ∗ ∗ =  T ((2Λx − Λy )Λx + (2Λz − 3Λy )Λz )  Λ−1   q T ((3Λx − Λz )Λ∗x + (3Λy − 2Λz )Λ∗y ). (19). 則 [(I − Q0 Q∗0 )T1 ]∗ [(I − Q0 Q∗0 )T1 ] =T1∗ (I − Q0 Q∗0 )T1 =[(2Λz − 3Λy )∗ + (3Λx − Λz ) ∗ (3Λx − Λz )∗ + (Λy − 2Λx ) ∗ (Λy − 2Λx )]Λ−1 q =Λ∗p Λp Λ−1 q. (20). 藉由定理 3.4 則 (23) 是滿秩的。定義  T ((Λy − 2Λx )Λ∗y − (3Λx − Λz )Λ∗z )   Q1 =  T ((2Λz − 3Λy )Λ∗z − (Λy − 2Λx )Λ∗x )  T ((3Λx − 1Λz )Λ∗x − (2Λz − 3Λx )Λ∗y ).   1  ∗  (Λp Λp Λq )− 2 . . (21). 則 Q1 是一個正交基底且 G∗ Q1 = 0 而又因 A = I3 ⊗ (G∗ G) − GG∗ 且. (G∗ G)T =(C1∗ C1 + C2∗ C2 + C3∗ C3 )T =C1∗ T Λx + C2∗ T Λy + C3∗ T Λz =T (Λ∗x Λx + Λ∗y Λy + Λ∗z Λz ) (22). =T Λq. 20.

(25) 所以 AQ1 = (I3 ⊗ (G∗ G))Q1   ∗ ∗ ∗ (G G)T ((Λy − 2Λx )Λy − (3Λx − Λz )Λz )   1   =  (G∗ G)T ((2Λz − 3Λy )Λ∗z − (Λy − 2Λx )Λ∗x )  (Λ∗p Λp Λq )− 2   ∗ ∗ ∗ (G G)T ((3Λx − Λz )Λx − (2Λz − 3Λx )Λy )   T Λq ((Λy − 2Λx )Λ∗y − (3Λx − Λz )Λ∗z )   1   =  T Λq ((2Λz − 3Λy )Λ∗z − (Λy − 2Λx )Λ∗x )  (Λ∗p Λp Λq )− 2   ∗ ∗ T Λq ((3Λx − Λz )Λx − (2Λz − 3Λx )Λy ) (23). = Q1 Λq . 因此 Q1 屬於 A 的值域。接下來可將其餘部分的值域定義為     ∗ ∗ T T (2Λz − 3Λy )     1 1     ∗ ∗ −2 ∗ ∗ Q2 = C  2T  (Λp λp ) =  T (3Λx − Λz )  (Λ∗p Λp )− 2     ∗ ∗ 3T T (Λy − 2Λx ) 又因 Q2 是么正矩陣且 G∗ C ∗ = 0，則. . ∗. G)T (2Λ∗z. 3Λ∗y ). (24). . (G −   1   − ∗ ∗ ∗ AQ2 = (I3 ⊗ (G G))Q2 =  (G G)T (3Λx − Λz )  Λp 2   (G∗ G)T (Λ∗y − 2Λ∗x )   ∗ ∗ T Λq (2Λz − 3Λy )   1   − ∗ ∗ =  T Λq (3Λx − Λz )  Λp 2   ∗ ∗ T Λq (Λy − 2Λx ) ∗. (25). = Q2 Λq 由於已經將 A 分割為 3 組互不相交的垂直基底，是故有以下定理定理 3.6. 令 Q0 , Q1 , Q2 定義如 (18)(21)(24)，且定義 Q = [Q0. Q1. Q3 ]. (26). 則 Q 是一個么正矩陣，且 Q∗ AQ = diag(0, Λq , Λq ) 21. (27).

(26) 至此得到了 A 的特徵值分解定理 3.7. 令 Q1 , Q2 如同 (26),(29) 之定義，並定義 Q2 ] ≡ (I3 × T )Λr ,. Qr = [Q1 則. Λr = diag(Λq , Λq ). (28). { 1} −1 2 span B Qr Λr = {x | Ax = λBx, λ > 0}. 1. 因此若將 A = A˜ = Qr Λr Q∗r , x = x˜ = B −1 Qr Λr2 y 則原始的廣義特徵值問題 Ax = λBx 則可以改寫為. ( 1 1 ) y (Qr Λr Q∗r ) B −1 Qr Λr2 = λQr Λr2 y 1 2. (29). 1. ⇒Λr Q∗r B −1 Qr Λr2 y. = λy. 因此將原本的 3n 的廣義特徵值問題轉化成了 22 的標準特徵值問題，若要用 CIRR 解決此問題，則所需解的線性系統問題 (A − ωj B)−1 V 則變成 (A − ωj I)−1 V )−1 ( 1 1 V ⇒ Λr2 Q∗r B −1 Qr Λr2 − ωj I ( 1 1 )−1 −1 ⇒ Λr2 (Q∗r B −1 Qr − ωj Λr 2 )Λr2 V −1. 進行 GCG 迭代時，所需之矩陣向量乘積為 Q∗r B −1 Qr V − ωj Λr 2 V 因此前項能經由 F F T 快速運算，後項僅為對角線矩陣乘積，之後再經過一次的矩陣加法運算即可算出。. 4. 周道積分法與三為光子晶體的馬克斯威爾方程已經分別於第二章介紹了要探討的演算法，並於第三章給定了針對的問題，. 因此除了個別分析其特性外，將在這個章節，介紹兩件由於兩者結合後才引發的探討。其第一件衍生出的探討事件為，由於我們所解的廣義特徵值問題是來自已知的微分方程離散後的矩陣書，所以必然可以作不同密度的分割來改變所需解的廣義特徵值問題的矩陣大小。其二，因為三維光子晶體的馬克斯威爾方程所離散之矩陣書於解線性系統時，有著當平移值越高，則解線性系統所需的迭代次數呈指數成長的現象，因此配合上我們有興趣的區間是低能階的光子能隙帶的原因，將借由這個特性選擇出一個合理的區間來使用 CIRR。 22.

(27) 4.1 Multi-level 的引入在有原始的馬克斯威爾方程的情況下，我們可以先以非常低廉的代價解出較粗分割 (rough grid) 的廣義特徵值問題，進而達到初步估計特徵值分佈的目的。以本文中的 4 個情形而言，較低的每邊分割可以取到目標分割數的 1/4 甚至 1/8 亦能得到與相當不錯的分布情形。所以估計的代價比起直接對目標分割數做，矩陣縮小了 43 甚至 83 倍。而且由於矩陣相對較小，使用瑞雷 -瑞茲投影後張出的網所覆蓋的區塊相對變大，因而得到的殘差亦會較小。由於選定 CIRR 的區間僅需要得到特徵值的分佈即可，所以不同於 Multi-level 大多數的研究重點都放在如何縮放特徵向量，進行高分割與低分割之間的切換，整體程序將會變成簡易且有效率的。雖然於 E. Polizzi[1] 中提到的迭代性質來說，提取特徵向量是有其用處的，但由於並非在同一個分割上所進行的迭代，而必須經過一個非常巨量的縮放，故除非其擁有好到不可思議的性質，否則無論是何種插值法都是沒有顯著效果的，故可以不用在這方面多做探討。 Algorithm 4 Two-level grid ˜ j , 1 ≤ j ≤ m. 1. 先以 CIRR 對低分割數的離散矩陣求解，得到 λ 2. 計算 γi , ρi , 1 ≤ i ≤ j − 1，其中 γi =. ˜ j +λ ˜ j+1 λ ,ρi 2. =. ˜ j −λ ˜ j+1 λ 2. 3. 使用 CIRR 以 γi , ρi 為中心和半徑來解出目標分割數的廣義特徵值問題。淺顯易懂的三個步驟，如果粗分割數的每個特徵值都逼近目標分割的特徵值且一一對應的話，那在於目標分割時所選出的新的範圍，所有的特徵值都會在於積分路徑的邊界上，且每個獨立的積分中都僅會有至多 2 個特徵值。在這種情況下的 CIRR 幾乎可以將其殘差 (residual) 精確度，直接提升至與解線性系統時的公差 (tolerance) 相同，甚至更低。該注意的是，為了避免重根和特徵簇的問題，需要與 MLCIRR 作相同的處理，在此不再列入。. 除了以自身為解低分割數的廣義特徵值解法外，顯然的，我們可以以其他的特徵值解法來求解，分別從理想和實際情況來探討兩者的優劣。以理想情況而言，CIRR 是一個無懈可擊的演算法，如果擁有趨近於無限的核心 (core) 以及近乎於 0 的資料傳輸時間，但如果真是如此的話，那也無需使用任何的估計方法，僅需對於連續且足夠小的區塊作 CIRR 即可。以其它的方法，此處以 Lanczos 為例，其不像 CIRR 需要選定一個範圍，亦即是其一定能得到有效的資訊，而且小的廣義特徵值問題以平行處理來說，由於資料傳遞和互相等待的關係，效 23.

(28) 能並不會比較好，除非是選定一個非常適當的範圍，否則 Lanczos 在有限的核心數下，其無論是速度還是精度通常都會比 CIRR 還要好。其有一個唯一個壞處，其無法預見解目標分割數時解線性系統的各個平移所要付出的代價 (於下一個小節中提到)，但這個問題其實是可以透過對實數軸進行一次平移值的掃描來解決。. 4.2. 線性系統與平移值之關聯. 如果沒有原始的馬克斯威爾方程，而是直接得到一組廣義特徵值問題的話，雖然我們亦可將矩陣策略性的區塊縮小或使用代數子空間法轉換為巢狀矩陣後後分別解較小的區塊求解，並且由於不提取特徵向量的緣故，並不會牽扯到代價較大的特徵向量還原。但前者的誤差通常非常的大，後者則難以找出完美的巢狀矩陣與其變換矩陣。所以在這此根據問題的特性，提出一個使用較小的代價去直接估算特徵簇的辦法. 回顧 CIRR 的參數 r 代表積分路徑的長短軸比，而由於我們現在的馬克斯威爾方程所離散出的 A, B ∈ Hermitian 所以我們的特徵值 λ ∈ R ，而我們又已知當一個線性系統將其特徵值作為平移值時，會造成不好的條件數 (ill-condtion)，進而使得線性系統的的迭代次數遽增. 因此，若將長短軸的比例改變，將路徑遠離實數軸，則解線性系統所需的迭帶次數將會是穩定且偏低的。因此可以藉由這個特性，降低使用 CIRR 作為初步估計時的代價。. 根據上述的特性，可以反向操作，直接對於實數軸上的數個等分點做平移，然後解線性系統，觀察其迭代次數，進而推估該區域是否可能存在眾多的特徵值。然而這個行為並不是對所有的問題都有效，但於我們的問題卻有明顯的正相關。雖然這是一個非常粗略的估計，但除了得到了預估可能有特徵值的範圍外，還可以選擇適當且合理的區間去避開過高代價，並且準確的預估程序結束的時間。因為已經提早知道了每個平移區間所需要的迭代次數。. 對於本文所解之問題，離散之矩陣書的分割數於平移上相對應的解線性系統所需的迭代代價是相近的，所以若我們對低分割的平移值做地毯式的掃描，則可以很快的為過高的平移值選擇一個停損點，進而選取一個適當的範圍。而與上小節中最後呼應的一點，我們可以將第一個廣義特徵值問題以 Lanczos 先解到一定量的特徵值後，第二個問題以 CIRR 選取ㄧ個近乎緊黏著實數軸的積分 24.

(29) 路徑，如此一來我們必可以於此得到等同對實數軸掃描所得到的各平移值對於解線性系統所需要的跌代次數估計值。. 5. 數值結果在此章節中，我們將呈現與上一章理論中相輔相成的實驗結果，而除了一. 次全面性的掃描 60 個廣義特徵值問題，用以證明是對整個三維光子晶體的馬克斯威爾方程都能使用外。其於的數據整理，皆為僅對其中一個廣義特徵值問題作分析。我們所使用的軟體是 MATLAB2012a，工作站為 HP SL230s Gen8， CPU:HP SL230s Gen8 E5-2640 Kit *2，硬碟:HP 300GB 6G SAS 10K 2.5in QR DP ENT HDD*2，記憶體 HP 8GB 2Rx8 PC3L-10600E-9 Kit*6 為了避免機體上的差異所造成的影響，我們在此對於計算代價的計量，多以迭代次數為基準，而少以時間作為評斷的考量。首先我們將先呈現使用 CIRR 所掃出的光子能隙，用以作為 CIRR 確實能解三維光子晶體的馬克斯威爾方程所離散出的廣義特徵值問題的證明。其中簡易立方晶格的參數為: • 晶格邊長 =1 • 球核半徑 =0.345 • 圓柱半徑 =0.11 • 介電系數比 =13 面心立方晶格的參數為: • 晶格邊長 =1 • 球核半徑 =0.12 • 橢柱短軸半徑 =0.11 • 介電系數比 =13. 5.1. 線性系統迭代代價關係. 接著我們要呈現的數據結果是線性系統的迭代次數與平移值之間的關係。可以觀察到各線性系統對於平移的穩定度，並不會因為分割數的多寡而產生顯著 25.

(30) 1 0.9 0.5 0.8 0.7 0.4 0.6 0.3. 0.5 0.4. 0.2 0.3 0.2 0.1 0.1 0 G. X. M. R. G. 0 X. U. L. G. X. W. K. Figure 2: 左為簡易立方晶格，右為面心立晶格之能帶隙圖的變化，而是當其越過或正好接近特徵值時，會有一定量的提升，因此我們可以在解低分割問題時，使用貼近實數軸的積分路徑，藉此同時得知初步估計用的特徵值，以及各平移值大致需求的跌代次數，用以設定將來的停損點，當然實際作法上，為了避免錯綜的操作，我們可以將其減化為以 Lanczos 作全部的初步特徵值估計，並且作一次沒有任何其它目的，僅為了觀測平移值對線性系統的穩定度得觀察，在有限的核心數下，這是一個比較實際的作法，當然如果試探討有限核心的情況下，那麼我們將會建議全部使用 Lanczos 作為廣義特徵值問題的解法。畢竟算一個光子能隙問題，就已經可以各自獨立解 40 個廣義特徵值問題了。而不能平行的 CIRR 是完全沒有效能可言的，其效能會比非常古老的轉向反冪法 (shift inverse power method) 還要低落，我們可以由後續的介紹觀察到。所以，雖然我們最終的目的是解出光子能隙，但在討論問題時，採取的考量方式是對單一的廣義特徵值問題探討。 NSFSC. SC 600. 450 400. n=24 n=96. 500. n=24 n=96. Iteration. 350. Iteration. 300 250. 400. 300. 200. 200. 150 100. 100 50. 0. 0 0. 2. 4. 6. 8. 10. 12. 14. 16. 18. 20. 0. 2. 4. 6. 8. 10. 12. 14. 16. 18. shift. shift. Figure 3: 左為原始簡易立方晶格 (SC)，右為面移去零空間型的簡易立方晶格 (NSFSC,Null Space Free simple cube). 26. 20.

(31) NSFFCC. FCC. 900. 1200. 800 n=24 n=96. 1000 n=24 n=96. 700 600. Iteration. Iteration. 800. 600. 500 400 300. 400. 200 200. 100 0. 0 0. 5. 10. 15. 20. 25. 30. 35. 40. 45. 0. 50. 5. 10. 15. 20. 25. 30. 35. 40. 45. shift. shift. Figure 4: 左為原始面心立方晶格 (FCC)，右為移去零空間型的面心立方晶格 (NSFCC,Null Space Free Face-Centered Cube). 5.2. 移去零空間的必要性. 承 3.3 小節所談到的，由於當求靠近零的區間的特徵對時，會受到其大量零根的影響，故移去零空間的轉換事必要的，而我們在這一小節呈現是否移去零空間之差異。如 Figure 6 所示於未移除零空間得原始矩陣中，為了避免讓零根 −2. 10. −3. 10. FCC NSFFCC. −4. residual. 10. −5. 10. −6. 10. −7. 10. −8. 10. −9. 10 −10. 0. 10. 20. 30. 40. 50. shift. Figure 5: 面心立方晶格有否經過分離零空間的比較. 直接進入選曲範圍，而以 γ = 25, ρ = 24.999 為選取區域，卻依然無法密免抓取大量的零根，進而導致整體的殘差下滑，並且抓取了一些假的特徵值。. 5.3. 參數 M 與 FEAST 的探討. 這小節將呈現為什麼將三個環環相扣的參數中以 M 為出基準的原因，我們將可以觀察到其展開之虛擬空間對於節省計算成本是多麼的有威力，以及當其. 27. 50.

(32) 適量時如何加速正交基底間的收斂，並且當其過量時又是如何的適得其反。在這裡固定了其它所有的變數，僅變動 M 的大小。 1. 2. 10. 10. 0. 10. M=2 M=4 M=8 M=16. −1. 10. M=4 M=8 M=16 M=32. 0. 10. −2. 10. −2. residual. residual. 10. −3. 10. −4. 10. −4. 10. −6. 10 −5. 10. −8. 10. −6. 10. −7. −10. 10. 10. 10. 15. 20. 25. 30. 35. 40. 45. 50. 10. 15. 20. 25. eigenvalue. 30. 35. 40. 45. 50. eigenvalue. Figure 6: 左圖固定參數 ρ = 25, N = 16, L = 16, r = 2.5 右圖固定參數 ρ = 25, N = 32, L = 8, r = 2.5 在 Figure 6 可以看出，隨著 M 的增加，殘差有著顯住的下降，但當 M 到達 16 而跟 N 相等時，殘差反而往回掉了一大節，這個例子很直接的告訴我們 FEAST 的演算法中，是完完全全的浪費了 M 這個強力的參數，注意此區間內有 19 個特徵對，故 M = 1 是完全不能運行的。在此附上 FEAST 經過迭代之 0. 2. 10. 10. First refine3 refine5 refine7 refine9. 0. 10. −2. 10. −4. 10. −6. −6. 10. −8. 10. 10. −8. −10. 10. 10. −10. 10. First refine2 refine3. −4. residual. residual. 10. −2. 10. 10. −12. 15. 20. 25. 30. 35. 40. 45. 50. eigenvalue. 10. 10. 15. 20. 25. 30. 35. 40. 45. eigenvalue. Figure 7: 左圖為與上述例子相同 r=2.5 的情況，右圖則是 r=0.001 貼近實數軸的情況後的效果，在此的參數為 N = 8, L = 32, r = 2.5，注意到三組參數在只解一次的代價都是相等的。在積分路徑同樣選擇遠離實數軸時，FEAST 的收斂速度非常不理想，其約進行了 7 次迭代後才有與 CIRR 適當 M 的精準度，而右圖即使將積分路徑緊貼實數軸，亦要等到 2 次迭代後，是故第一次的迭代，CIRR 的優勢地位是幾乎可以確定的。 28. 50.

(33) 5.4. 參數 r 與 M 的交互影響. 我們從上一小節可以明顯看出長短軸比 r 對於 FEAST 的收斂速率有著顯著的影響，那麼其對於 CIRR 又有何影響，我們將在這一個小節檢視。Firure 6中已經看到了 r = 2.5 的情況，我們在此呈現 r = 0.001 的相對應數據，可以看出 −2. −3. 10. 10 M=2 M=4 M=8 M=16. −3. 10. −4. 10. M=4 M=8 M=16 M=32. −4. 10. −5. residual. residual. 10 −5. 10. −6. 10. −6. 10. −7. 10. −7. 10. −8. 10. −8. 10. −9. 10. 10. −9. 15. 20. 25. 30. 35. 40. 45. 50. 10. 55. 10. 15. 20. 25. eigenvalue. 30. 35. 40. 45. 50. eigenvalue. Figure 8: 左圖固定參數 ρ = 25, N = 16, L = 16, r = 0.001 右圖固定參數 ρ = 25, N = 32, L = 8, r = 0.001 在 r 非常小時，M 的影響顯著的下滑，這也同時佐證了上小節中 FEAST 為什麼於 r = 0.001 時收斂速度遽增。. −3. −5. 10. 10. M=8 M=16 M=32. −4. 10. −7. residual. residual. −5. 10. −6. 10. −7. −8. 10. 10. −8. 10. 10. −9. 10. 10. M=8 M=16 M=32. −6. 10. −10. 15. 20. 25. 30. 35. 40. 45. 50. eigenvalue. 10. 10. 15. 20. 25. 30. 35. 40. 45. eigenvalue. Figure 9: 左圖固定參數 ρ = 25, N = 64, L = 4, r = 2.5 右圖固定參數 ρ = 25, N = 64, L = 4, r = 0.001 接著觀察其它數據，如 Figure 9 所顯示，我們可以觀察到，在左圖還有些為得改善，但於右圖時已經是完全沒有效果。但可以注意到，此二圖中最小的 M 都已經到達了 8，而回憶於 2.3.1 小節中提到 M 的建議值得其中一項限制為 M < 8，因此不排除是因為 M 本身已經過大而造成的，因此可以觀察 Figure 10 29. 50.

(34) −5. −5. 10. 10. M=1 M=2 M=4. −6. 10. M=1 M=2 M=4. −6. 10. residual. residual. −7. 10. −8. 10. −7. 10. −8. 10. −9. 10. −9. 10. −10. 10. 10. 15. 20. 25. 30. 35. 40. 45. 50. 10. 15. 20. 25. 30. 35. 40. 45. eigenvalue. eigenvalue. Figure 10: 左圖固定參數為 ρ = 25, N = 64, L = 32, r = 2.5, 右圖為 ρ = 25, N = 64, L = 32, r = 0.001 顯然得在 r = 2.5, 1 ≤ M ≤ 4 時 M 還是有良好並且顯著的影響，而於 r = 0.001 時在於 M = 1 及 M = 2 間依然是有明確的改善段層，但到了 M = 4 時，已經開始出現了交錯的情形，而不像左圖 r = 2.5 時的三段分明。所以在此可以充分的下一個結論，即是 M 在於 r 非常小的時候，影響較為低落，且 M > 8 後改善效果不佳，甚至造成反效果。. 5.5. N 與 L 的比較. 在 2.3.2 和 2.3.3 小節中提到，N 與 L 的同時都能給予更多的向量來用 M 進行組成，而 N 能進一步使積分的逼近更準確，且有定理 2.4 證明其對於正交基退化的影響，但 N 的擴展，並不能如同 L 給予整個系統上的穩定性。換句話說，N 雖然能在一次過程中收斂出高精度的特徵對，卻往往會因為過度的退化 (在 L 不足的情況下通常伴隨著較大的 M )，導致大量的特徵對遺漏。而在 N 不足，卻有足量的 L 的情況下，則較少發生特徵對遺漏的現象，雖然其殘差精準度通常較差。Figure 11 中，兩者擁所需解的線性系統量相同，且開拓出的投影空間亦相等，可以看出 N 和 M 在於加速收斂上驚人的效果，可是卻遺漏了 2 組特徵對，而 L 足量者，雖然僅作一次的成效不盡理想，卻穩建的將 19 組特徵對抓出。所以可以得到一個準則，在找尋大範圍概估特徵對分佈時，應採用較大的 L，而當以加強小範圍內特定幾個特徵對精準度為目標時，則應該取較大的 N 來加速收斂速率。. 30. 50.

(35) −3. 10. −4. 10. residual. −5. 10. −6. 10. −7. 10. N128M32L2 N8M2L32. −8. 10. 10. 15. 20. 25. 30. 35. 40. 45. 50. eigenvalue. Figure 11: N L 比較參考圖 80. 0. 10. 60. −2. 10 40. −4. MLCIRRfirst MLCIRRrefine1 Lanczos MLCIRRrefine2. 10. residual. 20 0. −6. 10. −8. 10. −20 −40. −10. 10. −60. −12. 10 −80. 0. 5. 10. 15. 20. 25. 30. 35. 40. 45. 10. 15. 20. 25. 30. 35. 40. 45. 50. eigenvalue. 50. Figure 12: 左圖為 MLCIRR 積分路徑參考圖，* 為起始的積分路徑，x 為第一次的迭代。右圖為各方法所解出的精度. 5.6 MLCIRR 與 Lancozs 的比較綜合上面數據的結果後，我們於本小節比較 Lanczos 與 MLCIRR 的精度與計算代價。其計算成本之數據如 Table 1，而精準度表現如 Figure 12，MLCIRR 的部分的參數為 N = 16, N2 = 8, M = 4, M2 = 4, L = 8, Γ = 25, ρ = 25, ρmin = 10−2 , r = 2.5, r2 = 1, δ = 10−4 , ϵ = 10−10 。可以看到 MLCIRR 經過一次迭代後就已經有了與 Lanczos 相當的殘差精準度，而再以各特徵值為圓心並帶入特徵向量為初始向量後，可以超越 Lanczos 的精準度。詳細的數據列於下表，若讓程序停止於相當 Lancozs 精度的ㄧ次迭代，則實測中 8 線程平行運算需要花費約 1.3 倍的時間，我們可以預估在 16 線程時超越 Lanczos。. 31.

(36) 程序. Table 1: Lanczso 與 MLCIRR 各階段的時間比較 CGS 迭代 (次) CGS 時間 (秒) 迭代倍數時間倍數. Lanczos. 2914. 6214. 1. 1. MLCIRRinitial. 10626. 21438. 3.65. 3.45. MLCIRRrefine2 13542. 28661. 4.64. 4.61. MLCIRRrefine2 14572. 30767. 5.01. 4.95. MLCIRRtotal. 80997. 13.29. 13.03. 38740. 實際 8 線程 MLCIRRinitial. 實際值 (理論值) 10602. 3449. 3.65. 0.55(0.46). MLCIRRrefine1 13517. 4642. 4.64. 0.75(0.57). MLCIRRrefine2 14538. 4826. 5.01. 0.78 (0.61). MLCIRRtotal. 12973. 13.29. 2.08(1.62). 6. 38657. 結論我們逐項的探討並且實際測試了 CIRR 中各項參數的交互影響，進一步比. 較了 CIRR 和 FEAST 的異同及優缺點，最後提出一個整合兩個方法的演算法 MLCIRR。此演算法的優勢在於採用 CIRR 的觀點，以參數 M 去創造虛擬的向量空間，得到了相當於加速 FEAST 第一次收斂的效果。之後的迭代中，我們重新劃分了路徑，使得同樣分割數的條件下，更逼近實際積分值，同時內部存在的特徵對相較於原本的問題，顯得更加的分散，更加重要的是，這些特徵對將落在不易遺漏的積分路徑旁，而非容易遺漏的中央區域。並且由於同樣擁有 M 參數的關係，我們進行迭代時的向量，將不像 FEAST 一直維持在估計量的 2 倍左右，而會逐漸與目標向量完全相等，達到完全不浪費計算量的情況。又因為積分路徑急速的收縮，導致外迭代收斂速度亦會是 FEAST 的數倍。而此演算法存在的風險為，當選取範圍內擁有大量特徵簇時，它並不能像 FEAST 慢慢的修正，而會因為積分路徑的重新劃分和起始向量不足得原因，造成無法挽回的特徵對遺漏，同樣因為積分路徑的重劃關係，當問題本身可以使用 LU 分解時，此演算法並不能像 FEAST 反覆的利用已解過的 LU 分解。關於 CIRR 的建議參數部分，我們在這邊下一個總結。 • N,M,L 的取捨當不進行迭代時，相同的線性系統計算成本下，N 與 L 之間的取捨將會是 32.

(37) 以 N 為重，原因在於若 N 不足將無法一次性的得到較好的殘差，並且同時限制了 M 的大小，進一步使得殘差得精準度降低，而在有迭代的情形下，增加 L 可以有效的確保特徵對不被遺漏。 • 長短軸比 r 的影響當的積分路徑 (平移值) 遠離特徵值時，能有效的降低迭代次數，但同時也使得無論是單一次的收斂速度或者是迭代時的收斂速度大幅下滑。並且此參數在開拓投影空間的 M 和 L 越小時，對迭代時收斂速度的影響越明顯。因此，除非僅是用以估計的情況下，否則盡量取小。 • 建議參數對於一個大範圍且尚未預估特徵對數的區域，我們建議的初始測試參數是，N = 16, M = 8, L = 16 這組參數一次容納約 50 組準確的特徵對，並且由於 L = 16 在 32 線程以下，將不會因為平移值的不同，因解線性系統時間不同而造成等帶的其情況，並且 N, M 亦都在理論上的最佳範圍內。若是對於一以知的特徵對求更精確的解，則為 N = 8, M = 4, L = 1 其中初始向量為已知的特徵向量。而關於整體 CIRR 方面: • CIRR 在大範圍下有著容易遺漏特徵對的問題 • CIRR 於良好的目標區域下 (線性系統迭代次數尚未大幅攀升, 內部特徵對稀疏) 表現將會比 Lanczos 來的好 • MLCIRR 大幅度的加速了 FEAST，但於起始的圈選範圍和預估的 M, L 有更嚴苛的要求 • 沒有任何的方式能保證 CIRR 不遺漏特徵對，故 CIRR 最安全的使用方式仍侷限於殘差的縮小於參考文獻中 [2, 20] 中可以發現，無論是 T. Sakura 或者 E. Polizzi 的團隊，於近期的研究中接有著墨於如何預估某特定範圍內的特徵對。因此，如何有效的估計範圍內的特徵對仍是 CIRR 當今最重要的課題，而本文中第四節可以靠著觀察平移值與迭代次數的關係，來避開大量的特徵簇存在，可以說是非常特列的作法。若是問題本身是經由微分方程離散而來者，靠著較粗的切割計算較小的矩陣書來估計特徵對的分佈也不失為一個好的方式。 33.

(38) References [1] J. Asakura, T. Sakurai, H. Tadano, T. Ikegami and K. Kimura, A numerical method for nonlinear eigenvalue problems using contour integrals, JSIAM Letters, Vol.1,52-55 (2009). [2] Y. Futamura, H. Tadano, and T. Sakurai, Parallel stochastic estimation method of eigenvalue distribution, JSIAM Letters, Vol. 2, 127-130, (2010). [3] T.-M. Huang, H.-E. Hsieh, W.-W. Lin and W. Wang, Eigendecomposition of the discrete double-curl operator with application to fast eigensolver for three dimensional photonic crystals, SIAM J. Matrix Anal. Appl., Vol. 34, 369-391(2013) [4] T.-M. Huang, H.-E. Hsieh, W.-W. Lin, and W. Wang.Fast lanczos eigenvalue solvers for band structures of three dimensional photonic Crystals with facecentered cubic lattice(preprint) [5] T.-M. HUANG, Y.-L. Huang.W.-W. Lin, and W.-C. Wang, A Null space free Jacobi-Davidson iteration for maxwell’s operator(preprint) [6] T.-M. Huang, W.-J. Chang, Y.-L. Huang, W.-W. Lin, W.-C. Wang and W. Wang, Preconditioning bandgap eigenvalue problems in three dimensional photonic crystals simulations, Journal of Computational Physics, Vol. 229,8684-8703(2010) [7] T.-M. Huang, Y.-C. Kuo and W. Wang, Computing extremal eigenvalues for three-dimensional photonic crystals with wave vectors near the Brillouin zone center, J. Sci. Comput., Vol. 55, 529-551(2013) [8] T. Ikegami and T. Sakurai, Contour integral eigensolver for non-Hermitian systems: a Rayleigh-Ritz-type approach, RANMEP2008Taiwanese J. Math., Vol. 14, pp. 825–837 (2010). [9] T. Ikegami, T. Sakurai, U. Nagashima, A filter diagonalization for generalized eigenvalue problems based on the Sakurai-Sugiura projection method,CS-TR08-13, Tsukuba (2008) 34.

(39) [10] A. Levin, D. Zhang, E. Polizzi, FEAST fundamental framework for electronic structure calculations: reformulation and solution of the muffintin problem computer physics communications, V183, I11, 2370-2375 (2012) [11] E. Polizzi, Density-matrix-based algorithms for solving eingenvalue problems Phys. Rev. B., Vol. 79, 115112 (2009) [12] E. Polizzi,A high-performance numerical library for solving eigenvalue problems: FEAST solver user’s guide v2.1s,(2013) [13] T. Sakurai, J. Asakura, H. Tadano, T. Ikegami and K. Kimura, A method for finding zeros of polynomial equations using a contour integral based eigensolver, Proc. Symbolic Numeric Computations 2009, Kyoto. 143-147 (2009) [14] T. Sakurai, J. Asakura, H. Tadano and T. Ikegami, Error analysis for a matrix pencil of Hankel matrices with perturbed complex moments, JSIAM Letters, Vol. 1,76-79 (2009). [15] T. Sakurai, H. Sugiura, A projection method for generalized eigenvalue problems, ISE-TR-02-189, Tsukuba,(2002) [16] T. Sakurai, H. Tadano, CIRR: a Rayleigh-Ritz type method with contour integral for generalized eigenvalue problems, Proc. The First China-JapanKorea Joint Conference on Numerical Mathematics, Vol. 36,745-757 (2007). [17] T. Sakurai, H. Tadano, T. Ikegami, U. Nagashima, A parallel eigensolver using contour integration for generalized eigenvalue problems in molecular simulation, CS-TR-08-14, Tsukuba(2008) [18] K. Senzaki, H. Tadano, T. Sakurai and Z. Bai, A method for profiling the distribution of eigenvalues using the AS method,RANMEP 2008 Taiwanese J. Math, Vol. 14,839-853 (2010). [19] M. Suzuki, I. Suzuk, Bloch theorem and Energy band, Lecture Note on Solid State Physics(2006) [20] P.-T. Tang, E. Polizzi,FEAST as a subspace iteration eigensolver Aacelerated by approximate spectral projection,SIMAX, (accepted) 35.