作者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/cn.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