筑波大学のスパコン「T2K-Tsukuba」で約2時間36分かけて5×5の魔方陣（≠魔法陣）の全解を求めたというニュースがありました。
「スパコン」ではなく「パソコン」だったらどうなのだろう、と思って試してみたら、10分でした（Core i7 4930K）。
コードはこちら（頭の悪そうなコードですが、スクリプトで枠組みを出力して手直して作れれば、そんなに大変ではないかと）。Visual Studio 2013 ExpressとGCC 4.6.3、Clang 3.4で動作を確認しています（このコードは、GCCとClangでは時間をちゃんと計れないので、コマンドtimeなどを使ってください）。ClangはOpenMPに対応させるのが面倒で、ちょっと試したところでは、GCCより遅かったです。
多重for文ですべての場合を探索するというのが基本方針ですが、それではさすがに終わらないので、次のような工夫をします。
- 各行・各列・対角線の和は65。
- （カンで）効率が良くなりそうなところから順番（下図の番号順）埋めていく。
- 回転と裏返しで同一になるものを重複して数えないように、「左上<右上<左下」かつ「左上<右下」のような制約を入れる。左上は22以下としてよい。
- a+b+c+d+eのうち、a,b,c,dが決まったら、e=65-a-b-c-dになる（下図の紫部分）。
- a+b+c+d+eのうち、a,b,cが決まったら、dはmax(1,65-a-b-c-25)からmin(25,65-a-b-c-1)の範囲で調べればよい。
| 6 | 10 | 11 | 12 | 8 |
| 20 | 1 | 18 | 2 | 21 |
| 22 | 13 | 5 | 16 | 23 |
| 24 | 3 | 19 | 4 | 25 |
| 9 | 14 | 17 | 15 | 7 |
魔方陣の解の列挙は並列化しやすい問題ですが、ここでの方針では、探索効率を上げるためには条件分岐が不可欠なので、GPGPUでうまくやる方法がわかりません。そこで、CPUに載っているコアのみで並列化します（Xeon Phiなら簡単なのでしょうか）。
一番外側の、0から(1<<25)-1まで変化する変数iのループをOpenMPで並列化します。変数iは上の図の緑の部分を5ビットで表現し、つなげたものです。マスに入りうる数は1から25までなので、5ビットというのはちょっと冗長ですが、とりあえずはこれでいいでしょう。
出力はバイナリ形式で、1つの解に25バイト使います（1つのマスに入る数を1バイトで表現する）。これもちょっと冗長ですが、テキスト形式よりは小さくて速いはずなので、とりあえずはこれでいいでしょう。
解は全部で2億7530万5224個、1つの解を25バイトで表現するので、（スレッド数と同数できる）出力ファイルのサイズの合計は275305224 x 25 = 6882630600バイト（約6.4GB）になります。
Core i7 4930K（12スレッド）で、出力先をRAMディスクにした結果は593秒（約10分）でした。
Core i7 4700MQ（8スレッド）で、出力先をSSDにした結果は1156秒（約20分）でした。
大森清美『魔方陣の世界』（日本評論社, 2013）（参考文献あり・索引あり）を見たら、a+b+c+d+e=65のとき、(26-a)+(26-b)+(26-c)+(26-d)+(26-e)=65なので、解の、すべての数xを(26-x)に置き換えたものもやはり解になるということが載っていました。この知識を使えばもっと速くなるでしょう（ここでは試しませんが）。
『魔方陣の世界』ではこの問題のためのコードも紹介されていますが、並列化されていないこともあって、解の数え上げで約5000秒、解の列挙で約3日と、あまりふるいません。
この問題は、パソコン・プログラミングに最適の題材である。(p.111)
今日では、5次方陣の総数検索問題は、パソコンに最適な題材となった。（p.124）
パソコンに最適とも、スパコンに不適とも思いませんが、興味深い題材であることは確かです。