CCS HPC seminar 2024

筑波大学計算科学研究センター
CCS HPCサマーセミナー「並列数値アルゴリズム II」
多田野 寛人tadano@cs.tsukuba.ac.jp
筑波大学計算科学研究センター

CCS HPC seminar 2024

講義内容
・連立一次方程式の求解法- 様々な行列に対する解法- 疎行列の扱い方- 基本的な線形演算の MPI での並列化

・複数本の右辺ベクトルをもつ連立一次方程式の解法- 解法の特色- 効率的な計算手法と OpenMP での並列化

・レポート課題について-1-

CCS HPC seminar 2024

連立一次方程式の求解法-2-

CCS HPC seminar 2024

自然現象・工学現象の解析
自然現象・工学現象

偏微分方程式の近似解
解析
方程式の求解モデル化
偏微分方程式の初期値・境界値問題

離散化

連立一次方程式
Ax = b

連立一次方程式は様々な分野における数値シミュレーションで現れ，計算時間の大部分が求解に費やされている-3-

CCS HPC seminar 2024

連立一次方程式
連立一次方程式：Ax = b<latexit sha1_base64="GGQ+DS17QDj05GpsS5fwGPYM7GE=">AAANSHicnVZLj9tUFD4TKLQmpR1gUYnNhemMWEyD4w70IUVqVRYs+2DaSuNp5MdNejV+YTuZpFb+QP8AC1YgsUCIX8EGsWfRn4DYIBWpQmLBuedeT5ImTjzYin197ne+87w3100CkeWm+WKj8cabZ956++w5453m+XcvXNx872EWD1KP73txEKePXSfjgYj4fi7ygD9OUu6EbsAfuUd35PyjIU8zEUdf5eOEH4ZOPxI94Tk5irqbGz/f7hh2wHv5gWG7vC+iwklTZzwpPLwmBnO6RXu3PWE7amTJke3HecZKUTSxbcJZJzhrEWdpnD0k4c7MwJ+TaLLohCxaJIsUmc0jX7tr2KnoP80PjV1mM5xxw2I0YR1WERsGNuq2mbQ16lpsxjElithaehfpq9ldYmdrLolbsO5WWze6F7fMlkkXWxy09WAL9HU33jxzG2zwIQYPBhAChwhyHAfgQIb3AbTBhARlh1CgLMWRoHkOEzBQd4AojggHpUf47OPXgZZG+C05M9L20EqAvxQ1GWybv5s/mi/NX82fzD/Mfyu5CuKQvozx7SpdnnQvPL/04NVarRDfOTydaq30OYceXCdfBfqekERG4Sn94bNvXj64eX+72DG/N/9E/78zX5i/YATR8G/vh3v8/rfIvq1/1X4FOJJZEORZuCKPy9CCvKnWyGGEUcRUx4yQ2yuwQtdnFaOrfaxGTCubIzo4sVpmf4ioFDUF+SUZJd5faTXBdw+fkjdf8HB+LkBmFy3I/hyjbg+/A7xdxOZwjE9ZxVXWlIcu+SvW1EQijuBZLcwVuIpxlp7W0RDUbwGi/JVrrMBnhDMl9xXScfB9Gv8Cstdf0wGycgPECequWeQsz2wFEvQpI7yqubJSRzMkrE+9k9TS+BrnY1qpGeyiVFoKFrq6SjvD/pSWSu0U7xh7Rn55NCf9qcfl0Xr1amZI9lpOfVeXX9CqVmuJo50yu9LX6U5X109OeXZoP16v45PFmCRTuy2ym6LGsj6q4urhyD+VxmyVWtp/2feDhb1hNYOva9uifWpQM/ara/eqkd6TymycZrea7p6q58t4XtdTHip8qPuHQQd/Ec4cU27KPUZJFDKkNeGjpq2zLlHyS40Zjj6qeJboVZzD/8VZ12pZkzJ6uT/4eMv+na7WT7GiMvoxyTuIM5Bp/l7WwR1inO/g3SW65V1iZdZV/l/f6Riel1rwGb6P9X82o7xMfS1gC09VEzpTLLM1WSpVdnLqvFnv56NS/xhTJmNN7ULMWIY7aEr5PgCLTj82sQnKeZtOfzbGIf+1Rzh7GZ4QRp2u5DpkOqLLFFN1XUtLMgMWZUCeKtSp9QZdTA2u7enBjfbJqfWh1Wp/3tq7t7d1y9Ln17PwIXwMn6Dta3ALvoS7sA9e44PGzcadxhfN35p/NV81/1HQxobWeR/mrvON/wC3i83/</latexit>

2
666 a1,1
666 a
666 2,1A = 666 .
666 ..
4
an,1

a1,2
a2,2
..
.

...
...
..
.

an,2

...

3
2
3
2
3
a1,n 77
x
b
6
7
6
777
1
1
666
777
666
777
666 x2 777
666 b2 7777
a2,n 7776
7
6
7.. 7777 , x = 6666 .. 7777 , b = 6666 .. 7777
666 . 777666 . 777
. 7775
4
5
4
5
an,n
xn
bn

連立一次方程式は様々な分野における数値シミュレーションで現れ，計算時間の大部分が求解に費やされている-4-

CCS HPC seminar 2024

直接解法と反復解法
直接解法• Gaussの消去法，LU分解法など• 有限回の演算で必ず解くことができる• 係数行列 A を変形するため，非零要素数が増大

反復解法• 定常反復法，クリロフ部分空間反復法• 必要な演算は係数行列 A とベクトルの積，内積など
係数行列の疎性をそのまま使える• 問題によっては多くの反復回数を要することがある-5-

CCS HPC seminar 2024

代表的な直接解法（Gaussの消去法）
Gaussの消去法では，連立一次方程式 Ax = b を変形し，上三角行列 U をもつ
連立一次方程式にした後で解ベクトル x を計算する．
変形後の方程式： U x = b0

連立一次方程式： Ax = b<latexit sha1_base64="KgJAuzzRmuuJK8oaomGR2m6sQaY=">AAAMHXicnVZNb9xEGH63paXsQpOWSxEXQxrooV28m8C2SJGKEFUPIPUrbaVsG/ljdmvFX9jeTVLLJ278AQ4cEEgcED8DIXFA3BDqT6g4FokLB555x453k12vU49sj2ee5/2e8Zih68SJrj9rnDj5yqnTr555rdl6/Y2zS8vnzt+Pg1FkiU0rcIPooWnEwnV8sZk4iSsehpEwPNMVD8ydT+X8g7GIYifw7yX7oXjkGUPfGTiWkWBoe/mrvi0G4LKk1HWGTxLTHYksvXnvi8+zdK3Xu3GjlzWnULYR7QwjIfwD1NXeZ91u1kz7h8VsZn3TS/cybUMrJku2nDKz97Pt5RW9rfOlHe108s4K5det4NypT6hPNgVk0Yg8EuRTgr5LBsVoW9QhnUKMPaIUYxF6Ds8LyqgJ7ggoAYSB0R08h/jaykd9fEuZMbMtaHFxR2BqtKr/qf+kv9B/03/Wn+v/zZWVsgxpyz7epuKKcHvp6wt3/13I8vBO6EnJqrQ5oQFdZVsd2B7yiPTCUvzx029e3P34zmr6nv6D/jfs/15/pv8CD/zxP9aPt8WdbyF9Nb/n2+WiJ6PgsGVeRRxnoR22Zj4joT14EXAeY0auVmCdPD9VEs3cxvmIMrMJ0O6B1iL6Y6AiMB22S0qUeLtSa4j3AE8pNzli4fScC8kmNMj63Ad3gG8XzQQ2oV08ZRartCkLTbbXWZATidihp7UwV2gNfhaW1mE4XG8uUHblGkvx9DFTyL7CHAPv49jnsr7hggqQmRsB53B1TSIn5UxmIIRNMeNVzpWWOkyPsTbXTliL8SXmA16pMV3GqNTkHqnqeewY9Sk1FewILUDNyC+L56Q99WRZvF6tmhGStZZw3dWV7/CqVmtJQE8RXWlrudPVtVNwnA3ejxdzbNYY8Eipt816IzBm1dE8WQP07GMxJrPUzu2XdT86sjdUS7Dz3LZ5nxrV9H1t4V61l+9JRTSOs1uVu6eq+cKfwzxlocJ7ef1otIHbx8wux6bYY9SIQnq8Jmww+3nUJUp+qb6G3jtzngW6Sub4pWTW1VrkpPBe7g82mqzfcrV+gIxK7/d5fAO4JiRNt1kVvMESpyv48gxu0QqsjLqK/+GdTsN5qU0f4r2b/7M1jktpa0orOFVlfKaYpSubOar0JFx5k9ZPe6X+GKWk5oLceYhYjB004nhvUZdPP32W5nDMO3z668MP+dfew+xFeswYdbqS61DLPbrIPs3Pa6FJRqDLEZCniiafWq/xpalObz3vXOscnFrvd9udj9rrt9dXruv5+fUMvU3v0iXo7tF1ukm3aBM+Pm8sNS403mp91/q19XvrDwU90cg5b9LU1frrf9K3fZY=</latexit><latexit sha1_base64="+M8e1hJrA8uuOE2vXa/7eMCI5jU=">AAALqnicnVbLbtNAFL2UV0l4FNggsTGUIhZtmLSFUqRKIDYsoaVQ1EAVP1JGtWNjO2lLlB/gB1jABiSQEJ/Bhh9gwScglkViw4Izd+wmaRPHxSPb4zvnvu9cjxm4MoqF+HFo5PCRo8eOj54oFE+eOn1m7Oy5x5HfCC1n2fJdP1wxq5HjyrqzHMvYdVaC0Kl6pus8MTfuqfUnTSeMpF9/FG8HzjOvul6XNWlVY5DWxsbuVkyvtdU2Fgw1MdtrY+OiJPgy9k/KyWSckuuBf/boXaqQTT5Z1CCPHKpTjLlLVYowVqlMggLQnlELtBAzyesOtakA3gZQDhBVUDfwXMfXakKt41vJjJjbghYXdwhOgybEd/FZ7Ihv4ov4Kf4OlNViGcqWbbxNzesEa2deX1j6M5TLwzumFx2uTJtjqtEttlXC9oApygtL8zdfvdlZur040boqPohfsP+9+CG+woN687f18aGz+BbSJ5J7sF0uZioKki3zMuLYDy3ZmsEcMW3BC5/zGDFyIgMrk/xkSTQTGwcjOpmNgXZ3tabRbwIVglOyXUqiwtuZWgO8a3gqufE+C3vXXEg2oUHV5zZ4a/h2MUxgY9rEU2UxS5u20GR75ZCcKMQGvcqFmaIZ+JlamodDcr25QNmZe6yFZx0rqewp5qnifRD7XNa3PqQCVOYawEmurm5kt5zuDASwKWK8zrnWkofTY6zNtRPk4niJdZ93akSToCpN7r6qHsQdoT6VppQ7xPBRM+rL4jVlTz5ZFu9XK2eEVK3FXHd55Uve1XovOdCTRlfZ2ul0ee10OM5V7sfDeWzW6DOlo7fEekNw9KujQbJqmNkH4ujOUimxX9V9Y19vyJZgJ7ktcZ9q5PR9Zmiv2kp6UhqNg3SrTvfUNZ/6s5dPW6jxXlI/Bi3grmNlk2OT9hhN0UiP94QNzkoSdYVSX3puYHZpwDNFZ8ls/pfMvFrTnKTeq/5gY6j67ezW68io8n6b6QvAFSCpd/Sr4AWW2FvBk31405FiVdR1/Pd2OgPnpRLdwHsz+WcbHJeOrS0ax6mqzWeKfrrafalaT8yV1219r1f6j9GRVBiSOw8Ri9BBQ473Kk3z6afC0iTHvMynvwr8UH/tLaxeoeeM0acrtQ+NxKMr7NPgvKaaVASmOQLqVFHgU+s8X4aezM0mk/ny7qn18XSpfLM0+3B2/I5Izq+jdJEu0zXonqM7dJ8e0DJ8bNI7+kifipPFxeLT4qqGjhxKeM5Tz1W0/wHcaE9V</latexit>

2
666 a1,1
666 a
666 2,1666 ..
666 .
4an,1<latexit sha1_base64="3hL/2DAZY7PFz7Wj8dJUH9g3Q5g=">AAANKHicnVbNj9tEFH8bKLQm0C5cULkYsrvisA2Ou9AWaaUiLhz7wbaV1tvIH5N0tP7CnmSTRvkHKnHmwAkkDsCfwYUjF4T2yg1xbCUuHHjzZrxJmtjxYiv2+M3v/d7nTMZLQ54LyzrdaLzy6oXXXr94yXij+eZbl69svv0gTwaZzw78JEyyR56bs5DH7EBwEbJHacbcyAvZQ+/4czn/cMiynCfxl2KcsqPI7ce8x31XoKi7ufG1E7KeODQcj/V5PHGzzB1PJz5eU8N0u5PObmdq7qiRLUdOkIjcLETx1HEIZ5/h7GWcrXHOkIQ7c4NgQaLJ4jOyeJksVmQOiwPtruFkvP9EHBmrY8FARt2OKblHXducc0SJYrOEbr+c0CNCc80lcUsGvXKDRvdKy2pbdJnLg44etEBfd5LNC5+BAwEk4MMAImAQg8BxCC7keB9CByxIUXYEE5RlOOI0z2AKBuoOEMUQ4aL0GJ99/DrU0hi/JWdO2j5aCfGXoaYJ29bv1o/Wc+tX62frL+vfUq4JcUhfxvj2lC5Lu5efvXv/n7VaEb4FPJlpVfosoAc3yVeOvqckkVH4Sn/49Jvn9z+9tz3Zsb63/kb/v7NOrV8wgnj4wv/hLrv3LbJv61+5XyGOZBY4eRZV5HEVmpM35RoCRhhFQnXMCbldgeW6PlWMnvaxHDGrrEB0eGa1yP4QURlqcvJLMkp8UGk1xXcPn5JXLHm4OBcis4cWZH+OUbeH3yHeHmIFnOBTVrHKmvLQI3/5mppIxDE8rYW5BtcxzsLTOhqc+i1EVFC5xib4jHGm4L5GOi6+z+NfSPb6azpAVm6AOE7dNY+c55mvQIo+5YRXNVdW6mhGhA2od9JaGl/hfEIrNYddlEpL4VJXl2nn2J/SUqGd4Z1gz8gvn+akP/W4fFqvfs0MyV4T1Hd1+TmtarWWGNopsit9ne10df1klGeX9uP1OgFZTEgys9smuxlqrOqjMq4ejoJzacxXqa39l30/WNobqhkCXds27VODmrFfX7tXjfSeVGTjPLvVbPdUPV/E87Ke8lDhI90/JuzjL8aZE8pNsccoiUJGtCYC1HR01iVKfqmxiaP3S54Fuopz+L8461otalJEL/eHAG/Zv7PV+hFWVEY/Jvk+4gxkWrxXdfA+MS528O4K3eIusDLrKv8v73Qmnpfa8DG+T/R/tkl5mfk6gRaeqqZ0plhla7pSquwI6rx57xejUv8YMyZjTe0izFiOO2hG+T4Em04/DrFxynmHTn8OxiH/tUc4uwWPCaNOV3IdmjqiLYqpvK6FJZkBmzIgTxXq1HqLLlMNbuzpwa3O2an1gd3ufNLeu7vXum3r8+tFeA8+gA/R9g24DV/AHTgAf+NF42qj1dhq/tT8rflH81RBGxta5x1YuJp//gehscYS</latexit>

a1,2
a2,2
..
.

...
...
..
.

an,2

...

32
3 2
3
a1,n 77 66 x1 77 66 b1 77
76
7 6
7
a2,n 7777 6666 x2 7777 6666 b2 7777
76
7 6
7.. 7777 6666 .. 7777 = 6666 .. 7777
. 7775 6664 . 7775 6664 . 7775
an,n
xn
bn

2
666 u1,1
666
666
666
666
4

u1,2
u2,2

<latexit sha1_base64="DBiEsxIN0Jmrkews6dw9h7H0YvM=">AAAQU3icvVbLjttUGP47TKAEyszABokFhukUFtPgeKakRRqpCLXqgktv01YaD5EvJxlrfMM+ySS1/AK8AAtWILFAvAB7NrxAF30ExLJILGDBf/5jT5KJ7TiVqK3Yx+f/vv9+To4Zuk7MVfXpuZWXVhsvv3L+1eZrr194Y219480HcTCILLZvBW4QPTKNmLmOz/a5w132KIyY4Zkue2gefybkD4csip3Av8/HITv0jL7v9BzL4DjV3Vj9VbdZD8mkKnGd/hE33QFLk1v3v/g8TXY6nZs3O2lzBmUb0XE/Ysw/RV3t3NA0RPnspBf4PNHNfprwUaTEluEyW0erMWfhFYS4rMcPlBd7NXWT9R0/MaLIGKeJhVf6AqwqiX42rYNu0t5up8L6pTKxhuIioW4HPE4riH6a6nqz2qkSsva8VrWaVpd0Sbel2ULZUMqE1XlposeeER8l+pEZjEQjKmqapulCH5ZOmZ9FrjPfzjqrqUcCcZh1+dm2S5vKqNsmv0ddjd4ylmzKV0rU7VVn9/+5yoPI0zHZB8xu+4PZgkzLtFJZUSmnmb5kFuWku76ptlS6lPlBOxtsQnbdDjYan4IONgRgwQA8YOADx7ELBsR4H0AbVAhx7hASnItw5JCcQQpN5A4QxRBh4OwxPvv4dZDN+vgtdMbEttCKi78ImQpsqU/Un9Vn6u/qL+of6r+luhLSIXwZ49uUXBZ21759+97fC1kevjkcTViVPnPowVXy1UHfQ5oRUViSP3z83bN7n9zdSi6pP6p/ov8/qE/V3zACf/iX9dMddvd71L6V/cr9cnEksuCQZ15FHovQDnlTzuAwwigCqmNMyK0KrJPVp0qjmflYjphUliPaPbWaZ3+IqAiZDvklNAq8XWk1xHcPn0Ivn/NwVuaiZhMtiP4cI7eH3y7eJmI5nOBTVLHKmvTQJH+dBTURiGN4XAtzGXYwztzTOgyH+s1FlF25xhJ8+ijJdV8mjoHvZfxzyV5/QQeIyg0Q51B3TSOn9UxXIESfYsLLmksrdZgeYW3qnbAW4xuUB7RSY9jGWWHJnevqMnaM/Sks5ewI7wB7RnxZJBP+1NNl0Xq1amZI9Bqnvqur36FVLdcSQzt5doWvk52urp+M8mzQfryYY5PFgGYmdltkN0JGUR+V6erhyF6KMV2lVua/6PvB3N5QrcHOatuifWpQM/adhXvVKNuT8mwss1tNdk/Z83k8Z3nSQ4n3sv5RYA9/PkpOKDf5HiNnJNKjNWEjU8+yLlDiS44VHL1X8szRVTqHz6WzrtW8Jnn0Yn+w8Rb9O1mtH2FFRfRjmt9DXBM1zd5FHbxHGmc7eLuAm985VmRd5v/sTqfgeakFV/B9kv1nK5SXia8JbOKpKqUzRZGttHBW2uHUedPez0Yl/zEmmpoLaudhxmLcQSPK9wFodPrRSZtDOW/T6U/HOMS/9gilF+FrwsjTlViHShbRRYqpvK65JZEBjTIgThXy1HqNLkUOOrvZ4Fr79NT6QGu1P27t3tndvK5l59fz8A68Dx+i7Q5ch1twG/bBarzbuNH4svHV2pO1f9ZX1lcldOVcxnkLZq71C/8B062ZTw==</latexit>

変形

0

...
...
..
.

32
3
3 2
u1,n 77 66 x1 77 66 b01 77
76
7
7 6
u2,n 7777 6666 x2 7777 6666 b02 7777
76
7
7 6.. 7777 6666 .. 7777 = 6666 .. 7777
. 7775 6664 . 7775 6664 . 7775
un,n
b0n
xn

変形後は下記の後退代入計算により，解ベクトル x の成分を計算する．
⌘
1 ⇣ 0
b u j, j+1 x j+1 u j, j+2 x j+2 · · · u j,n xn , j = n, n 1, . . . , 1.
xj =u j, j j<latexit sha1_base64="BnRJ1AZuGqK0hsSGbGOqnnvdcog=">AAAMHXicnVZLj9tUFD4pFEoCzBQ2rdhcmA4UkQl2OrQUaaQiNl32wbSVxiWKH5neGb+wnUymUVbs+AMsWCCQWCB+BkLqArFDqD+hYlkkNl30u+fak6STOB58Zfv63O+8zz2+duzLNDOMx7VTL718+pVXz7xWb7z+xpsrq2ffupNG/cTxtp3Ij5J7djf1fBl625nMfO9enHjdwPa9u/b+F2r97sBLUhmFX2aHsXc/6O6GsiedbgZSZ/WbYWdPbAmrl3SdkTke9TujvebeeGz5Xi+7aHf2PhAbQhM/MsdiiJl6H9HaOa2taJYbZWmxFqqVUFiJ3H2QfdgUFoZSFTZFuGE2NbgpzJborK4ZLYMvcXxi5pM1yq8b0dnTn5NFLkXkUJ8C8iikDHOfupRi7JBJBsWg3acRaAlmktc9GlMdvH2gPCC6oO7juYuvnZwa4lvJTJnbgRYfdwJOQevGX8YvxlPjkfGr8cR4tlDWiGUoWw7xtjWvF3dWvj13+7+lXAHeGT2YcJXanFGPPmVbJWyPmaK8cDT/4OF3T29/dmt99L7xk/EP7P/ReGz8Bg/Cwb/Ozze9W99D+np+L7bLx0xFQbJlQUkc56ElW7OYI6MhvIg4jykj10uwMs9PmUQ7t3ExYpLZDGj/SGsR/QFQCTgl26UkKrxbqjXGu4enkpsds3B2zYdkGxpUfR6Ct4dvH8MGNqMDPFUWy7RpC222Vy7JiULs08NKmA26BD8LS6twSK43Hyi3dI+N8AyxUsjeYJ4u3iexz2d9u0sqQGWuD5zk6ppGTsuZzkAMm1LG65xrLVU4A8a6XDtxJY6vsR7xTk2pCarS5B+r6kXcKepTaSq4E4wINaO+HF5T9lST5fB+dSpGSNVaxnVXVb7kXa33kgc9RXSVrZNOV9VOj+Pc5X68nMdljRFTJnpbrDcBx7w6WiSrh5l7Io7pLLVy+1Xd94/1hnIJbp7bFvepfkXfLy3tVcO8JxXROEm3mnRPXfOFPy/yaQs1PsjrR9AW7hArBxybosdoikYGvCdccFp51BVKfem5wOzdBc8CXSZz8L9kVtVa5KTwXvUHF0PV72S3foyMKu8Pmb4FXB2SZse8Ct5iibMV3JzDW4wCq6Ku4/9ipxM4L7XoE7wP8n+24LhMbB3RGk5VYz5TzNM1nkvVejKuvGnrZ73Sf4yJpPqS3AWIWIoOmnC8d6jNpx+LpUmOucmnPwt+qL/2EKsX6CvG6NOV2oci9+gC+7Q4r4UmFYE2R0CdKup8ar3Kl9CTK5v55Kp5dGq9026Zl1ubNzfXrrXz8+sZeofeo4vQfYWu0XW6Qdvw8UltpXaudr7xQ+P3xh+NPzX0VC3neZtmrsbfzwEB5nR6</latexit>

上三角行列をもつ方程式への変形過程の計算量は O(n ) ，後退代入の計算量は O(n2 ) ．
<latexit sha1_base64="QuR3KIyk/BYIZrI2fjbtTHo05S0=">AAALpXicnVZJb9NAFH5lLQlLCxckLoZSFqkEpy2UIlUCceGCKEsLUlMqbymjeqvtpEvUP8Af4MAFkDggfgYX/gAHfgLiWCQuHPjmjU0SmjguHtkev/ne/uZ5zNAVcaLr34YOHDx0+MjR4WOl8vETJ0+NjJ5ejINGZDkLVuAG0XPTiB1X+M5CIhLXeR5GjuGZrvPMXLsn1581nSgWgf802QqdZc9Y9UVdWEYC0nKrZhmu9nDniv9i6urKyJhe0fnS9k6q6WSM0ms+GD18l2pkU0AWNcgjh3xKMHfJoBhjiaqkUwjaMrVAizATvO7QDpXA2wDKAcIAdQ3PVXwtpVQf31JmzNwWtLi4I3BqNK5/1T/qu/oX/ZP+Xf/dV1aLZUhbtvA2Fa8Trpx6dfbJr4FcHt4JvWxz5dqcUJ1usa0CtodMkV5Yir+5/Xr3ye3H461L+nv9B+x/p3/TP8MDv/nT+vDIefwG0sfTu79dLmYyCoIt83Li2Ast2Jr+HAltwouA8xgzcjwHK9L85Ek0Uxv7I9qZTYB2/2rNot8EKgKnYLukRIm3c7WGeNfxlHKTPRZ2r7mQbEKDrM8t8Nbx7WKYwCa0gafMYp42ZaHJ9ooBOZGINdouhLlGU/Azs7QIh+B6c4Gyc/dYC08fK5nsa8xj4L0f+1zWtzqgAmTmGsAJrq5OZKeczgyEsClmvMq50lKE02OszbUTFuJYx3rAOzWmCVClJndPVffjjlGfUlPGHWEEqBn5ZfGatKeYLIv3q1UwQrLWEq67ovIF72q1lxzoyaIrbW13uqJ2Ohxng/vxYB6bNQZMaeutsN4IHL3qqJ+sOmb2vjg6s1RJ7Zd139jTG/Il2GluK9ynGgV9nxrYqzbTnpRFYz/dqt09Vc1n/vzLpyxUeC+tH43mcPtY2eDYZD1GURTS4z1hg7OWRl2i5Jeaa5id7/PM0Hkym/8ls6jWLCeZ97I/2Biyftu79ToyKr3fYvoccCVI6h69KniOJXZX8EQP3mxkWBl1Ff9/O52G81KFbuC9kf6zNY5L29YWjeFUtcNnil66dnpSlZ6EK6/T+m6v1B+jLak0IHceIhajg0Yc7yWa5NNPjaUJjnmVT381+CH/2ptYvUgvGKNOV3IfaqlHF9mn/nnNNMkITHIE5KmixKfWWb40NZmZTiez1b+n1sXJSvVmZfrR9NidyfT8Okzn6AJdge4ZukP3aZ4W4OM6vaa39K58ufyg/LS8qKAHhlKeM9R1lVf+AJ8pTW0=</latexit>

3<latexit sha1_base64="VF9fHwoLJUSr+fNjUDIGHpB8ROc=">AAALpXicnVbLbtNAFL08CiXhVdggsTGU8pAgOGl5SpVAbNggSksLUlMqv1JG9au2k7ZE/QF+gAUbisQC8Rls+AEWfAJiWSQ2LDhzxyYJTRwXj2yP75z7vnM9ZuiKONH1b/v2Hzg4cujw6JFS+eix4ydOjp1aiINmZDnzVuAG0XPTiB1X+M58IhLXeR5GjuGZrvPMXH0g15+1nCgWgf802QydJc9Y8UVDWEYC0lK7bhmu9njrsv+idmX55Lhe0fnSdk+q6WSc0msmGBu5T3WyKSCLmuSRQz4lmLtkUIyxSFXSKQRtidqgRZgJXndoi0rgbQLlAGGAuornCr4WU6qPbykzZm4LWlzcETg1mtC/6h/1Hf2L/kn/rv8eKKvNMqQtm3ibitcJl0+8PjP3ayiXh3dCLztcuTYn1KDbbKuA7SFTpBeW4m+9erMzd3d2on1Rf6//gP3b+jf9MzzwWz+tD0+c2beQPpHeg+1yMZNREGyZlxPHfmjB1gzmSGgDXgScx5iREzlYkeYnT6KZ2jgY0clsArT7V2sW/RZQETgF2yUlSrydqzXEu4GnlJvssrB3zYVkExpkfW6Ct4FvF8MENqF1PGUW87QpC022VwzJiUSs0qtCmGs0CT8zS4twCK43Fyg7d4+18fSxksm+xjwG3nuxz2V9K0MqQGauCZzg6upGdsvpzkAIm2LGq5wrLUU4PcbaXDthIY41rAe8U2O6CqrU5O6q6kHcMepTasq4I4wANSO/LF6T9hSTZfF+tQpGSNZawnVXVL7gXa32kgM9WXSlrZ1OV9ROh+NscD8ezmOzxoApHb0V1huBo18dDZLVwMzeE0d3liqp/bLum7t6Q74EO81thftUs6Dvk0N71Ubak7Jo7KVbdbqnqvnMn3/5lIUK76X1o9E0bh8r6xybrMcoikJ6vCdscNbTqEuU/FJzDbNzA54ZOk9m679kFtWa5STzXvYHG0PWb2e3XkdGpfebTJ8GrgRJvaNfBU+zxN4KvtqHNxsZVkZdxf/fTqfhvFShG3ivp/9sjePSsbVN4zhVbfGZop+urb5UpSfhyuu2vtcr9cfoSCoNyZ2HiMXooBHHe5FqfPqpszTBMa/y6a8OP+RfewOrF+gFY9TpSu5DLfXoAvs0OK+ZJhmBGkdAnipKfGq9w5emJrem0smd6t9T60KtUr1ZmXoyNX6vlp5fR+ksnafL0H2L7tFDmqF5+LhGb+gdbZcvlR+Vn5YXFHT/vpTnNPVc5eU/k/NNbA==</latexit>-6-

CCS HPC seminar 2024

代表的な直接解法（LU分解法）
LU分解法は係数行列 A を下三角行列 L と上三角行列 U の積に分解することで，Ax = b を解く方法である．係数行列だけが変形され，右辺ベクトルは変わらない．
変形後の方程式： LU x = b<latexit sha1_base64="m1hD+SqDR7T+bW/qW1vTy3L77XI=">AAAMGXicnVZNb9xEGH63UCi7QFM+JCQuhjSIQ7vMJoGlSJGKgqoeQOpX2krZEq3t2XQUf2F7N0ktizt/gENPIHFA/ArEpeq14tAfwAFxLBIXDjzzjt3dTXa9Dh7ZHs88z/s947EjTyWpEE8bp1548fRLL595pdl69bXXzy6de+N2Eg5jR245oRfGd+1+Ij0VyK1UpZ68G8Wy79uevGPvber5OyMZJyoMbqWHkbzn93cDNVBOP8XQztK3PVcOwGVJmRMrPwmDPLt666sv8+yLzc762mbenMJ4avd+antDWaLWut0rV7p5M+tNC8FU7yhlK+/ZfnaQWxuW7tj5ztKyaAu+rOOdTtFZpuK6Fp47/Tn1yKWQHBqST5ICStH3qE8J2jZ1SFCEsXuUYSxGT/G8pJya4A6BkkD0MbqH5y6+tovRAN9aZsJsB1o83DGYFq2I38XP4pl4JH4Rf4p/58rKWIa25RBv23BltHP2u3du/rOQ5eOd0v0xq9LmlAb0KduqYHvEI9oLx/BHD75/dvOzGyvZB+JH8Rfs/0E8Fb/Bg2D0t/PTdXnjIaSvFPd8uzz0dBQUW+ZXxHEWWrE18xkpHcCLkPOYMHKlAquK/FRJtAsb5yPGmU2B9p5rLaM/AioGU7FdWqLGu5VaI7wHeGq56TELp+c8SLahQdfnIbgDfHtoNrAp7eOps1ilzVhos71qQU40Yo8e1MJcpDX4WVpah6G43jyg3Mo1luEZYKaUfZE5fbxPYp/H+nYXVIDO3BA4xdU1iZyUM5mBCDYljDc5N1rqMH3Gulw7US3GN5gPeaUmdAGjWpN3rKrnsRPUp9ZUsmO0EDWjvxye0/bUk+XwenVqRkjXWsp1V1e+4lVt1pKEnjK62tbxTlfXTslx7vN+vJjjssaQR8Z626w3BmNWHc2TNUDPPRFjMkvtwn5d98Nje0O1BLfIbZv3qWFN39cW7lUHxZ5URuMku9V49zQ1X/pzlGcsNHi/qB+LNnAHmNnn2JR7jBkxSJ/XhAtmr4i6Rukv07fQe2/Os0RXyRz9L5l1tZY5Kb3X+4OLput3vFo/Qka194c8vgFcE5Km26wK3mCJ0xV8YQa3bCVWR93E/+hOZ+G81KaP8d4v/tkWx2Vsa0bLOFXlfKaYpSufOWr0pFx5k9ZPe2X+GGNJzQW58xGxBDtozPHeplU+/fRYmuKYd/j014Mf+q99gNnz9DVjzOlKr0Or8Og8+zQ/r6UmHYFVjoA+VTT51HqJL8t0uutF51Ln+an19mq780l7/fr68mVRnF/P0Lv0Pn0I3V26TFfpGm3Bxz8arcZbjbdbD1u/th61HhvoqUbBeZOmrtaT/wDoLHw6</latexit>

2
666 1
666 l
666 2,1666 ..
666 .
4
ln1

0

<latexit sha1_base64="lrPCwFbGpmRqexf4svWMj4j2riw=">AAAQgHicnVZLj9xEEK4sGQjDI7PkgsTFsNlVDpuJ7d3NkEgrBS2KcghSXptEWi8jP3pmrfVjsD2vGP8B/gAHTiBxQPwMLvwBDvkBHBDHIHHhQHW1Hc/D9njSlu12V9VXX1dVt9sYOHYYyfLLCxtvXWy8/c6ld5vvvf/Bh5dbmx89Df1hYLJj03f84Lmhh8yxPXYc2ZHDng8CpruGw54Z50dc/mzEgtD2vSfRdMBOXb3v2T3b1CMc6m5e/FOzWA+NCSo2A9sNfS+J7z356n4Sf3mk7O8dJc05Hcfun0WGM2SZ1l6nc/duB7U8Nu75XhRrRj+Jo0kghabuMEtDn2HEBgeo4rBedNLUDNa3vVgPAn2axCa2pCnF2jwHJZFE25FKW7loZwkv1kJXD89izfHHLFDb6gGbaGeGP+F8JTnBJmnaMg+nG6u7CgqXIWtQ3CnE1EaWH4VJEc0qkSVExZjI05OKeXLRrloo0kq9KZgu5llplppawPN+WiuFeYkMu7FCwZv1sCCeJ5YLl6ktGHpJwgNR2UqM1Tf1qtb0uialPLUFsrQiZpM+swizsl4o5VUc1g6ZJ4mZ1y8KrIhJVyHek65KbzGXdMiTSuAOywGNFNBYBjTKAbutLbktU5OWO0ra2YK0PfA3G1+ABhb4YMIQXGDgQYR9B3QI8ToBBWQY4NgpxDgWYM8mOYMEmmg7RC2GGjqOnuOzj18n6aiH3xwzJGsTvTh4B2gpwbb8h/yL/Er+Xf5V/kv+rxQrJgzOZYpvQ9iyQffydx8//nellYvvCM5yq0rOEfTgc+JqI/cBjfBZmMJ+9OL7V49vP9qOd+Sf5L+R/4/yS/k3nIE3+sf8+SF79AOib6d3OS8HezwKNjFzK+JYpG0Tm3KLCCY4C5/yGJLmdoWuneanCtFIOZZr5JmNUNt57TWL/gi1ArS0iRdH5PpWpdcBvnv45LjREsN5mYPIBnrg9TlF2x5+O3gZqBvBGJ88i1XeBEOD+NorcsI1zuFFLZ3rsIfzzJjWsbCp3hzUsirXWIxPDyUZ9nWy0fG9Dj+H/PVXVADP3BD1bKquWc1ZnNkMDJBTSPoi58JLHUuXdC2qnUEti29Q7tNKDWEXR7knZ6mqy6xDrE/uKbMO8PKxZviXSTLOpx6WSevVrBkhXmsR1V1dfJtWtVhLDP1k0eVc852uLk9GcdZpP15tY5FHn0Zyv23yG6BFUR2VYfWwZ61lMZuldsqf1/1waW+oRrDS3LZpnxrWnPveyr1qku5JWTTW2a3y3VPUfDafRTvBUOi7af1IcIi3h5IxxSbbY8SI0HRpTVhoqaVR51r8S/Ql7H1a8sy0qzBHb4RZ12uWk2z2fH+w8OL1m6/WG5hRPvspjR+iXhOR5q+iCj4kxPkK3i2wza5Ml0ddxH9xp5PwvNSGA3yP03+2RHHJucawhaeqhM4URb6SwlHhJ6LKm2U/Pyvxx8iRmity52LEQtxBA4r3Cah0+tEIzaaYK3T603Ae/K89QelV+Jp0xOmKr0MpndFVmlN5XjNPPAIqRYCfKsSp9RY1SXQ6+2nnlvL61PpUbSs32/sP97fuqOn59RJ8Ap/BNfTdgTtwDx7AMZiN+42gETe+bW20rrVutBShunEhtbkCc611+38VLOND</latexit>

1
..
.
ln,2

..

.
...

1

解ベクトルの計算手順

32
777 666 u1,1
777 666
777 666
777 666
777 666
54

0

u1,2
u2,2

...
...
..
.

32
3 2
3
u1,n 77 66 x1 77 66 b1 77
76
7 6
7
u2,n 7777 6666 x2 7777 6666 b2 7777
76
7 6
7.. 7777 6666 .. 7777 = 6666 .. 7777
. 7775 6664 . 7775 6664 . 7775
unn
xn
bn

1. y = U x とおいて Ly = b を前進代入で解き，y を求める．
2. U x = y を後退代入で解く．
<latexit sha1_base64="uh72FqwhLq3RCOyQ5W+ygpAr7I8=">AAAMBXicnVY9j9RGGH6XhEB2QzhIEymNk+MQBWy8d5dsiHQS6CJEQSS+DpBu4bS2Z4/R+Su2d+8Wy3Uk/kCKVIlEgSI6REuRJm2KFPyEKCVINBQ8847N7nK7Xl88sj2eeZ73e8Zjha6ME9N8UTv0wYeHPzpy9ON645Njnx5fOHHyVhz0I1ts2IEbRHesbixc6YuNRCauuBNGoutZrrht7ayr+dsDEcUy8G8mw1Dc9brbvuxJu5tgaGthp+OIHrgsKbUj6cWBn6WXb/54JUt/WG+trqxn9QmMK7fvJ5bbFwVqpd2+dKmd1dPOpJArWcfy0mFmrBmqY2VbC4tm0+TL2N9p5Z1Fyq+rwYnDF6lDDgVkU588EuRTgr5LXYrRNqlFJoUYu0spxiL0JM8LyqgObh8oAUQXozt4buNrMx/18a1kxsy2ocXFHYFp0JL5j/nYfGn+Zf5h/mu+mSkrZRnKliHeluaKcOv4w89vvJ7L8vBO6P6IVWpzQj36jm2VsD3kEeWFrfmDB7+8vPH99aX0tPm7+R/s/818Yf4JD/zBK/vRNXH9V0hfyu/ZdrnoqShItswrieM0tGRrZjMS2oMXAecxZuRSCVbm+SmTaOU2zkaMMpsA7b7TWkR/AFQEpmS7lESFd0q1hnj38FRyk30WTs65kGxBg6rPIbg9fLtoFrAJ7eKpslimTVtosb1yTk4UYoceVMKcoxX4WVhahSG53lygnNI1luLpY6aQfY45XbwPYp/L+rbnVIDKXB84ydU1jhyXM56BEDbFjNc511qqMD3GOlw7YSXGT5gPeKXGdBajSpO7r6pnsWPUp9JUsCO0ADWjvmyeU/ZUk2XzerUrRkjVWsJ1V1W+5FWt15KAniK6ytbRTlfVTsFx7vJ+PJ/jsMaAR0Z6m6w3AmNaHc2S1UPPORBjPEvN3H5V9/19e0O5BCfPbZP3qX5F31fm7lV7+Z5UROMgu9Vo99Q1X/jzPk9bqPFeXj8GreH2MbPLsSn2GD2ikR6vCQfMTh51hVJfum+g9+WMZ4Eukzn4XzKrai1yUniv9gcHTdXvaLV+jYwq74c8vgZcHZIm27QKXmOJkxV8dgq3aAVWRV3H//2dzsB5qUnf4L2b/7MNjsvI1pQWcarK+EwxTVc2dVTrSbjyxq2f9Er/MUaS6nNy5yFiMXbQiOO9Sct8+umwNMkxb/HprwM/1F97D7On6B5j9OlKrUMj9+gU+zQ7r4UmFYFljoA6VdT51HqeL0N32qt553zr3an11nKz9W1z9drq4oXl/Px6lL6gr+gMdLfpAl2mq7QBH/+m1zWq1Ro/N540njaeaeihWs75jCauxvO3jipz9A==</latexit><latexit sha1_base64="L/XOHRVq1xb0JT6GP5lqisorKu8=">AAAMB3icnVZNb9RGGH4TSkt3oYRyqdSL2xDUAyzeJO2WSpFAqRCHVuIrgJSFaG3PhlH8he3dZLH8A+APcOgJpB6qHnqreq166b3iwE+oeqRShcShz7xjs7vJrtfBI9vjmed5v2c8VujKODHNl3PzR947+v4Hxz6s1Y+f+OjkwqmPb8dBL7LFhh24QXTX6sTClb7YSGTiirthJDqe5Yo71s66mr/TF1EsA/9WMgjFPa+z7cuutDsJhrYWvLYjuuCypNSOpBcHfpZevfX9d1n67XpzdWU9q41hXLn9ILHcnihQK63WlSstoCwvHWTGmpG29yM3MjW5l20tLJoNky/jYKeZdxYpv64Fp45epjY5FJBNPfJIkE8J+i51KEbbpCaZFGLsHqUYi9CTPC8ooxq4PaAEEB2M7uC5ja/NfNTHt5IZM9uGFhd3BKZBS+YL8yfzlfmn+bP5t/lmqqyUZShbBnhbmivCrZNPPrn530yWh3dCD4asUpsT6tLXbKuE7SGPKC9sze8/evrq5jc3ltKz5nPzH9j/zHxp/gEP/P6/9o/XxY0fIH0pv6fb5aKnoiDZMq8kjpPQkq2ZzkhoD14EnMeYkUslWJnnp0yilds4HTHMbAK0+1ZrEf0+UBGYku1SEhXeKdUa4t3FU8lNDlg4PudCsgUNqj4H4Hbx7aJZwCa0i6fKYpk2baHF9soZOVGIHXpUCXOeVuBnYWkVhuR6c4FyStdYiqePmUL2eeZ08D6MfS7r255RASpzPeAkV9coclTOaAZC2BQzXudca6nC9BjrcO2ElRgPMR/wSo3pHEaVJvdAVU9jx6hPpalgR2gBakZ92Tyn7Kkmy+b1aleMkKq1hOuuqnzJq1qvJQE9RXSVrcOdrqqdguPc4f14NsdhjQGPDPU2WG8ExqQ6miari55zKMZolhq5/aruewf2hnIJTp7bBu9TvYq+r8zcq/byPamIxmF2q+HuqWu+8Gc/T1uo8V5ePwat4fYxs8uxKfYYPaKRHq8JB8x2HnWFUl+6b6D32ZRngS6T2X8nmVW1FjkpvFf7g4Om6ne4Wi8go8r7AY+vAVeDpPE2qYLXWOJ4BZ+bwC1agVVR1/Hfv9MZOC816Eu8d/N/tsFxGdqa0iJOVRmfKSbpyiaOaj0JV96o9eNe6T/GUFJtRu48RCzGDhpxvDdpmU8/bZYmOeZNPv214Yf6a+9h9gzdZ4w+Xal1aOQenWGfpue10KQisMwRUKeKGp9aL/Jl6E5rNe9cbL49td5ebjS/aqxeX128tJyfX4/Rp/Q5fQHdLbpEV+kabcDHv+j13Pzckfrj+i/1X+u/aej8XM45TWNX/ff/AS1JdOw=</latexit><latexit sha1_base64="cUSiLBOdy1wLuFqo6myKgVj5pBk=">AAAMB3icnVZNb9RGGH4ToKW7UEJ7qcTFJQRxgK03SbulUiRQKsSBSnwFkLIQre3ZMIq/ans3WSz/gPYPcOiplXpAPfRW9Yp66b3iwE+oeqRShcSBZ96x2d1k1+vUI9vjmed5v2c8VujKODHNl3PzR44ee+/94x/U6idOfnhq4fRH9+KgF9liww7cIHpgdWLhSl9sJDJxxYMwEh3PcsV9a2ddzd/viyiWgX83GYTiodfZ9mVX2p0EQ1sLXtsRXXBZUmpH0osDP0uv3/3mRpZ+vd5cXVnPamMYV24/Tiy3JwrUSqt17Vorq6Xt/YCNrG156V5mrBmqM8i2FhbNhsmXcbDTzDuLlF83g9PHrlKbHArIph55JMinBH2XOhSjbVKTTAox9pBSjEXoSZ4XlFEN3B5QAogORnfw3MbXZj7q41vJjJltQ4uLOwLToCXzhfnMfGX+af5i/m2+mSorZRnKlgHeluaKcOvU95/c+W8my8M7ocdDVqnNCXXpS7ZVwvaQR5QXtub3nzx9deer20vpefMn8x/Y/6P50vwDHvj9f+2fb4nbP0D6Un5Pt8tFT0VBsmVeSRwnoSVbM52R0B68CDiPMSOXSrAyz0+ZRCu3cTpimNkEaPed1iL6faAiMCXbpSQqvFOqNcS7i6eSmxywcHzOhWQLGlR9DsDt4ttFs4BNaBdPlcUybdpCi+2VM3KiEDv0pBLmEq3Az8LSKgzJ9eYC5ZSusRRPHzOF7EvM6eB9GPtc1rc9owJU5nrASa6uUeSonNEMhLApZrzOudZShekx1uHaCSsxvsV8wCs1posYVZrcA1U9jR2jPpWmgh2hBagZ9WXznLKnmiyb16tdMUKq1hKuu6ryJa9qvZYE9BTRVbYOd7qqdgqOc4f349kchzUGPDLU22C9ERiT6miarC56zqEYo1lq5Paruu8d2BvKJTh5bhu8T/Uq+r4yc6/ay/ekIhqH2a2Gu6eu+cKf/TxtocZ7ef0YtIbbx8wux6bYY/SIRnq8Jhww23nUFUp96b6B3qdTngW6TGb/f8msqrXISeG92h8cNFW/w9X6GTKqvB/w+BpwNUgab5MqeI0ljlfwxQncohVYFXUd//07nYHzUoM+x3s3/2cbHJehrSkt4lSV8Zlikq5s4qjWk3DljVo/7pX+Ywwl1WbkzkPEYuygEcd7k5b59NNmaZJj3uTTXxt+qL/2HmbP0SPG6NOVWodG7tE59ml6XgtNKgLLHAF1qqjxqfUyX4butFbzzuXmu1PrveVG84vG6q3VxSvL+fn1OJ2hs3QBult0ha7TTdqAj3/R67n5uSP17+q/1n+r/66h83M552Mau+rP3wI7zHTs</latexit><latexit sha1_base64="kI2r+Q3PVGXdmuxwlS9Z0ebbmOg=">AAALunicnVY7b9NQFD68IQFaYEFiMZQikNrgpOUpKhWxsFEeBaQGKr9SruoXtpO2RFkZ+AMMTCAxIFjhD7DwBxj4CYgRJBYGvnuuTZI2cVx8Zfv63O+8zz2+ZuiKONH1bzt27tq9Z+++/QdK5YOHDo+NHzl6Pw6akeUsWoEbRA9NI3Zc4TuLiUhc52EYOYZnus4Dc/WGXH/QcqJYBP69ZCN0HnnGii8awjISkJbHtXbdMlztVqd91n88c64z1Ueonessj0/oFZ0vbeukmk4mKL0WgiN7rlOdbArIoiZ55JBPCeYuGRRjLFGVdApBe0Rt0CLMBK871KESeJtAOUAYoK7iuYKvpZTq41vKjJnbghYXdwROjSb1r/o7/af+RX+vf9f/DJXVZhnSlg28TcXrhMtjL47f/T2Sy8M7oSddrlybE2rQZbZVwPaQKdILS/G3nr38effqncn2Gf2N/gP2v9a/6Z/hgd/6Zb297dx5BemT6T3cLhczGQXBlnk5cRyEFmzNcI6E1uFFwHmMGTmZgxVpfvIkmqmNwxHdzCZAu/+0ZtFvARWBU7BdUqLE27laQ7wbeEq5yRYL+9dcSDahQdbnBngb+HYxTGATWsNTZjFPm7LQZHvFiJxIxCo9K4SZphn4mVlahENwvblA2bl7rI2nj5VM9jTzGHhvxz6X9a2MqACZuSZwgqurF9krpzcDIWyKGa9yrrQU4fQYa3PthIU4nmI94J0a0xSoUpO7paqHcceoT6kp444wAtSM/LJ4TdpTTJbF+9UqGCFZawnXXVH5gne12ksO9GTRlbZ2O11ROx2Os8H9eDSPzRoDpnT1VlhvBI5BdTRMVgMze1scvVmqpPbLum9u6Q35Euw0txXuU82Cvs+M7FXraU/KorGdbtXtnqrmM3828ykLFd5L60ejOdw+VtY4NlmPURSF9HhP2OCsp1GXKPml5hpmJ4c8M3SezNZ/ySyqNctJ5r3sDzaGrN/ubj2PjErvN5g+B1wJkvrHoAqeY4n9FTw1gDcbGVZGXcV/c6fTcF6q0AW819J/tsZx6drapgmcqjp8phikqzOQqvQkXHm91vd7pf4YXUmlEbnzELEYHTTieC9RjU8/dZYmOOZVPv3V4Yf8a69j9TQ9Zow6Xcl9qKUenWafhuc10yQjUOMIyFNFiU+tV/jS1OTSbDq5Uv13ar1fq1QvVmZvz07M19Lz6346QafoLHRfonm6SQu0CB+f0wf6SJ/K18pmWZRXFXTnjpTnGPVd5eQvElxU2g==</latexit>

32

行列 A のLU分解と解ベクトルの計算量はそれぞれ，O(n ), O(n ) である．
右辺ベクトルが複数ある場合は，LU分解法を用いた方が効率良く求解できる．

-7-

CCS HPC seminar 2024

疎行列の直接解法ソフトウェア– MUMPS（MPI）• http://mumps.enseeiht.fr/– SuperLU（SuperLU_DIST で MPI並列も可）• http://crd-legacy.lbl.gov/~xiaoye/SuperLU/– Pardiso（MPI, スレッド並列も可）• Intel Math Kernel Library (MKL) にも入っている• http://www.pardiso-project.org/– WSMP（MPI, スレッド並列も可）• http://researcher.watson.ibm.com/researcher/view_project.php?id=1426-8-

CCS HPC seminar 2024

クリロフ部分空間とは？
クリロフ部分空間
行列 A と非ゼロベクトル v から生成されるベクトル列で張られる部分空間：!

k−1

Kk (A, v) ≡ span v, Av, . . . , A"v

をクリロフ部分空間という．

クリロフ部分空間反復法

アレクセイ・クリロフ（Wikipediaより）• クリロフ部分空間を用いて連立一次方程式 Ax = b の近似解を生成する方法をクリロフ部分空間反復法という．
• 演算の主要部は行列 A とベクトルの積で，計算量は αn2 .
（ α は係数行列の非零要素数の割合）-9-CCS HPC seminar 2024

クリロフ部分空間反復法• 任意の初期解 x0 を与え，反復によって近似解を更新し，真の解 x を探索• r0 は初期残差と呼ばれ，r0 = b – Ax0 で計算される• 一般に真の解 x は分からないため，残差 rk = b – Axk のノルムを用いて真の解への収束をチェックする
x0 + Kk (A; r0 )x0 + Kk−1 (A; r0 )x0 + K3 (A; r0 )x0 + K2 (A; r0 )x0 + K1 (A; r0 )
•x1
x0•

•

x2

...

• x3

クリロフ部分空間反復法の概念図

xk

•
•
x

x

•

k−1

- 10 -CCS HPC seminar 2024

エルミート行列に対する解法
1. 係数行列がエルミート行列（A = AH）の場合
・共役勾配法（Conjugate Gradient method: CG法）
・共役残差法（Conjugate Residual method: CR法）
・最小残差法（Minimal Residual Method: MINRES法）

係数行列のエルミート性を使うことで短い漸化式（計算量が少ない）のアルゴリズムが導出できる

補足：エルミート行列

A = AH = ĀT(ai j = ā ji )- 11 -CCS HPC seminar 2024

共役勾配法（CG法）のアルゴリズム
x0 is an initial guess,
Compute r0 = b − Ax0 ,Set p0 = r0 ,For k = 0, 1, . . . , until ∥rk ∥2 ≤ εTOL ∥b∥2 do :qk = Apk ,行列ベクトル積(rk , rk )
αk =
,
内積
( pk , qk )xk+1 = xk + αk pk ,ベクトルの定数倍と加算
rk+1 = rk − αk qk ,(rk+1 , rk+1 )
βk =
,
内積
(rk , rk )pk+1 = rk+1 + βk pk ,ベクトルの定数倍と加算
End For- 12 -CCS HPC seminar 2024

相対残差, ∥rk ∥2 /∥b∥2

CG法の相対残差履歴のイメージ

このイメージでは∥rk ∥2 /∥b∥2 ≤ 10−12を満たしたら停止

反復回数, k・反復の過程で，相対残差（ ∥rk ∥2 /∥b∥2 ）をチェックする．
・自分で定めた条件を満たしたら反復をやめ，xk を解として採用．

- 13 -CCS HPC seminar 2024

非エルミート行列に対する解法
2. 係数行列が 非エルミート行列（ A ! AH ）の場合
短い漸化式を用いる解法
・BiCG法・CGS法・BiCGSTAB法
・GPBiCG法 ・BiCGSTAB( ℓ )法 ・IDR(s)法など

計算量は少ないが，残差ノルムは単調減少しない
長い漸化式を用いる解法
・GCR法

・GMRES法

残差ノルムは単調減少するが，計算量が多い- 14 -CCS HPC seminar 2024

BiCG法のアルゴリズム
x0 is an initial guess,
Compute r0 = b − Ax0 ,Choose r∗0 such that (r∗0 , r0 ) ! 0,
Set p0 = r0 and p∗0 = r∗0 ,For k = 0, 1, . . . , until ∥rk ∥2 ≤ εTOL ∥b∥2 do:
qk = Apk ,q∗k = AH p∗k ,(r∗k , rk ),αk = ∗( pk , qk )xk+1 = xk + αk pk ,rk+1 = rk − αk qk , r∗k+1 = r∗k − ᾱk q∗k ,(r∗k+1 , rk+1 )βk =,(r∗k , rk )pk+1 = rk+1 + βk pk , p∗k+1 = r∗k+1 + β̄k p∗k ,End For

行列ベクトル積
内積
ベクトルの定数倍と加算- 15 -CCS HPC seminar 2024

BiCGSTAB法のアルゴリズム
x0 is an initial guess,
Compute r0 = b − Ax0 ,Choose r∗0 such that (r∗0 , r0 ) ! 0,
Set p0 = r0 ,For k = 0, 1, . . . , until ∥rk ∥2 ≤ ∥b∥2 do:
qk = Apk ,(r∗0 , rk )αk = ∗,(r0 , qk )t k = rk − αk qk ,sk = At k ,(sk , t k ),ζk =(sk , sk )xk+1 = xk + αk pk + ζk t k ,rk+1 = t k − ζk sk ,αk (r∗0 , rk+1 )βk =·,ζk (r∗0 , rk )pk+1 = rk+1 + βk (pk − ζk qk ),End For

行列ベクトル積
内積
ベクトルの定数倍と加算- 16 -CCS HPC seminar 2024

GCR法のアルゴリズム
x0 is an initial guess,
Compute r0 = b − Ax0 ,Set p0 = r0 and q0 = s0 = Ar0 ,For k = 0, 1, . . . , until ∥rk ∥2 ≤ εTOL ∥b∥2 do :(qk , rk )αk =,(qk , qk )xk+1 = xk + αk pk ,rk+1 = rk − αk qk ,sk+1 = Ark+1 ,(qi , sk+1 )βk,i = −, (i = 0, 1, . . . , k)(qi , qi )kpk+1 = rk+1 + ∑ βk,i pi ,i=0k

qk+1 = sk+1 + ∑ βk,i qi ,End Fori=0・行列ベクトル積は1反復あたり1回・長い漸化式を使うため多くのメモリが必要
・リスタートにより演算量とメモリ量を削減- 17 -CCS HPC seminar 2024

反復法の収束特性

相対残差, ∥rk ∥2 /∥b∥2

行列 A：poisson3Da（n = 13,514，nnz(A) = 352,762, SuiteSparse Matrix Collection）
右辺 b ：b = [1, 1, . . . , 1]T

反復回数, k反復法の相対残差履歴■：BiCG法，■：CGS法，■：BiCGSTAB法，■：GMRES法- 18 -CCS HPC seminar 2024

複素対称行列に対する解法
TH

3. 係数行列が複素対称行列（ A = A ! A ）の場合
共役直交共役勾配法（Conjugate Orthogonal Conjugate Gradient: COCG法）

係数行列が複素対称行列の場合は，

補足：複素対称行列

1反復あたり1回の行列ベクトル積で，かつ短い漸化式で計算ができる

A = A T ! AH
(ai j = a ji ! ā ji )- 19 -CCS HPC seminar 2024

COCG法のアルゴリズム
x0 is an initial guess,
Compute r0 = b − Ax0 ,Set p0 = r0 ,For k = 0, 1, . . . , until ∥rk ∥2 ≤ εTOL ∥b∥2 do :qk = Apk ,行列ベクトル積(r̄k , rk )
αk =
,
内積
( p̄k , qk )xk+1 = xk + αk pk ,ベクトルの定数倍と加算
rk+1 = rk − αk qk ,(r̄k+1 , rk+1 )
βk =
,
内積
(r̄k , rk )pk+1 = rk+1 + βk pk ,ベクトルの定数倍と加算
End For- 20 -CCS HPC seminar 2024

疎行列が現れる例
2次元 Poisson 問題

⎧ 2
2
⎪
∂
u
∂
u
⎪
⎪
⎨
+ 2 = f , in Ω
2
⎪
∂x
∂y
⎪
⎪
⎩ u = ū,on ∂Ω
f, ū は既知の関数．
W をx, y 方向に (M+1) 等分し5点中心差分で離散化

y
1

O

1

x

M M次行列をもつ連立一次方程式に帰着
行列の要素数：M4 個
非零要素の数：5M2 - 4M 個
- 21 -CCS HPC seminar 2024

疎行列の格納形式（CRS形式）
Compressed Row Storage (CRS) 形式【補足】Compressed Sparse Row (CSR) と呼ぶこともあります．
val: 非零要素を格納する配列⎡⎤⎢⎢⎢ a11 0 a13 0 a15 ⎥⎥⎥⎥⎥⎢⎢⎢ 0 a22 0 a24 a25 ⎥⎥⎥⎥⎢⎢⎢A = ⎢⎢⎢⎢ a31 a32 a33 0 0 ⎥⎥⎥⎥⎥⎥⎥⎢⎢⎢ 0 0 a a
0
43
44
⎥⎦
⎢⎣
0 a52 0 a54 a55

col_ind: 非零要素の列番号を格納
row_ptr: 各行の先頭の非零要素が格納されている場所を指す配列

val: a11 a13 a15 a22 a24 a25 a31 a32 a33 a43 a44 a52 a54 a55
col_ind:

1

3

5

2

4

row_ptr:

1

4

7

10

12

5

1

2

3

3

4

2

4

5

15

最後は非零要素数+1 の値- 22 -CCS HPC seminar 2024

疎行列の格納形式（CCS形式）
Compressed Column Storage (CCS) 形式【補足】Compressed Sparse Column (CSC) と呼ぶこともあります．
val: 非零要素を格納する配列⎡⎤⎢⎢⎢ a11 0 a13 0 a15 ⎥⎥⎥⎥⎥⎢⎢⎢ 0 a22 0 a24 a25 ⎥⎥⎥⎥⎢⎢⎢A = ⎢⎢⎢⎢ a31 a32 a33 0 0 ⎥⎥⎥⎥⎥⎥⎥⎢⎢⎢ 0 0 a a
0
43
44
⎥⎦
⎢⎣
0 a52 0 a54 a55

row_ind: 非零要素の行番号を格納
col_ptr: 各列の先頭の非零要素が格納されている場所を指す配列

val: a11 a31 a22 a32 a52 a13 a33 a43 a24 a44 a54 a15 a25 a55
row_ind:

1

3

2

3

5

1

col_ptr:

1

3

6

9

12

15

3

4

2

4

5

1

2

5

最後は非零要素数+1 の値- 23 -CCS HPC seminar 2024

CRS形式の行列ベクトル積
CRS形式では，行列 A の非零要素は行方向に連続にアクセスされる．即ち，y の i 番目の成分 yi は，行列 A の i 行目とベクトル x の内積で求められる．

行列 A とベクトル x の積 y = Ax
⎡
⎤ ⎡
⎤⎡
⎤
⎢⎢⎢ y1 ⎥⎥⎥ ⎢⎢⎢ a11 a12 . . . a1n ⎥⎥⎥ ⎢⎢⎢ x1 ⎥⎥⎥⎢⎢⎢ y ⎥⎥⎥ ⎢⎢⎢ a⎥⎢⎥⎢⎢⎢ 2 ⎥⎥⎥ ⎢⎢⎢ 21 a22 . . . a2n ⎥⎥⎥⎥⎥ ⎢⎢⎢⎢⎢ x2 ⎥⎥⎥⎥⎥⎢⎢⎢ .. ⎥⎥⎥ = ⎢⎢⎢ ..
..
.. ⎥⎥⎥ ⎢⎢⎢ .. ⎥⎥⎥..
⎥⎥⎥ ⎢⎢⎢ . ⎥⎥⎥⎢⎢⎢ . ⎥⎥⎥ ⎢⎢⎢ .
.
.
.
⎣
⎦ ⎣
⎦⎣
⎦
yn
an1 an2 . . . ann
xn

Fortran のコード例
do i=1,ny(i) = 0.0D0
do j=row_ptr(i), row_ptr(i+1)-1
y(i) = y(i)+val(j)*x(col_ind(j))
end doend do- 24 -CCS HPC seminar 2024

CCS形式の行列ベクトル積
CCS形式では，行列 A の非零要素は列方向に連続にアクセスされる．即ち，行列ベクトル積 y = Ax は，行列 A の列ベクトルの線型結合として求められる．
行列 A とベクトル x の積 y = Ax⎡⎤x⎢⎢⎢ 1 ⎥⎥⎥⎢⎢⎢ x ⎥⎥⎥n⎢⎢⎢ 2 ⎥⎥⎥
y = [a1 , a2 , . . . , an ] ⎢⎢⎢ . ⎥⎥⎥ = ∑ ai xi⎢⎢⎢ .. ⎥⎥⎥ i=1⎣⎦
xn

Fortran のコード例
do i=1,ny(i) = 0.0D0
end dodo j=1,ndo i=col_ptr(j),col_ptr(j+1)-1
k = row_ind(i)y(k) = y(k) + val(i) * x(j)end doend do- 25 -CCS HPC seminar 2024

CRS形式の転置行列ベクトル積
T

CRS形式の転置行列ベクトル積 z = A x<latexit sha1_base64="DA9HTDPd+cnSdnod+R+B3K3Klxg=">AAALt3icnVZJb9NAFH5lJ2EpcEFCSIa2iEMbnFKWIlVqxYUjWxepoZW3tKPasbGddIly4sYf4MAJJA4VN/gJXPgDHPgJiCNIXDjwzRu7SdrEcfHI9vjN9/Y3z2MGrohiXf8+dOToseMnTp46XSieOXvu/PCFiwuRXw8tZ97yXT9cMo3IcUXNmY9F7DpLQegYnuk6i+bGQ7m+2HDCSPi15/F24LzwjLWaqArLiEFaHb5aMb3mTkub0eZWmhXPiNdDr/m81ZLkrdbq8Ihe0vnSDk7KyWSEkuuxf+H4HFXIJp8sqpNHDtUoxtwlgyKMZSqTTgFoL6gJWoiZ4HWHWlQAbx0oBwgD1A081/C1nFBr+JYyI+a2oMXFHYJTozH9m76r/9K/6h/1H/rfvrKaLEPaso23qXidYPX868vP/gzk8vCOab3NlWlzTFW6z7YK2B4wRXphKf7Gzptfzx48HWve0N/rP2H/O/27/gUe1Bq/rQ9PnKdvIX0sufvb5WImoyDYMi8jjr3Qgq3pzxHTFrzwOY8RI8cysCLJT5ZEM7GxP6Kd2Rhod09rGv0GUCE4BdslJUq8nak1wLuKp5QbH7Cwe82FZBMaZH1ug7eKbxfDBDamTTxlFrO0KQtNtlcMyIlEbNBOLswE3YafqaV5OATXmwuUnbnHmnjWsJLKnmAeA+/D2OeyvrUBFSAzVwdOcHV1IjvldGYggE0R41XOlZY8nB5jba6dIBfHS6z7vFMjGgdVanIPVHU/7gj1KTWl3CGGj5qRXxavSXvyybJ4v1o5IyRrLea6yytf8K5We8mBnjS60tZ2p8trp8NxNrgfD+axWaPPlLbeEusNwdGrjvrJqmJmH4qjM0ulxH5Z9/UDvSFbgp3ktsR9qp7T99sDe9VW0pPSaBymW7W7p6r51J/9fMpChfeS+tFoBncNK5scm7THKIpCerwnbHBWkqhLlPxScw2za32eKTpLZuO/ZObVmuYk9V72BxtD1m97t95CRqX320yfAa4ASd2jVwXPsMTuCh7vwZuOFCujruK/v9NpOC+V6A7em8k/W+O4tG1t0ghOVS0+U/TS1epJVXpirrxO67u9Un+MtqTCgNx5iFiEDhpyvJdpkk8/FZYmOOZlPv1V4If8a29hdZRWGKNOV3IfaolHo+xT/7ymmmQEJjkC8lRR4FPrNF+amtybSibT5b1T68JkqXy3NPVkamR2Mjm/nqIrdJ1uQvc9mqVH9Jjm4eMr2qVP9Lk4XVwtVovrCnpkKOG5RF1X8eU/rxtVeA==</latexit><latexit sha1_base64="lHVwHQOAfhYzye+XVzwNFYTVhqM=">AAAQaHicxVZPb9xEFH8tXSjmTzZwAMTFNE3FIV28Xm/bIEVqNyHqAaT+SdqKOF35z2xixfYutjfZjeUvwBfgwAkkDoiPwYULRw79CMCxSFw48OZ5nN3t2l6nQmKs9Yzn/X7vvXnvzeyYA9cJI0V5duHiK5dqr752+XXpjTffenupvvzOo7A/DCy2a/XdfvDENELmOj7bjZzIZU8GATM802WPzaNNLn98zILQ6fs70XjA9j3jwHd6jmVEONVdvvSrbrMekklV7BlR4IyaSXx354vPk3h7u9PsqImUg1EnmM82O+1cTCvDdNTtra2tXIx2htE67e3NXEw7w2zeWr+x3UmkVUk3vfg0kTckfejbLDADw2KxpLusF+2hkB04fmwEgTFOYou3RJJ1i7nu7DKNbtxcaybytXmhmihy2nKELc5sFTC1Mmb7TKjruS6VcFVuVV1Tc60Kl/KFC1xCZpszC1yiKLVeOkr5TI0LtQLhfxAlrSBKJUzhUj5TREkrj1K7IErkUr5wgUuLE8f3HvNtUe6SHjgHh9F+0o3vPI11RB6GvXgnSRL5/2pV9ijuz1G3SbEddVXRt0Svib7N+/y18vNglCTd+orSUKjJ84OmGKyAaPf6y7U7oIMNfbBgCB4w8CHCsQsGhPjsQRMUGODcPsQ4F+DIITmDBCTkDhHFEGHg7BG+D/BrT8z6+M11hsS20IqLvwCZMqwqvyk/Ks+VX5SflN+Vfwp1xaSD+zLG3ky5bNBd+vr9h38vZHnYR3A4YZX6HEEPbpGvDvo+oBm+CivlH59+8/zhpw9W42vK98qf6P93yjPlZ1yBf/yX9cN99uBb1L4qfsV+uTjiUXDIM68kjnloh7wpZkQwwlX0KY8hIVdLsI7IT5lGU/hYjJhkNkK0e2Y1i/4xogJkOuQX18jxdqnVAfY9fHO90ZyHszIXNZtogdfnGLk9/HbxMREbwQm+eRbLrKUemuSvsyAnHHEEp5Uw16GF68w8rcJwqN5cRNmleyzGt4+STPd14hjYn8c/l+wdLKgAnrkh4hyqrmnktJ7pDAzQp5Dwac5TK1WYHmFtqp1BJcZXKO/TTg1hDWe5JXeuqovYIdYnt5SxA3z6WDP8yyIZ96eaLov2q1UxQrzWIqq7qvod2tXpXmJoJ4su93Vy0lX1k1GcDTqPF3NsstinmYndBtkNkJFXR0W6ejiyz8WYzlJD+M/rfjh3NpRrsEVuG3RODSuuvbXwrBqJMymLxnlOq8npmdZ8tp4XeamHKd4T9SPDBv58lJxQbLIzJp1JkR7tCRuZuog6R/GvdCzj6KOCd4Yu03n8UjqrWs1ykq2enw82Prx+J7v1E8woX/2Y5jcQJ6Gm2SevgjdI42wFr+VwsyfD8qin8X/xpJPxvtSANvYn4j9bprhMfI1hBW9VCd0p8mwlubOpnYgqb9r72VWl/xgTTdKC3HkYsRBP0IDivQcq3X500uZQzJt0+9NxHfxfe4TSq/CUMOntiu9DWazoKq2pOK+ZJR4BlSLAbxUS3VrXqcnp4KYmBuvNs1vrI7XRvNHQ7msrt1Vxf70MH8IV+Bht34TbcBfuwS5YNa32Zc2q2Ut/1Ov19+ofpNCLFwTnXZhp9Sv/Alzv3G0=</latexit>

2
666 a1,1666 0666
z = 6666 a1,3
666 0
64
a1,5
|

0
a2,2
0
a2,4
a2,5

a3,1
a3,2
a3,3
0
0
{z
AT

0
0
a4,3
a4,4
0

32
3
0 77 66 x1 77
76
7
a5,2 7777 6666 x2 7777
76
7
0 7777 6666 x3 7777
76
7
a5,4 77775 66664 x4 77775
a5,5x5} | {z }x

Fortran のコード例
do i=1,nz(i) = 0.0D0
end dodo j=1,ndo i=row_ptr(j),row_ptr(j+1)-1
k = col_ind(j)z(k) = z(k) + val(i) * x(i)end doend do• CRS形式で格納される行列の非零要素は，行方向に連続にアクセスされる．
• この行列の転置をとると，非零要素は列方向に連絡にアクセスされることになる．
• 即ち，CRS形式で格納された行列の転置行列ベクトル積を行うには，CCS形式の行列ベクトル積コードを実行すればよい，ということになる．
- 26 -CCS HPC seminar 2024

CCS形式の転置行列ベクトル積
T

CCS形式の転置行列ベクトル積 z = A x<latexit sha1_base64="DA9HTDPd+cnSdnod+R+B3K3Klxg=">AAALt3icnVZJb9NAFH5lJ2EpcEFCSIa2iEMbnFKWIlVqxYUjWxepoZW3tKPasbGddIly4sYf4MAJJA4VN/gJXPgDHPgJiCNIXDjwzRu7SdrEcfHI9vjN9/Y3z2MGrohiXf8+dOToseMnTp46XSieOXvu/PCFiwuRXw8tZ97yXT9cMo3IcUXNmY9F7DpLQegYnuk6i+bGQ7m+2HDCSPi15/F24LzwjLWaqArLiEFaHb5aMb3mTkub0eZWmhXPiNdDr/m81ZLkrdbq8Ihe0vnSDk7KyWSEkuuxf+H4HFXIJp8sqpNHDtUoxtwlgyKMZSqTTgFoL6gJWoiZ4HWHWlQAbx0oBwgD1A081/C1nFBr+JYyI+a2oMXFHYJTozH9m76r/9K/6h/1H/rfvrKaLEPaso23qXidYPX868vP/gzk8vCOab3NlWlzTFW6z7YK2B4wRXphKf7Gzptfzx48HWve0N/rP2H/O/27/gUe1Bq/rQ9PnKdvIX0sufvb5WImoyDYMi8jjr3Qgq3pzxHTFrzwOY8RI8cysCLJT5ZEM7GxP6Kd2Rhod09rGv0GUCE4BdslJUq8nak1wLuKp5QbH7Cwe82FZBMaZH1ug7eKbxfDBDamTTxlFrO0KQtNtlcMyIlEbNBOLswE3YafqaV5OATXmwuUnbnHmnjWsJLKnmAeA+/D2OeyvrUBFSAzVwdOcHV1IjvldGYggE0R41XOlZY8nB5jba6dIBfHS6z7vFMjGgdVanIPVHU/7gj1KTWl3CGGj5qRXxavSXvyybJ4v1o5IyRrLea6yytf8K5We8mBnjS60tZ2p8trp8NxNrgfD+axWaPPlLbeEusNwdGrjvrJqmJmH4qjM0ulxH5Z9/UDvSFbgp3ktsR9qp7T99sDe9VW0pPSaBymW7W7p6r51J/9fMpChfeS+tFoBncNK5scm7THKIpCerwnbHBWkqhLlPxScw2za32eKTpLZuO/ZObVmuYk9V72BxtD1m97t95CRqX320yfAa4ASd2jVwXPsMTuCh7vwZuOFCujruK/v9NpOC+V6A7em8k/W+O4tG1t0ghOVS0+U/TS1epJVXpirrxO67u9Un+MtqTCgNx5iFiEDhpyvJdpkk8/FZYmOOZlPv1V4If8a29hdZRWGKNOV3IfaolHo+xT/7ymmmQEJjkC8lRR4FPrNF+amtybSibT5b1T68JkqXy3NPVkamR2Mjm/nqIrdJ1uQvc9mqVH9Jjm4eMr2qVP9Lk4XVwtVovrCnpkKOG5RF1X8eU/rxtVeA==</latexit><latexit sha1_base64="4mh9hYZxNcucb7o2+gThhau2/G4=">AAAP73icnVbNbttGEJ6kSpuyP4rbS4Fe2DoOenAUiqKSuICBRHaNHFogP3YSwHQE/qxkwiSlkJQtmeAL9AVyKFC0BXoI+gg99tIX6CHHHtscU6CXHjo7S1pSTFJ0SIhc7nzfzLezs6s1h64TRory4tz5t2oX3n7n4rvSe+9/8GH90tJHD8PBKLDYjjVwB8Fj0wiZ6/hsJ3Iilz0eBszwTJc9Mg82uP3RIQtCZ+BvR5Mh2/OMvu/0HMuIsKu7VNvWbdZDMrmKPSMKnHEzie9sf/N1Em9tdZodNZFyMOoU89VGp52LaWWYjrq1ubmZi9FOMFqnvbWRi2lnmI2ba9e3Oom0IummFx8n8rqkj3ybBWZgWCyWdJf1ol00sr7jx0YQGJMktviVSLJuMdedH6bRjZurzUS+kmdUZHHlGpHZejNmZtT1HElqGVflUdVVNTeqKiQVGRe5bXNmrqSWyFIr13GrzHFLSCpiolFb6DZXklYWVRNZ0nIdL2ZqhUzKklaQpbbIUjuX2xaSiowlkhYYSRJfe8y303KX9MDp70d7STe+/STWEbkf9uLtJEmqLBVcJuNuk4Y47qrpu5W+tfTd5u/8kHxZjpOke2lZaSh0yacbzbSxDOl1d7B04TboYMMALBiBBwx8iLDtggEh3rvQBAWG2LcHMfYF2HLIziABCbkjRDFEGNh7gM8+fu2mvT5+c58hsS2M4uIvQKYMK8ofynPllfK78ovyl/Jfoa+YfHAtE3ybgsuG3fq3nzz4dyHLw3cE+1NWqeYIenCTtDqofUg9fBSW4B8eP3v14Mv7K/EV5SflJer/UXmh/IYj8A//sX6+x+5/h95X0l+xLhdbPAsOKfNK8piHdkhNMSOCMY5iQPMYEnKlBOuk81Pm0Uw1FiOmMxsh2j2JmmX/EFEBMh3SxT1yvF0adYjvHj653+iUwnmbi55NjMDrc4LcHn67eJuIjeAIn3wWy6IJhSbpdRbMCUccwHElzFVo4TgzpVUYDtWbiyi7dI3F+PTRkvm+ShwD32fR51K8/oIK4DM3QpxD1TWLnPUzOwND1BQSXsy5iFKF6RHWptoZVmI8RfuAVmoIq9jLI7mnqrqIHWJ98kgZO8B7gDXDvyyycT3VfFm0Xq2KGeK1FlHdVfXv0KoWa4lhnCy7XOt0p6uqk1GeDdqPF3NsijignmncBsUNkJFXR0W+etiyz8SYnaVGqp/X/ejU3lDuwU7ntkH71Kji2FsL96pxuidl2TjLbjXdPUXNZ+N5nScUCryX1o8M6/jz0XJEucn2GNEjkB6tCRuZepp1juJfoi1j67OCZ4Yu83n4Rj6rRs3mJBs93x9svHn9TlfrNZxRPvoJ9a8jTkJP83deBa+Tx/kKXs3hZneG5VkX+X99p5PxvNSANr6P0v9smfIy1RrDMp6qEjpT5MVKcntFnIgqb1b9/KjEP8bUk7Rg7jzMWIg7aED53gWVTj86eXMo5006/ek4Dv6vPUbrZXhCGHG64utQTkd0mcZUPK9ZJJ4BlTLATxUSnVrX6JJF44aWNtaaJ6fWh2qjeb2h3dOWb6np+fUifAqfwxcY+wbcgjtwF3bAqv1a+7P2d+1l/Wn9Wf37+g8Cev5cyvkY5q768/8BNH/Qiw==</latexit>

2
666 a1,1666 0666
z = 6666 a1,3
666 0
64
a1,5
|

0
a2,2
0
a2,4
a2,5

a3,1
a3,2
a3,3
0
0
{z
AT

0
0
a4,3
a4,4
0

32
3
0 77 66 x1 77
76
7
a5,2 7777 6666 x2 7777
76
7
0 7777 6666 x3 7777
777 666
777
a5,4 775 664 x4 775
a5,5x5} | {z }Fortran のコード例
do i=1,nz(i) = 0.0D0
do j=row_ptr(i), row_ptr(i+1)-1z(i) = z(i)+val(j)*x(col_ind(j))
end doend dox• CRS形式では行列の非零要素は列方向に連続にアクセスされる，転置をとると非零要素は行方向に連絡にアクセスされることになる．
• 先程と同様に，CCS形式で格納された行列の転置行列ベクトル積を行うには，
CRS形式の行列ベクトル積コードを実行すればよい．
- 27 -CCS HPC seminar 2024

行列ベクトル積の並列化
・CRS形式の場合の y = Ax の計算
Proc. 0
Proc. 1
Proc. 2=*

Proc. 3Axy

各プロセスで分散してデータを保持する

全てのプロセス
で保持する

MPI_Gather でProc. 0 に集める
- 28 -CCS HPC seminar 2024

行列ベクトル積の並列化

Proc. 3

Proc. 2

Proc. 1

Proc. 0・CCS形式の場合の y = Ax の計算

=

*

+

+

+

A

x

y

各プロセスで分散してデータを保持する

分散して
保持する

MPI_Reduce でProc. 0に結果を足し合わせて送る- 29 -CCS HPC seminar 2024

Proc. 1=x jy jj=1

Proc. 2

Proc. 3=

Proc. 0=(x, y) =n!

=

内積の並列化
xy
tmp_sum tmp_sum

tmp_sum tmp_sum

sum
MPI_Reduce で Proc.0 に集める
- 30 -CCS HPC seminar 2024

内積計算の MPI コード例
n!
program main(x, y) =x jy jinclude 'mpif.h'j=1...
call mpi_init(ierr)call mpi_comm_size(mpi_comm_world, nprocs, ierr)call mpi_comm_rank(mpi_comm_world, myrank, ierr)...
tmp_sum = 0.0D0
do i=istart(myrank+1), iend(myrank+1)
tmp_sum = tmp_sum + x(i) * y(i)end docall mpi_reduce(tmp_sum, sum, 1, mpi_double_precision,
mpi_sum, 0, mpi_comm_world, ierr)...
call mpi_finalize(ierr)- 31 -CCS HPC seminar 2024

ベクトルの定数倍と加算の並列化
y = y + ax （x, y：ベクトル，a：スカラー）a を MPI_Bcast で全プロセスに送る
Proc. 0aProc. 1

Proc. 2

Proc. 3

a

a

a

a

+

+

+

+

x

y
- 32 -CCS HPC seminar 2024

複数本の右辺ベクトルをもつ
連立一次方程式の解法- 33 -CCS HPC seminar 2024

複数右辺ベクトルをもつ方程式
右辺が s 本の連立一次方程式
AX = B
ここで，A : n次行列，""!
!
(1) (2)(s)(1) (2)(s)
X = x , x ,...,x ,B = b , b ,...,b
直接法による解法
・係数行列の完全分解（A = LU など）が必要
・完全分解できれば，s 回の前進・後退代入で OK・分解には多くの計算量，メモリ量が必要
反復法で s 本の方程式を効率よく解けないか？
- 34 -CCS HPC seminar 2024

ブロッククリロフ部分空間反復法
ブロッククリロフ部分空間反復法の種類

・Block BiCG法
・Block GMRES法・Block QMR法・Block BiCGSTAB法
・Block BiCGGR法

O’Leary (1980)
Vital (1990)
Freund (1997)
El Guennouni (2003)
Tadano (2009)

複数本の右辺ベクトルをまとめて扱うことで効率よく解を求めることができる- 35 -CCS HPC seminar 2024

ブロッククリロフ部分空間反復法「効率よく」とはどういうことか？

相対残差, ∥Rk ∥F /∥B∥F

1本ずつ解いた場合より少ない反復回数で残差が収束する（ことがある）

反復回数, kBlock BiCGSTAB法の相対残差履歴．
ここで，■：s = 1，■：s = 2，■：s = 4．

- 36 -CCS HPC seminar 2024

Block CG法
X0 ∈ Cn×s is an initial guess,
Compute R0 = B − AX0 ,CG法と異なるところ

Set P0 = R0 ,For k = 0, 1, . . . , until ∥Rk ∥F ≤ εTOL ∥B∥F do:
Qk = APk ,H
Solve (PHk Qk )αk = Rk Rk for αk ,1. 行列・ベクトル積の本数が1本から s 本に増加．
2. αk , βk が s 次行列になった．

Xk+1 = Xk + Pk αk ,Rk+1 = Rk − Qk αk ,H
Solve (RH
k Rk )βk = Rk+1 Rk+1 for βk ,3. ベクトルの定数倍の計算が行列・行列積になった．

Pk+1 = Rk+1 + Pk βk ,End For- 37 -

行列・ベクトル積の効率化

CCS HPC seminar 2024

・行列は CRS 形式で格納されているとする．
・Y = AX を計算．Y, X は n 行 s 列の配列．
do k=1,s! s についてのループが一番外側
do i=1,ndo j=row_ptr(i), row_ptr(i+1)-1
Y(i,k) = Y(i,k) + val(j) * X(col_ind(j),k)end doend doend do[ 問題点 ]・X についてメモリの連続アクセスができていない（Fortran は行方向に連続にデータが並んでいる）
・行列データを s 回読まなければならない．

- 38 -

行列・ベクトル積の効率化

CCS HPC seminar 2024

[ 解決策 ]・X, Y を転置した形で保持させる
do i=1,ndo j=row_ptr(i), row_ptr(i+1)-1
do k=1,s! s についてのループが一番内側

Y(k,i)= Y (k,i) + val(j) * X(k,col_ind(j))
end doend doend do・X について（少なくとも s 回は）連続アクセスができる．
・行列データは1回しか読み込まない．
・Y についても連続アクセスが保たれている．

- 39 -CCS HPC seminar 2024

n s 行列・s s 行列積の計算
・行列・ベクトル積の効率化のためにベクトルを転置した．
・n s 行列と s s 行列の積の計算も工夫が必要．
転置Xk+1 = Xk + Pk αk

T
Xk+1= XkT + αTk PTk

do j=1,ndo i=1,s
do k=1,sX(k,j) = X(k,j) + Alpha(k,i) * P(i,j)
end doend doend doベクトルと αk を全て転置して保持することで連続アクセスが可能に．

Alpha は転置済みとする．
- 40 -

s n 行列・n s 行列積の計算

CCS HPC seminar 2024

・ak, bk を求めるために必要．
H・ Ck = Pk Qk を計算することを考える．
do j=1,ndo i=1,s
do k=1,sC(k,i) = C(k,i) + dconjg(P(k,j)) * Q(i,j)
end doend doend do- Ck の計算も，連続メモリアクセスを保持できる．
- dconjg：複素数の共役を求める関数．
- 41 -CCS HPC seminar 2024

OpenMP による並列化
・共有メモリ向けの並列化インターフェース．
・既存のプログラムに数行加えるだけで並列化ができる．
!$OMP PARALLEL【 プログラム 】!$OMP END PARALLELと書くと，スレッドが立ち上がり，スレッド毎に別々の処理ができるようになる．
（以下のコードは，!$OMP PARALLEL と!$OMP END PARALLEL で囲まれているとします．）- 42 -CCS HPC seminar 2024

OpenMP による並列化
1. 行列・ベクトル積の並列化!$OMP DO PRIVATE(j,k)
do i=1,ndo j=row_ptr(i), row_ptr(i+1)-1
do k=1,sY(k,i)=Y(k,i)+val(j)*X(k,col_ind(j))
end doend doend do!$OMP END DO

最初の do ループの前に !$OMP DO ... を加えるだけ．
最後の !$OMP END DO は省略可能．

- 43 -CCS HPC seminar 2024

OpenMP による並列化
2. n s行列・s s行列の積の並列化!$OMP DO PRIVATE(i,k)
do j=1,ndo i=1,s
do k=1,sX(k,j) = X(k,j) + Alpha(k,i) * P(i,j)
end doend doend do!$OMP END DO

一番外の do ループの前に !$OMP DO ... を加えるだけ．
最後の !$OMP END DO は省略可能．
- 44 -CCS HPC seminar 2024

OpenMP による並列化
3. s n行列と n s行列の積の並列化
H
以下を実行して Ck = Pk Qk を計算する．
!$OMP SINGLE
do j=1,sdo i=1,sC(i,j) = 0.0D0
end doend do!$OMP END SINGLE!$OMP DO PRIVATE(i,k) REDUCTION(+:C)
do j=1,ndo i=1,s
do k=1,sC(k,i) = C(k,i) + dconjg(P(k,j)) * Q(i,j)
end doend doend do!$OMP END DO- 45 -CCS HPC seminar 2024

実行時間 / s [秒]行列ベクトル積の性能
改良版
単純な実装

右辺ベクトル数, s

右辺ベクトル数の変化に対する計算時間変化（OpenMP 12並列）．
・行列：格子量子色力学計算（QCD）で現れる行列
・行列サイズ：1,572,864，非ゼロ要素数：80,216,064．
・実験環境：CPU : Intel Xeon E5-2620v3 2.4GHz 2，
・コンパイラ：gfortran ver. 5.4，オプション：-O3 –fopenmp- 46 -CCS HPC seminar 2024

まとめ・連立一次方程式の解法である
クリロフ部分空間反復法を取り上げた．
・疎行列に対する行列・ベクトル積の実装
方法とその並列化について述べた．
・ブロッククリロフ部分空間反復法とコードの最適化，及び OpenMP での並列化について述べた．

- 47 -CCS HPC seminar 2024

レポート課題
次の行列 A と右辺ベクトル b をもつ連立一次方程式 Ax = b を，BiCG法とBiCGSTAB法
で解くプログラムを作りなさい．但し，行列 A はCRS形式で格納すること．
行列サイズ n は任意とする．
⎡
⎢⎢⎢ 2
⎢⎢⎢ γ
⎢⎢⎢
⎢
A = ⎢⎢⎢⎢
⎢⎢⎢
⎢⎢⎢
⎣

0

1
2
..
.

1
..
.
γ

0
..

.
2
γ

⎡ ⎤
⎤
⎢⎢⎢ 1 ⎥⎥⎥⎥⎥⎥⎢⎢⎢ 1 ⎥⎥⎥⎥⎥⎥⎢⎢⎢ ⎥⎥⎥
⎥⎥⎥
⎢ . ⎥
⎥⎥⎥
⎥⎥⎥ , b = A ⎢⎢⎢⎢⎢ .. ⎥⎥⎥⎥⎥⎢⎢⎢ ⎥⎥⎥⎥
1 ⎥⎥⎥⎥⎦⎢⎢⎣ 1 ⎥⎥⎦
1
2

12
• 行列 A のパラメータ は 0 < < 1 とする．反復は条件：krk k2 /kbk2  10 を満たしたら停止すること．初期解 x0 はゼロベクトルとする．
• 複数のパラメータ について実験し，反復過程における相対残差 krk k2 /kbk2 をグラフにプロットするとともに考察を述べなさい．
• プログラムリストも提出すること．プログラミング言語は何を用いてもよい．
<latexit sha1_base64="NHraSC1NY5jp8aEFEdp7VfQeWUg=">AAALyXicnVa9b9tGFH/5bCq1iZMuBbowtV10iOWT4iRNAQMpsgTIksS1HcByBPFD7sGkSJOU/KFoypY96JApBToE+RM6tgX6D3Twn1B0dIEuHfq7d0dLiiWKDg8kj+9+7/vd49mRL5NUiMMzZ8+dv3Dxo0sfl8qffHr5yszVa2tJ2Ikdb9UJ/TB+ajcTz5dtbzWVqe89jWKvGdi+t25v31fr610vTmTY/j7dj7zNoLnVli3pNFOQGjOL9edW3Q56cb+xbdWfN2rWomVIdl8T6r63Y1XFs95CtdZvzMyKiuDLOjmpmsksmetRePXCd1Qnl0JyqEMBedSmFHOfmpRgbFCVBEWgbVIPtBgzyese9akE3g5QHhBNULfx3MLXhqG28a1kJsztQIuPOwanRfPiT/FWHIk/xDvxl/hvoqwey1C27ONta14valx5+fnKv1O5ArxT+mHAlWtzSi36hm2VsD1iivLC0fzdgx+PVr59Mt/7Svwk/ob9b8Sh+BUetLv/OD8/9p68hvR5c0+2y8dMRUGyZUFOHMehJVszmSOlPXgRch4TRs7nYKXJT55E29g4GTHIbAq0f6w1i34XqBicku1SEhXezdUa4d3CU8lNT1g4uuZDsg0Nqj73wdvCt49hA5vSLp4qi3natIU22yun5EQhtumgEGaBbsLPzNIiHJLrzQfKzd1jPTzbWMlkLzBPE+/T2Oezvq0pFaAy1wFOcnUNI4flDGcggk0J43XOtZYinAFjXa6dqBDHDtZD3qkJ3QBVafJPVPUk7gT1qTRl3DFGiJpRXw6vKXuKyXJ4vzoFI6RqLeW6Kypf8q7We8mDniy6ytZBpytqp8dxbnI/ns7jssaQKQO9FdYbg2NcHU2S1cLMPRXHcJYqxn5V950TvSFfgmtyW+E+1Sno+82pvWrP9KQsGqfpVoPuqWs+8+d9Pm2hxgemfixaxt3Gyi7HJusxmqKRAe8JF5x1E3WFUl96bmF2fcIzQ+fJ7H6QzKJas5xk3qv+4GKo+h3s1kVkVHm/z/Rl4EqQNDrGVfAySxyt4BtjeLORYVXUdfzf73QWzksVuoX3rvlnWxyXga09msWpqs9ninG6+mOpWk/KlTds/ahX+o8xkFSakrsAEUvQQWOO9wbV+PRTZ2mSY17l018dfqi/9h5W5+gZY/TpSu1Dy3g0xz5NzmumSUWgxhFQp4oSn1rv8mXpyZ0lM7lbPT61rtUq1duVpcdLs/dq5vx6ib6gL+lr6L5D9+gBPaJV+PiKfqHf6Pfyw/JOea98oKFnzxiez2jkKr/4H4u9Wik=</latexit><latexit sha1_base64="IeXBZTxMk5bEeo9705P8HLiK+tE=">AAALpnicnVbLbtNAFL28S8KbDRIbQwliAWFSCgVEJRAbVohXW6SmgF8po9qxsZ20IeoP8AMgsWolFojPYMMPsOATEEuQ2LDgzB27SWjiuHhke3zn3Ped67FCT8aJEN927d6zd9/+AxMHS+VDh48cPXb8xHwctCLbnbMDL4ieWmbserLpziUy8dynYeSavuW5C9bKXbW+0HajWAbNJ0kndJd8c7kpG9I2E5CeCeOWUV82fd/EpPb82KSoCr6M7ZNaOpmk9HoQHN93h+rkUEA2tcgnl5qUYO6RSTHGItVIUAjaEnVBizCTvO7SOpXA2wLKBcIEdQXPZXwtptQmvpXMmLltaPFwR+A0qCK+io/ip/giPonv4s9IWV2WoWzp4G1pXjd8fvTNqce/x3L5eCf0sseVa3NCDbrOtkrYHjJFeWFr/vbrtz8f33xU6Z4Xm+IH7N8Q38RneNBs/7I/PHQfvYf0SnqPtsvDTEVBsmV+ThyHoSVbM5ojoTV4EXAeY0ZWcrAyzU+eRCu1cTSil9kEaG9Laxb9NlAROCXbpSQqvJOrNcS7gaeSm2yzcHDNg2QLGlR9dsDbwLeHYQGb0CqeKot52rSFFtsrx+REIVbodSHMJboCPzNLi3BIrjcPKCd3j3XxbGIlk32JeUy8d2Kfx/qWx1SAylwLOMnV1Y/sl9OfgRA2xYzXOddainD6jHW4dsJCHK+wHvBOjekiqEqTt62qR3HHqE+lKeOOMALUjPqyeU3ZU0yWzfvVLhghVWsJ111R+ZJ3td5LLvRk0VW29jpdUTtdjrPJ/Xg8j8MaA6b09FZZbwSOYXU0SlYDM2dHHP1Zqqb2q7pvbesN+RKcNLdV7lOtgr5fGdur1tKelEVjJ92q1z11zWf+/MunLdR4P60fg2ZxN7GyyrHJeoymaKTPe8IBZz2NukKpLz03MDsz4pmh82S2/0tmUa1ZTjLvVX9wMFT99nbrZWRUed9h+ixwJUgaHMMqeJYlDlbwxSG82ciwKuo6/v92OgPnpSpdxXs1/WcbHJeerV2axKlqnc8Uw3StD6VqPQlXXr/1g17pP0ZPUmlM7nxELEYHjTjeizTFp586S5Mc8xqf/urwQ/2117B6jp4xRp+u1D40Uo/OsU+j85ppUhGY4gioU0WJT603+DL0ZGY6ndyobZ1a56eqtWvV6YfTk7en0vPrBJ2ms3QBumfoNt2jBzQHHyN6Rxu0Wb5Qvl+eKy9o6O5dKc9JGrjKL/4Cs1xM7A==</latexit><latexit sha1_base64="ip2HZYk9xOvEY/hJWaiMOl7ZqxM=">AAALnnicnVbLbtNAFL3lWRKgLWyQ2BhKEQsanLZQilQJhITYIB6lpVJTKr9ShtqxsZ30EfUHWLFDwAokFojPYMMPsOATEMsisWHBmTs2SdrEcfHI9vjOue8712MGrohiXf8+cODgocNHjg4eKxSPnzg5NDxyaiHy66HlzFu+64eLphE5rqg587GIXWcxCB3DM13nibl2W64/aThhJPza43gzcJY9Y7UmqsIyYpAWKquG5xkrw6N6SedL2zspJ5NRSq4H/sjhW1Qhm3yyqE4eOVSjGHOXDIowlqhMOgWgLVMTtBAzwesObVMBvHWgHCAMUNfwXMXXUkKt4VvKjJjbghYXdwhOjcb0b/onfUf/qn/Wf+h/espqsgxpyybepuJ1gpWhl2fmfvfl8vCO6VmLK9PmmKp0nW0VsD1givTCUvyNrdc7czcejTUv6h/0n7D/vf5d/wIPao1f1seHzqN3kD6W3L3tcjGTURBsmZcRx25owdb05ohpA174nMeIkWMZWJHkJ0uimdjYG9HKbAy0+09rGv0GUCE4BdslJUq8nak1wLuKp5Qb77Gwc82FZBMaZH1ugreKbxfDBDamdTxlFrO0KQtNtlf0yYlErNFWLsw4TcLP1NI8HILrzQXKztxjTTxrWElljzOPgfd+7HNZ32qfCpCZqwMnuLrake1y2jMQwKaI8SrnSkseTo+xNtdOkIvjBdZ93qkRXQZVanL3VHUv7gj1KTWl3CGGj5qRXxavSXvyybJ4v1o5IyRrLea6yytf8K5We8mBnjS60tZWp8trp8NxNrgf9+exWaPPlJbeEusNwdGtjnrJqmJm74ujPUulxH5Z9/U9vSFbgp3ktsR9qp7T98m+vWoj6UlpNPbTrVrdU9V86s9uPmWhwntJ/Wg0i7uGlXWOTdpjFEUhPd4TNjgrSdQlSn6puYbZuR7PFJ0ls/FfMvNqTXOSei/7g40h67e1W68go9L7TabPAleApM7RrYJnWWJnBV/uwpuOFCujruK/u9NpOC+V6Cre68k/W+O4tGxt0ihOVdt8puima7srVemJufLare/0Sv0xWpIKfXLnIWIROmjI8V6iCT79VFia4JiX+fRXgR/yr72B1Qv0lDHqdCX3oZZ4dIF96p3XVJOMwARHQJ4qCnxqneFLU5PpqWQyU/53al2YKJWvlaYeTo3enEjOr4N0ls7TJeieppt0lx7QPHx8Tq/oDb0tasU7xXvF+wp6YCDhOU0dV3HxLxCXS0M=</latexit><latexit sha1_base64="xu1hIYeTsipgotNzmEsjuvuhd34=">AAALu3icnVZJb9NAFH7sJWEpcEHi4lKKOEBw0rJKlUBcOJalgNRAFC8pQ+3Y2E66pDkj8Qc4cAKJA0Jc4Qdw4Q9w4CcgjiBx4cA3b8YkoYnj4pHt8Zvv7W+exwo9ESem+XXHzl279+zdN7G/UDxw8NDhySNH78VBK7LdRTvwguiBVY9dTzTdxUQknvsgjNy6b3nufWvlhly/33ajWATNu8l66D7068tN0RB2PQGpNjlV3TSqlt+JurUVo7pZqxjnDU2yukyoTU6bJZMvY+ukrCfTpK+F4Mie61QlhwKyqUU+udSkBHOP6hRjLFGZTApBe0gd0CLMBK+71KUCeFtAuUDUQV3BcxlfS5raxLeUGTO3DS0e7gicBs2YX8y35g/zs/nO/Gb+HimrwzKkLet4W4rXDWuHnx+/82ssl493Qo97XJk2J9Sgy2yrgO0hU6QXtuJvb7z4cefq7ZnOafO1+R32vzK/mp/gQbP9035zy739EtJn9D3aLg8zGQXBlvkZcRyGFmzNaI6E1uBFwHmMGTmTgRU6P1kSLW3jaEQvswnQ3l+tafTbQEXgFGyXlCjxTqbWEO8GnlJussXCwTUPki1okPW5Dt4Gvj0MC9iEVvGUWczSpiy02F4xJicSsUIbuTDnaBZ+ppbm4RBcbx5QTuYe6+DZxEoq+xzz1PHejn0e61seUwEycy3gBFdXP7JfTn8GQtgUM17lXGnJw+kz1uHaCXNxPMV6wDs1prOgSk3elqoexR2jPqWmlDvCCFAz8svmNWlPPlk271c7Z4RkrSVcd3nlC97Vai+50JNGV9ra63R57XQ5znXux+N5HNYYMKWnt8R6I3AMq6NRshqYOdvi6M9SSdsv6761pTdkS3B0bkvcp1o5fZ8d26vWdE9Ko7GdbtXrnqrmU3/+5VMWKryv68egedxNrKxybNIeoygK6fOecMBZ1VGXKPml5gZmUyOeKTpLZvu/ZObVmuYk9V72BwdD1m9vt55HRqX360yfB64ASYNjWAXPs8TBCj47hDcdKVZGXcX/305n4LxUogt4r+p/tsFx6dnaoWmcqrp8phimqzuUqvQkXHn91g96pf4YPUmFMbnzEbEYHTTieC9RhU8/VZYmOOZlPv1V4Yf8a69h9RQ9Yow6Xcl9aGiPTrFPo/OaapIRqHAE5KmiwKfWK3wZanJpTk+ulP+eWu9VSuWLpblbc9PXKvr8OkEn6CSdge5LdI1u0gItwsdn9J4+0MfifNEuPil6Crpzh+Y5RgNXsfUHyp9VeA==</latexit><latexit sha1_base64="ip2HZYk9xOvEY/hJWaiMOl7ZqxM=">AAALnnicnVbLbtNAFL3lWRKgLWyQ2BhKEQsanLZQilQJhITYIB6lpVJTKr9ShtqxsZ30EfUHWLFDwAokFojPYMMPsOATEMsisWHBmTs2SdrEcfHI9vjOue8712MGrohiXf8+cODgocNHjg4eKxSPnzg5NDxyaiHy66HlzFu+64eLphE5rqg587GIXWcxCB3DM13nibl2W64/aThhJPza43gzcJY9Y7UmqsIyYpAWKquG5xkrw6N6SedL2zspJ5NRSq4H/sjhW1Qhm3yyqE4eOVSjGHOXDIowlqhMOgWgLVMTtBAzwesObVMBvHWgHCAMUNfwXMXXUkKt4VvKjJjbghYXdwhOjcb0b/onfUf/qn/Wf+h/espqsgxpyybepuJ1gpWhl2fmfvfl8vCO6VmLK9PmmKp0nW0VsD1givTCUvyNrdc7czcejTUv6h/0n7D/vf5d/wIPao1f1seHzqN3kD6W3L3tcjGTURBsmZcRx25owdb05ohpA174nMeIkWMZWJHkJ0uimdjYG9HKbAy0+09rGv0GUCE4BdslJUq8nak1wLuKp5Qb77Gwc82FZBMaZH1ugreKbxfDBDamdTxlFrO0KQtNtlf0yYlErNFWLsw4TcLP1NI8HILrzQXKztxjTTxrWElljzOPgfd+7HNZ32qfCpCZqwMnuLrake1y2jMQwKaI8SrnSkseTo+xNtdOkIvjBdZ93qkRXQZVanL3VHUv7gj1KTWl3CGGj5qRXxavSXvyybJ4v1o5IyRrLea6yytf8K5We8mBnjS60tZWp8trp8NxNrgf9+exWaPPlJbeEusNwdGtjnrJqmJm74ujPUulxH5Z9/U9vSFbgp3ktsR9qp7T98m+vWoj6UlpNPbTrVrdU9V86s9uPmWhwntJ/Wg0i7uGlXWOTdpjFEUhPd4TNjgrSdQlSn6puYbZuR7PFJ0ls/FfMvNqTXOSei/7g40h67e1W68go9L7TabPAleApM7RrYJnWWJnBV/uwpuOFCujruK/u9NpOC+V6Cre68k/W+O4tGxt0ihOVdt8puima7srVemJufLare/0Sv0xWpIKfXLnIWIROmjI8V6iCT79VFia4JiX+fRXgR/yr72B1Qv0lDHqdCX3oZZ4dIF96p3XVJOMwARHQJ4qCnxqneFLU5PpqWQyU/53al2YKJWvlaYeTo3enEjOr4N0ls7TJeieppt0lx7QPHx8Tq/oDb0tasU7xXvF+wp6YCDhOU0dV3HxLxCXS0M=</latexit>- 48 -