[問題] 其實是電腦問題，我的桌機有多強？

[wohtp (會喵喵叫的大叔)]

話說我這個用了半輩子紙筆做研究的人，現在碰到非要數值解硬爆
Ginzburg-Landau equation的狀況。數值方法什麼的其實都還可以，
但是我對我的電腦有多少能耐完全沒概念...

然後這麼low的問題拿到程式相關板去問感覺會被噓爆，所以我來
看溫馨的物理板能不能接納我。


配備是普通的i5，8Gb ram，用scipy做計算。



我只是想很簡單的解個線性方程組 Ax = b。

A 是 N*N 對稱矩陣，而且非常sparse，只有 9N 個元素不是零。
但是 N 大概是幾十萬。

然後我要一直 loop 到牛頓法收斂，所以這個要做很多次...



所以我想幹的事：

 1) 大得太誇張了，洗洗睡吧
 2) 電腦可能夠力，但是要有心理準備會很慢
 3) 都9102年了，誰家的電腦跑不動這種小問題？

是以上哪個？



--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 123.192.93.64 (臺灣)
※ 文章網址: https://www.ptt.cc/bbs/Physics/M.1575042503.A.6A2.html
推 j0958322080: 算法？lib？語言？
→ j0958322080: 算下來好像要10GB RAM才可以的樣子
不就說scipy。

Sparse matrix本身應該不到100mb，但是我真的不知道solver會吃到多大。
（剛剛少算一個零，還是很小就是了）


推 Eriri: 內存應該不夠大
→ Eriri: 我做過N等於幾萬的 我在桌電跑跑到crash
→ Eriri: 後來丟到cluster上開平行解決
→ Eriri: 所以答案是 （4） 都9102年了 學最最基本平行運算跟cluster
→ Eriri: 應用 是理論物理學家必須技能
小時候沒學好，來不及啦。
※ 編輯: wohtp (123.192.93.64 臺灣), 11/30/2019 01:13:53
→ j0958322080: 是說丟國網中心記憶體感覺是無限的XD
推 j0958322080: 線性方程組我我是用
※ 編輯: wohtp (123.192.93.64 臺灣), 11/30/2019 01:19:58
推 Eriri: matrix本身當然吃不到那麼大內存 可是不確定做的運算需要吃
→ Eriri: 多少
→ Eriri: 你還沒直接丟到電腦跑看看會不會出事嗎?
relaxation怎麼調都不穩定，所以才剛想要換方法。

有人說32Gb的內存連一萬都跑不動，也有人宣稱他在8Gb筆記本上面跑到幾百萬...
看來這真的很吃人品。

※ 編輯: wohtp (123.192.93.64 臺灣), 11/30/2019 01:50:57
推 j0958322080: N^2 = 10^10 = 10^7 K = 10^4 M = 10G 這樣吧？
→ j0958322080: 不過用LU分解的話大概是 N^2 +2N 記憶體
→ j0958322080: 如果你的A是方陣的話，回歸可以用QR或SVD解
→ j0958322080: 不過你只存非0的元素或許可以很小
→ Eriri: 他是sparse matrix
→ Eriri: 但就算是sparse, N到了十萬 運算起來大概也不會很快的
推 Eriri: 我剛在matlab隨便測了一下 用matlab內建方式直接解稀疏矩陣
→ Eriri: 線性方程組 N=十萬的時候被告知記憶體不夠 我現在的桌電記
→ Eriri: 意體是16GB 看有沒有更聰明能更善用稀疏性質的演算法囉
→ Eriri: 但理論上 matlab應該自己也會調用能善用稀疏性質的方式 所
→ Eriri: 以可能你需要一些更特別或專門的演算法
→ Eriri: 其實與其煩惱這些...直接學怎麼做平行運算跟cluster操作就
→ Eriri: 好了...
推 iHaveAPen: 如果是sparse matrix,就存不是0的元素就好。另外再多
→ iHaveAPen: 創一個矩陣存這些非0元素在原矩陣的位置，記憶體就不存
→ iHaveAPen: 在不夠的問題
→ iHaveAPen: 然後這種東西哪有吃人品....資源的upper bound絕對是
→ iHaveAPen: 算得出來的
→ Eriri: 存矩陣當然夠啊= =
→ iHaveAPen: 然後如果是用迭代找數值解，基本上要多試幾種，沒有一
→ iHaveAPen: 種方法可以保證你的數值解在收斂過程中不會爆掉
→ Eriri: 人家矩陣創出來 要拿來解方程組的 我是沒特別研究一般解線
→ Eriri: 性方程組能夠優化到甚麼程度 但就以最最無腦的解法 直接對A
→ Eriri: 求逆 一個稀疏矩陣的逆矩陣通常就不會是稀疏的 自然到了十
→ Eriri: 萬以上記憶體就會不太夠
→ iHaveAPen: 就已經是數值解了，為什麼還要用反矩陣去思考？
→ Eriri: 題外話 避免誤解 我前面matlab測試當然不是直接求逆矩陣算
→ iHaveAPen: 數值迭代Jacobian, Gauss-Sediel, successive over rel
→ iHaveAPen: axation, etc
→ Eriri: 呃 我只是要舉個很好懂了例子 說明為什麼求稀疏線性方程組
→ Eriri: 解 記憶體可能會不夠 即使矩陣式稀疏的 也沒能保證算法需要
→ iHaveAPen: Conjugate gradient 系列的方法也很多
→ Eriri: 的記憶體可以壓縮到很小 我沒特別研究最好的算法可以做到多
→ Eriri: 好就是了
→ iHaveAPen: 但其實就是有方法可以讓稀疏矩陣碟待不會出現記憶體不
→ iHaveAPen: 夠的情況
→ Eriri: 反正實際上就是matlab直接求A\b做不到 即使A跟b都是sparse
→ iHaveAPen: 像是最基本的Jacobian
→ Eriri: type 與其花時間搞懂最優算法 自己重建演算法 不如直接開平
→ Eriri: 行
→ iHaveAPen: 迭代不是只有matlab可以用好嗎...可以自己寫
→ iHaveAPen: 又不是多複雜 自己寫一個不用幾行的程式就能算了
→ iHaveAPen: 用matlab就只是偷懶而已
→ Eriri: 搞平行丟cluster也不用幾行阿...= =
→ Eriri: 這又不是我的工作 我當然偷懶 我只是在幫他測試記憶體大概
→ Eriri: 要多少而已
→ Eriri: 解決問題的方法當然不只一個 但這個問題又不是我要做的問題
→ iHaveAPen: 綜合上述 我的答案是2
→ iHaveAPen: 然後我建議自己寫一個c++ fortran來解就好
推 Eriri: 很多方法 減少內存使用 但卻可能是透過增加執行時間達成 線
→ Eriri: 性代數運算lib都是很成熟很優化過的 我相信scipy在線性代數
→ Eriri: 上已經不算太慢了 但要做到N=10萬在單機上跑 我相信不管怎
→ Eriri: 樣一定都至少很吃力要等很久 就算用C跟Fortran也是 尤其我
→ Eriri: 想 通常不會只是求到解就好 如果真的要算甚麼物理量的話 那
→ Eriri: 又常常要拿特徵解矩陣做運算 這種情況矩陣運算也未必是稀疏
→ Eriri: w大也是做凝態理論 反正我的建議是 學會平行運算跟cluster
→ Eriri: 操作 只會有好處 未來絕對還會用 而且w大非常聰明 也是名校
→ Eriri: PhD畢業 沒有道理幾天內搞不定
→ Eriri: 如果真的不想搞定 那其實也可以交給小弟 小弟"保證"可以搞
→ Eriri: 定 只是到時候文章二作記得放上小弟吧XD
→ iHaveAPen: 補充一點，我跟E大的建議其實沒衝突 一樣可以自己用for
→ iHaveAPen: tran c寫 然後平行有很多level 從純cpu的openmp, mpi,
→ iHaveAPen:  到單張gpu多張gpu 在筆電上最容易執行的是openmp 在
→ iHaveAPen: 原程式加入不多的語法就可以很容易操滿筆電的計算資源
→ rex0707: 我用一般i7解Poisson eq.到上千萬網格都沒問題啊 RAM插
→ rex0707: 24G 印象中應該很夠用
→ rex0707: 我是用PCG
推 maplefff: 我記得scipy還是numpy有bend storage的方式
→ maplefff: 你的矩陣時間上只需要9xN的size, 然後再直接調用求解函
→ maplefff: 數即可
→ maplefff: *時間->空間
→ maplefff: from scipy.linalg import solve_banded
→ maplefff: 建議可以嘗試搜尋bend storage和上面這個函數嘗試看看
→ maplefff: https://reurl.cc/VaeDWY  <- Scipy的ref
→ maplefff: https://i.imgur.com/opjC9SY.png
推 j0958322080: 他的0應該是隨機在矩陣中吧
→ wohtp: scipy的sparse matrix solver直接叫umfpack，我相信任何人
→ wohtp: 自己用fortran重寫都不會好多少
推 HeterCompute: 都2019了還有人叫人重造輪子，網路上package這麼多
→ HeterCompute: 可以用，自己寫平行能寫得贏我是不信啦
→ HeterCompute: mumps pardiso先試用這兩個，當初念書也是用這2個
→ HeterCompute: 我的答案是2019年會用工具的話就是3，瞎試會變成1
→ Eriri: 那些線性代數lib 無論穩定性跟優化程度都是非常高的 你所能
→ Eriri: 想像的到優化方法他們大概都自動加入了 甚至只會做得更好
→ Eriri: 以我的經驗 我也不相信一般人寫solver可以輕易超過那些lib
→ Eriri: 反正想花時間嘗試有沒有單機可行的調用方法就花時間 但在
→ Eriri: cluster上開平行 效能直接用滿 根本就不必煩惱這些問題...
→ Eriri: 恩 我想應該不會有人誤會= = 但預防萬一 還是澄清一下
→ Eriri: 我說的平行當然不是叫你自己用平行運算寫solver 而是除了
→ Eriri: solver以外的一般運算 有必要的話可以用簡單的平行 剩下的
→ Eriri: 在cluster上裝lib 照用這些lib的函式就好了
→ Eriri: 你這個問題 根本不困難 我相信用不到MPI
→ saltlake: 就算有公認可靠的程式庫可呼叫，自己也得花時間了解，
→ HeterCompute: 至少他這case不需要用cluster，cluster從零開始到
→ HeterCompute: 能work對新手也不容易
→ saltlake: 挑選真正適用的函式，畢竟各種演算法可能有其限制的適用
→ saltlake: 範圍，像是有些僅適用對稱矩陣。挑對之後也得弄清楚
→ saltlake: 正確的呼叫方式，如變數型態以及各個變數的意義。
→ saltlake: 搞錯變數型態(如單倍精)或者陣列資料排列格式等，這種
→ Eriri: cluster跟openmp 一兩天內就可搞定 學完之後再也不必這樣
→ Eriri: 見招拆招
→ saltlake: 看起來很荒謬的錯誤，實務上都是會發生的。
→ Eriri: 老實講 他這個case需不需要用cluster 還要看他還打算做什
→ Eriri: 麼
→ saltlake: 還有，程式庫內涵式的撰寫固然往往已經考慮各種情況
→ saltlake: 的優化，卻也因此讓程式變得複雜，進而讓瞭解正確叫用
→ saltlake: 方式變得複雜。
→ saltlake: 前面所述了解函示庫哪些個函式才適用自己問題的過程，
→ saltlake: 就表示用戶必須有有一定的數值計算理論基礎，否則可能
→ saltlake: 根本沒意識到要如何作挑選。
→ saltlake: 甚至可能不知道應該要針對問題某些特性作挑選
→ iHaveAPen: 還是要看依你的矩陣好壞來決定該用什麼方法，然後自己
→ iHaveAPen: 寫lib有自己寫的好處，因為往往不是只要答案，還要做
→ iHaveAPen: 其他事。所以也不用嘴為什麼我建議自己寫，我也沒有要
→ iHaveAPen: 嘴直接使用lib的，依人所好
→ saltlake: 用既有的函式還可能遇到該函式為了含括各種狀況而肥大
推 sunev: 要看wohtp對程式有多不熟了，無論是自己寫還是用lib
→ sunev: 都有許多「眉角」，有熟的人在旁可以問的話，可以省不少事
→ saltlake: ，而且計算速度可能會低於針對某特定問題寫的函式
→ saltlake: 另外，真要求嚴謹的話，不會"直接"呼叫別人的函式
→ saltlake: 必須經過一定的測試程序之後才用
→ HeterCompute: 自己寫的也需要測試，用別人的lib也需要測試，這應
→ HeterCompute: 該是常識
→ HeterCompute: 另外如果已經有能力自己從頭寫得好，也和別人的lib
→ HeterCompute: 效能差不多，你當然可以做，但對於初學者建議重新造
→ HeterCompute: 輪這想法我不同意，但我沒有要嘴人的意思，有能力
→ HeterCompute: 自己造輪那當然自己造輪，好處我想大家都知道
推 j0958322080: 看他急不急吧
→ wohtp: 其實我現在傾向於相信我邊界條件下錯了
→ wohtp: 然後relaxation才會壞掉。如果可以修好，我也不用搞牛頓法
→ wohtp: 了。
→ saltlake: 不管用誰的程式前都要測試沒錯，是否常識且不論，只針對
→ saltlake: 用他人程式狀況提醒，乃因取此圖者一般更傾向省事
→ Eriri: 在線性代數中 自己造輪子絕大多數都不會是好想法 這不是在
→ Eriri: 嘴 這是對任何真的在做研究的人都需了解的
→ Eriri: 這些傳統線性代數lib 都是由最專業的應用數學家跟電腦科學
→ Eriri: 家 花費大量心力認真建構而成的 上至NASA上太空 下到華爾
→ Eriri: 街跟矽谷 無數專業的科學家跟工程師 大量使用這些線性代數l
→ Eriri: ib數十年以上 其優化跟穩定程度 一般人是比不上的
→ Eriri: 當然 如果只是做個很簡單 大學生作業程度的事情 自己寫輪子
→ Eriri: 那或許有額外好處 但是越複雜越龐大 就越該減少使用自己的
→ Eriri: 輪子（以線性代數來說）
→ Eriri: 例如 某個史丹佛教授是這麼說的
→ Eriri: https://imgur.com/a/q8gpYz8
→ Eriri: 越是了解數值分析 其實就會越了解這些線性代數lib的偉大
→ Eriri: 我發現我昨天matlab的矩陣設錯了 看來這終究是能用執行時
→ Eriri: 間來換取調整記憶體使用量 我的記憶體全部都被佔滿 但終究
→ Eriri: 沒有out 只是我出去吃個飯回來都還沒跑完 要做事就先關了
→ Eriri: 如果你時間多可以等 又只要找出解就好 之後不會做其他規模
→ Eriri: 太大的事 （就像我以前玩個幾萬的矩陣就crash 是因為我還
→ Eriri: 做其他大矩陣運算 甚至二維積分 而且我的不是sparse)  就
→ Eriri: 像前面大家說的 或許還是有一般電腦就可以搞定的方法 你自
→ Eriri: 己決定吧～
推 Eriri: 畢竟我也沒有真正specificly測過解線性方程在一般電腦的極
→ Eriri: 限 我只知道一般N到了數萬以上的矩陣問題或操作 不管是用
→ Eriri: 什麼語言或lib 在單台個人電腦很容易就等相當久 如果不是
→ Eriri: 不能做的話
噓 yzfr6: 演算法  記憶體
→ dhtsai: 這篇想到之前搞翻intel的數值計算,發現intel的硬體bug
→ dhtsai: Pentium FDIV bug,intel回收一整批處理器,花$475 million
→ recorriendo: 既然友人都丟史丹佛了 我也丟一個
→ recorriendo: http://stanford.edu/class/msande318/notes.html
→ recorriendo: http://stanford.edu/group/SOL/software.html 現成
→ recorriendo: 函式 不過沒有打包成漂亮的library 將就用
→ wohtp: 回來回報一下，我好像可以用relaxation了
→ wohtp: 所以不需要解牛頓法喔耶
→ wohtp: 喵的沒人告訴我做relaxation的時候不能先把gauge定下來
噓 ostracize: