作者wohtp (會喵喵叫的大叔)
看板Physics
標題[問題] 其實是電腦問題，我的桌機有多強？
時間Fri Nov 29 23:48:20 2019
話說我這個用了半輩子紙筆做研究的人，現在碰到非要數值解硬爆
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://webptt.com/m.aspx?n=bbs/Physics/M.1575042503.A.6A2.html
1F：推 j0958322080: 算法？lib？語言？ 11/30 00:43
2F：→ j0958322080: 算下來好像要10GB RAM才可以的樣子 11/30 00:45
不就說scipy。
Sparse matrix本身應該不到100mb，但是我真的不知道solver會吃到多大。
（剛剛少算一個零，還是很小就是了）
3F：推 Eriri: 內存應該不夠大 11/30 00:54
4F：→ Eriri: 我做過N等於幾萬的 我在桌電跑跑到crash 11/30 00:55
5F：→ Eriri: 後來丟到cluster上開平行解決 11/30 00:58
6F：→ Eriri: 所以答案是 （4） 都9102年了 學最最基本平行運算跟cluster 11/30 01:05
7F：→ Eriri: 應用 是理論物理學家必須技能 11/30 01:05
小時候沒學好，來不及啦。
※ 編輯: wohtp (123.192.93.64 臺灣), 11/30/2019 01:13:53
8F：→ j0958322080: 是說丟國網中心記憶體感覺是無限的XD 11/30 01:14
9F：推 j0958322080: 線性方程組我我是用 11/30 01:16
※ 編輯: wohtp (123.192.93.64 臺灣), 11/30/2019 01:19:58
10F：推 Eriri: matrix本身當然吃不到那麼大內存 可是不確定做的運算需要吃 11/30 01:20
11F：→ Eriri: 多少 11/30 01:20
12F：→ Eriri: 你還沒直接丟到電腦跑看看會不會出事嗎? 11/30 01:21
relaxation怎麼調都不穩定，所以才剛想要換方法。
有人說32Gb的內存連一萬都跑不動，也有人宣稱他在8Gb筆記本上面跑到幾百萬...
看來這真的很吃人品。
※ 編輯: wohtp (123.192.93.64 臺灣), 11/30/2019 01:50:57
13F：推 j0958322080: N^2 = 10^10 = 10^7 K = 10^4 M = 10G 這樣吧？ 11/30 01:29
14F：→ j0958322080: 不過用LU分解的話大概是 N^2 +2N 記憶體 11/30 01:30
15F：→ j0958322080: 如果你的A是方陣的話，回歸可以用QR或SVD解 11/30 01:31
16F：→ j0958322080: 不過你只存非0的元素或許可以很小 11/30 01:32
17F：→ Eriri: 他是sparse matrix 11/30 01:34
18F：→ Eriri: 但就算是sparse, N到了十萬 運算起來大概也不會很快的 11/30 01:36
19F：推 Eriri: 我剛在matlab隨便測了一下 用matlab內建方式直接解稀疏矩陣 11/30 02:36
20F：→ Eriri: 線性方程組 N=十萬的時候被告知記憶體不夠 我現在的桌電記 11/30 02:37
21F：→ Eriri: 意體是16GB 看有沒有更聰明能更善用稀疏性質的演算法囉 11/30 02:38
22F：→ Eriri: 但理論上 matlab應該自己也會調用能善用稀疏性質的方式 所 11/30 02:44
23F：→ Eriri: 以可能你需要一些更特別或專門的演算法 11/30 02:44
24F：→ Eriri: 其實與其煩惱這些...直接學怎麼做平行運算跟cluster操作就 11/30 02:46
25F：→ Eriri: 好了... 11/30 02:47
26F：推 iHaveAPen: 如果是sparse matrix,就存不是0的元素就好。另外再多 11/30 04:34
27F：→ iHaveAPen: 創一個矩陣存這些非0元素在原矩陣的位置，記憶體就不存 11/30 04:34
28F：→ iHaveAPen: 在不夠的問題 11/30 04:34
29F：→ iHaveAPen: 然後這種東西哪有吃人品....資源的upper bound絕對是 11/30 04:35
30F：→ iHaveAPen: 算得出來的 11/30 04:35
31F：→ Eriri: 存矩陣當然夠啊= = 11/30 04:36
32F：→ iHaveAPen: 然後如果是用迭代找數值解，基本上要多試幾種，沒有一 11/30 04:37
33F：→ iHaveAPen: 種方法可以保證你的數值解在收斂過程中不會爆掉 11/30 04:37
34F：→ Eriri: 人家矩陣創出來 要拿來解方程組的 我是沒特別研究一般解線 11/30 04:39
35F：→ Eriri: 性方程組能夠優化到甚麼程度 但就以最最無腦的解法 直接對A 11/30 04:41
36F：→ Eriri: 求逆 一個稀疏矩陣的逆矩陣通常就不會是稀疏的 自然到了十 11/30 04:42
37F：→ Eriri: 萬以上記憶體就會不太夠 11/30 04:42
38F：→ iHaveAPen: 就已經是數值解了，為什麼還要用反矩陣去思考？ 11/30 04:42
39F：→ Eriri: 題外話 避免誤解 我前面matlab測試當然不是直接求逆矩陣算 11/30 04:43
40F：→ iHaveAPen: 數值迭代Jacobian, Gauss-Sediel, successive over rel 11/30 04:44
41F：→ iHaveAPen: axation, etc 11/30 04:44
42F：→ Eriri: 呃 我只是要舉個很好懂了例子 說明為什麼求稀疏線性方程組 11/30 04:44
43F：→ Eriri: 解 記憶體可能會不夠 即使矩陣式稀疏的 也沒能保證算法需要 11/30 04:45
44F：→ iHaveAPen: Conjugate gradient 系列的方法也很多 11/30 04:45
45F：→ Eriri: 的記憶體可以壓縮到很小 我沒特別研究最好的算法可以做到多 11/30 04:46
46F：→ Eriri: 好就是了 11/30 04:46
47F：→ iHaveAPen: 但其實就是有方法可以讓稀疏矩陣碟待不會出現記憶體不 11/30 04:46
48F：→ iHaveAPen: 夠的情況 11/30 04:46
49F：→ Eriri: 反正實際上就是matlab直接求A\b做不到 即使A跟b都是sparse 11/30 04:47
50F：→ iHaveAPen: 像是最基本的Jacobian 11/30 04:47
51F：→ Eriri: type 與其花時間搞懂最優算法 自己重建演算法 不如直接開平 11/30 04:47
52F：→ Eriri: 行 11/30 04:47
53F：→ iHaveAPen: 迭代不是只有matlab可以用好嗎...可以自己寫 11/30 04:48
54F：→ iHaveAPen: 又不是多複雜 自己寫一個不用幾行的程式就能算了 11/30 04:49
55F：→ iHaveAPen: 用matlab就只是偷懶而已 11/30 04:50
56F：→ Eriri: 搞平行丟cluster也不用幾行阿...= = 11/30 04:50
57F：→ Eriri: 這又不是我的工作 我當然偷懶 我只是在幫他測試記憶體大概 11/30 04:51
58F：→ Eriri: 要多少而已 11/30 04:51
59F：→ Eriri: 解決問題的方法當然不只一個 但這個問題又不是我要做的問題 11/30 04:51
60F：→ iHaveAPen: 綜合上述 我的答案是2 11/30 04:54
61F：→ iHaveAPen: 然後我建議自己寫一個c++ fortran來解就好 11/30 04:56
62F：推 Eriri: 很多方法 減少內存使用 但卻可能是透過增加執行時間達成 線 11/30 05:09
63F：→ Eriri: 性代數運算lib都是很成熟很優化過的 我相信scipy在線性代數 11/30 05:10
64F：→ Eriri: 上已經不算太慢了 但要做到N=10萬在單機上跑 我相信不管怎 11/30 05:11
65F：→ Eriri: 樣一定都至少很吃力要等很久 就算用C跟Fortran也是 尤其我 11/30 05:13
66F：→ Eriri: 想 通常不會只是求到解就好 如果真的要算甚麼物理量的話 那 11/30 05:15
67F：→ Eriri: 又常常要拿特徵解矩陣做運算 這種情況矩陣運算也未必是稀疏 11/30 05:17
68F：→ Eriri: w大也是做凝態理論 反正我的建議是 學會平行運算跟cluster 11/30 05:18
69F：→ Eriri: 操作 只會有好處 未來絕對還會用 而且w大非常聰明 也是名校 11/30 05:19
70F：→ Eriri: PhD畢業 沒有道理幾天內搞不定 11/30 05:20
71F：→ Eriri: 如果真的不想搞定 那其實也可以交給小弟 小弟"保證"可以搞 11/30 05:21
72F：→ Eriri: 定 只是到時候文章二作記得放上小弟吧XD 11/30 05:23
73F：→ iHaveAPen: 補充一點，我跟E大的建議其實沒衝突 一樣可以自己用for 11/30 05:28
74F：→ iHaveAPen: tran c寫 然後平行有很多level 從純cpu的openmp, mpi, 11/30 05:28
75F：→ iHaveAPen: 到單張gpu多張gpu 在筆電上最容易執行的是openmp 在 11/30 05:28
76F：→ iHaveAPen: 原程式加入不多的語法就可以很容易操滿筆電的計算資源 11/30 05:28
77F：→ rex0707: 我用一般i7解Poisson eq.到上千萬網格都沒問題啊 RAM插 11/30 09:29
78F：→ rex0707: 24G 印象中應該很夠用 11/30 09:29
79F：→ rex0707: 我是用PCG 11/30 09:30
80F：推 maplefff: 我記得scipy還是numpy有bend storage的方式 11/30 11:44
81F：→ maplefff: 你的矩陣時間上只需要9xN的size, 然後再直接調用求解函 11/30 11:44
82F：→ maplefff: 數即可 11/30 11:44
83F：→ maplefff: *時間->空間 11/30 11:44
84F：→ maplefff: from scipy.linalg import solve_banded 11/30 11:49
85F：→ maplefff: 建議可以嘗試搜尋bend storage和上面這個函數嘗試看看 11/30 11:50
88F：推 j0958322080: 他的0應該是隨機在矩陣中吧 11/30 12:18
89F：→ wohtp: scipy的sparse matrix solver直接叫umfpack，我相信任何人 11/30 14:49
90F：→ wohtp: 自己用fortran重寫都不會好多少 11/30 14:49
91F：推 HeterCompute: 都2019了還有人叫人重造輪子，網路上package這麼多 11/30 15:57
92F：→ HeterCompute: 可以用，自己寫平行能寫得贏我是不信啦 11/30 15:57
93F：→ HeterCompute: mumps pardiso先試用這兩個，當初念書也是用這2個 11/30 15:59
94F：→ HeterCompute: 我的答案是2019年會用工具的話就是3，瞎試會變成1 11/30 16:00
95F：→ Eriri: 那些線性代數lib 無論穩定性跟優化程度都是非常高的 你所能 11/30 16:03
96F：→ Eriri: 想像的到優化方法他們大概都自動加入了 甚至只會做得更好 11/30 16:04
97F：→ Eriri: 以我的經驗 我也不相信一般人寫solver可以輕易超過那些lib 11/30 16:05
98F：→ Eriri: 反正想花時間嘗試有沒有單機可行的調用方法就花時間 但在 11/30 16:08
99F：→ Eriri: cluster上開平行 效能直接用滿 根本就不必煩惱這些問題... 11/30 16:09
100F：→ Eriri: 恩 我想應該不會有人誤會= = 但預防萬一 還是澄清一下 11/30 16:16
101F：→ Eriri: 我說的平行當然不是叫你自己用平行運算寫solver 而是除了 11/30 16:17
102F：→ Eriri: solver以外的一般運算 有必要的話可以用簡單的平行 剩下的 11/30 16:17
103F：→ Eriri: 在cluster上裝lib 照用這些lib的函式就好了 11/30 16:18
104F：→ Eriri: 你這個問題 根本不困難 我相信用不到MPI 11/30 16:19
105F：→ saltlake: 就算有公認可靠的程式庫可呼叫，自己也得花時間了解， 11/30 17:02
106F：→ HeterCompute: 至少他這case不需要用cluster，cluster從零開始到 11/30 17:02
107F：→ HeterCompute: 能work對新手也不容易 11/30 17:02
108F：→ saltlake: 挑選真正適用的函式，畢竟各種演算法可能有其限制的適用 11/30 17:03
109F：→ saltlake: 範圍，像是有些僅適用對稱矩陣。挑對之後也得弄清楚 11/30 17:03
110F：→ saltlake: 正確的呼叫方式，如變數型態以及各個變數的意義。 11/30 17:04
111F：→ saltlake: 搞錯變數型態(如單倍精)或者陣列資料排列格式等，這種 11/30 17:05
112F：→ Eriri: cluster跟openmp 一兩天內就可搞定 學完之後再也不必這樣 11/30 17:06
113F：→ Eriri: 見招拆招 11/30 17:06
114F：→ saltlake: 看起來很荒謬的錯誤，實務上都是會發生的。 11/30 17:06
115F：→ Eriri: 老實講 他這個case需不需要用cluster 還要看他還打算做什 11/30 17:07
116F：→ Eriri: 麼 11/30 17:07
117F：→ saltlake: 還有，程式庫內涵式的撰寫固然往往已經考慮各種情況 11/30 17:08
118F：→ saltlake: 的優化，卻也因此讓程式變得複雜，進而讓瞭解正確叫用 11/30 17:09
119F：→ saltlake: 方式變得複雜。 11/30 17:09
120F：→ saltlake: 前面所述了解函示庫哪些個函式才適用自己問題的過程， 11/30 17:11
121F：→ saltlake: 就表示用戶必須有有一定的數值計算理論基礎，否則可能 11/30 17:11
122F：→ saltlake: 根本沒意識到要如何作挑選。 11/30 17:12
123F：→ saltlake: 甚至可能不知道應該要針對問題某些特性作挑選 11/30 17:13
124F：→ iHaveAPen: 還是要看依你的矩陣好壞來決定該用什麼方法，然後自己 11/30 19:03
125F：→ iHaveAPen: 寫lib有自己寫的好處，因為往往不是只要答案，還要做 11/30 19:03
126F：→ iHaveAPen: 其他事。所以也不用嘴為什麼我建議自己寫，我也沒有要 11/30 19:03
127F：→ iHaveAPen: 嘴直接使用lib的，依人所好 11/30 19:03
128F：→ saltlake: 用既有的函式還可能遇到該函式為了含括各種狀況而肥大 11/30 19:29
129F：推 sunev: 要看wohtp對程式有多不熟了，無論是自己寫還是用lib 11/30 19:29
130F：→ sunev: 都有許多「眉角」，有熟的人在旁可以問的話，可以省不少事 11/30 19:30
131F：→ saltlake: ，而且計算速度可能會低於針對某特定問題寫的函式 11/30 19:30
132F：→ saltlake: 另外，真要求嚴謹的話，不會"直接"呼叫別人的函式 11/30 19:32
133F：→ saltlake: 必須經過一定的測試程序之後才用 11/30 19:32
134F：→ HeterCompute: 自己寫的也需要測試，用別人的lib也需要測試，這應 11/30 22:43
135F：→ HeterCompute: 該是常識 11/30 22:44
136F：→ HeterCompute: 另外如果已經有能力自己從頭寫得好，也和別人的lib 11/30 22:46
137F：→ HeterCompute: 效能差不多，你當然可以做，但對於初學者建議重新造 11/30 22:48
138F：→ HeterCompute: 輪這想法我不同意，但我沒有要嘴人的意思，有能力 11/30 22:49
139F：→ HeterCompute: 自己造輪那當然自己造輪，好處我想大家都知道 11/30 22:50
140F：推 j0958322080: 看他急不急吧 11/30 23:25
141F：→ wohtp: 其實我現在傾向於相信我邊界條件下錯了 12/01 01:19
142F：→ wohtp: 然後relaxation才會壞掉。如果可以修好，我也不用搞牛頓法 12/01 01:20
143F：→ wohtp: 了。 12/01 01:21
144F：→ saltlake: 不管用誰的程式前都要測試沒錯，是否常識且不論，只針對 12/01 02:25
145F：→ saltlake: 用他人程式狀況提醒，乃因取此圖者一般更傾向省事 12/01 02:26
146F：→ Eriri: 在線性代數中 自己造輪子絕大多數都不會是好想法 這不是在 12/01 02:34
147F：→ Eriri: 嘴 這是對任何真的在做研究的人都需了解的 12/01 02:34
148F：→ Eriri: 這些傳統線性代數lib 都是由最專業的應用數學家跟電腦科學 12/01 02:37
149F：→ Eriri: 家 花費大量心力認真建構而成的 上至NASA上太空 下到華爾 12/01 02:37
150F：→ Eriri: 街跟矽谷 無數專業的科學家跟工程師 大量使用這些線性代數l 12/01 02:37
151F：→ Eriri: ib數十年以上 其優化跟穩定程度 一般人是比不上的 12/01 02:37
152F：→ Eriri: 當然 如果只是做個很簡單 大學生作業程度的事情 自己寫輪子 12/01 02:40
153F：→ Eriri: 那或許有額外好處 但是越複雜越龐大 就越該減少使用自己的 12/01 02:40
154F：→ Eriri: 輪子（以線性代數來說） 12/01 02:40
155F：→ Eriri: 例如 某個史丹佛教授是這麼說的 12/01 02:40
157F：→ Eriri: 越是了解數值分析 其實就會越了解這些線性代數lib的偉大 12/01 02:49
158F：→ Eriri: 我發現我昨天matlab的矩陣設錯了 看來這終究是能用執行時 12/01 10:44
159F：→ Eriri: 間來換取調整記憶體使用量 我的記憶體全部都被佔滿 但終究 12/01 10:44
160F：→ Eriri: 沒有out 只是我出去吃個飯回來都還沒跑完 要做事就先關了 12/01 10:44
161F：→ Eriri: 如果你時間多可以等 又只要找出解就好 之後不會做其他規模 12/01 10:52
162F：→ Eriri: 太大的事 （就像我以前玩個幾萬的矩陣就crash 是因為我還 12/01 10:52
163F：→ Eriri: 做其他大矩陣運算 甚至二維積分 而且我的不是sparse) 就 12/01 10:53
164F：→ Eriri: 像前面大家說的 或許還是有一般電腦就可以搞定的方法 你自 12/01 10:53
165F：→ Eriri: 己決定吧～ 12/01 10:53
166F：推 Eriri: 畢竟我也沒有真正specificly測過解線性方程在一般電腦的極 12/01 10:55
167F：→ Eriri: 限 我只知道一般N到了數萬以上的矩陣問題或操作 不管是用 12/01 10:55
168F：→ Eriri: 什麼語言或lib 在單台個人電腦很容易就等相當久 如果不是 12/01 10:55
169F：→ Eriri: 不能做的話 12/01 10:55
170F：噓 yzfr6: 演算法 記憶體 12/02 00:05
171F：→ dhtsai: 這篇想到之前搞翻intel的數值計算,發現intel的硬體bug 12/03 11:07
172F：→ dhtsai: Pentium FDIV bug,intel回收一整批處理器,花$475 million 12/03 11:10
173F：→ recorriendo: 既然友人都丟史丹佛了 我也丟一個 12/03 19:28
176F：→ recorriendo: 函式 不過沒有打包成漂亮的library 將就用 12/03 19:35
177F：→ wohtp: 回來回報一下，我好像可以用relaxation了 12/06 19:46
178F：→ wohtp: 所以不需要解牛頓法喔耶 12/06 19:46
179F：→ wohtp: 喵的沒人告訴我做relaxation的時候不能先把gauge定下來 12/06 19:48
180F：噓 ostracize: 02/03 00:34