行列データとベンチマーク

まず、比較的小規模で予め逆行列や固有値が記載された例題集として
　 Collection of Matrices for Testing Computational Algorithms(Robert T. Gregory , David L. Karney)
が便利です。

大規模な行列なら、テキストファイルやプログラム形式でダウンロード可能な
　The Matrix Market
があります。

※以下、計算は Intel CORE i7-2600S 2.8GHz 8GB メモリ マシンでの結果です。

行列のダウンロードと読み込み

そのなかで、大規模な例として E40R5000: Driven cavity （一般行列）と、ZENIOS: Air-traffic Control Model （対称行列）を用いて連立方程式と固有値の計算を実行してみます（データファイル書式は共にMatrixMarket format形式の e40r5000.mtx.gz と zenios.mtx.gz をダウンロード）。

データは圧縮されているので解凍する必要があります。
解凍ソフトが手元にない場合は 　The GZip home page から ZIP 形式で圧縮された gzip124xN.zip か自己解凍形式の gzip124xN.exeをダウンロードして gzip.exe をコマンドプロンプト上で実行できるようにしてください。

例えば、コマンドプロンプト画面から zenios.mtx.gz のあるフォルダ(D:\Download\gzip> とします)に移動してから

　D:\Download\gzip>gzip -d zenios.mtx.gz

%%MatrixMarket matrix coordinate real symmetric
2873 2873 15032
1 1  0.0000000000000e+00
2 2  0.0000000000000e+00
10 2  2.1347330876700e-01
14 2  2.0394863343300e-01
18 2  6.4994637116100e-02
21 2  1.2923523430100e-01
25 2  4.2872933495400e-01
...............

   string[] ReadText(StreamReader sr)
   {
    Retry:
       string st = sr.ReadLine();
       if (st.StartsWith("%")) goto Retry;
       if (st.StartsWith(";")) goto Retry;
       string[] sts = st.Split(' ');
       int j = 0;
       for (int i = 0; i < sts.Length; ++i)
       {
           if (sts[i].Trim() == "") continue;
           sts[j++] = sts[i];
       }
       return sts;
   }

   //
   // データをファイル（path）から読み込み T で指定された行列を作成後、要素データをセットしてリターンします。
   _Matrix MatrixFromFile<T>(string path)
   {
       StreamReader sr = new StreamReader(path);
       string[] st = ReadText(sr);
       int N = int.Parse(st[0]);
       int M = int.Parse(st[1]);
       int L = int.Parse(st[2]);
       Type t = typeof(T);
       _Matrix A = null;
       if (t.AssemblyQualifiedName.Contains("General"))
       {
          // 一般行列の作成
          A = (_Matrix)Activator.CreateInstance(t, N, M);
       }
       else
       {
          // 正方行列の作成
          A = (_Matrix)Activator.CreateInstance(t, N);
       }
       for (int i = 1; i <= L; ++i)
       {
          st = ReadText(sr);
          int I = int.Parse(st[0]);
          int J = int.Parse(st[1]);
          double re = double.Parse(st[2]);
          if (A is _MatrixDouble)
          {
              if(re!=0.0) ((_MatrixDouble)A)[I, J] = re;
          }
          else
          {
              double im = double.Parse(st[3]);
              if(re!=0.0 || im!=0.0) ((_MatrixComplex)A)[I, J] = new Complex(re,im);
          }
       }
       sr.Close();
　     return A;
  }

各データは共にサイズは大きいものの疎な行列です。このようなデータを密な行列に格納することは無駄どころか、計算不可能になり得ます。 例えば、e40r5000.mtx の行数（＝列数）は 17281 ですので全要素数は 298632961（=17281*17281)個となります。各要素は double の１６バイトですから必要バイト数は 4778127376 バイトと途方もない大きさになります。

連立方程式の解:e40r5000.mtx

1. A を因子（LU）分解(LuSolver lu = LuSolver.Create(A);)して、新規に作成した LuSolver オブジェクトに結果を格納する。
2. 作成された LuSolver オブジェクトの Solve メソッドに b を与えて x を計算(x = lu.Solve(b);)する。

※一度因子分解すれば、異なる b を与えてx=lu.Solve(b);を何度も呼び出すことができます。b に行列を指定することもできます。

特に疎な行列を指定した場合、スイスのバーゼル大学（Universität Basel）で開発された有名な pardiso(Parallel Direct and Iterative Solvers) が 使用されます。 MKL ではインテルCPU用に pardiso を最適化したものが提供されています（pardiso mkl 等で検索すれば色々な情報を得ることができます）。 読み込みと実行プログラムは以下のようになります。

  DateTime dt = DateTime.Now;
  MatrixSparseGeneralDouble A = (MatrixSparseGeneralDouble)MatrixFromFile<MatrixSparseGeneralDouble>("D:\\Download\\gzip\\e40r5000.mtx");
  TimeSpan ts = DateTime.Now - dt;
  T.WL("読み込み時間=" + ts + " 行列サイズ=" + A.RowCount + " 非ゼロ要素数="+A.Count);
  VectorDenseDouble b = new VectorDenseDouble(A.ColumnCount);
  // 定数ベクトルの全要素を 1.0 にする.
  _Array.Fill(b,new Complex(1.0,0.0));
  dt = DateTime.Now;
  LuSolver lu = LuSolver.Create(A);
  ts = DateTime.Now - dt;
  T.WL("因子分解時間=" + ts);
  dt = DateTime.Now;
  VectorDenseDouble x = lu.Solve(b);
  ts = DateTime.Now - dt;
  T.WL("解法時間=" + ts);
  // 残差 d = Ax-b の計算
  VectorDenseDouble d = (VectorDenseDouble)b.Clone();
  _Blas.LeMv(d, A, x, _Mol.MATRIX_OPERATION.AS_IS, 1.0, -1.0);
  T.WL("Norm=" + d.NormE());
  T.WL("最大残差=" + _Array.MaxAbs(d));

　読み込み時間=00:00:26.2704461 行列サイズ=17281 非ゼロ要素数=558361
　因子分解時間=00:00:00.2808005
　解法時間=00:00:00.0156001
　Norm=1.71579434333068E-08
　最大残差=2.41149766821991E-09

固有値問題の解:zenios.mtx

  DateTime dt = DateTime.Now;
  // MatrixSparseSymmetricDouble A = (MatrixSparseSymmetricDouble)MatrixFromFile("D:\\Download\\gzip\\zenios.mtx");
  MatrixDenseSymmetricDouble A = (MatrixDenseSymmetricDouble)MatrixFromFile("D:\\Download\\gzip\\zenios.mtx");
  TimeSpan ts = DateTime.Now - dt;
  T.WL("読み込み時間=" + ts + " 行列サイズ=" + A.SquareSize + " 非ゼロ要素数=" + A.Count);
  VectorDenseDouble        eigen_values  = new VectorDenseDouble(A.SquareSize); // 固有値
  MatrixDenseGeneralDouble eigen_vectors = new MatrixDenseGeneralDouble(A.SquareSize, A.SquareSize); // 固有ベクトル
  dt = DateTime.Now;
  // 範囲 [-1.0,+1.0] の固有値を求める。 
  // EgvSolver egv = EgvSolver.Solve(eigen_values, eigen_vectors, (MatrixSparseSymmetricDouble)A.Clone(), -1.0, 1.0);
  EgvSolver egv = EgvSolver.Solve(eigen_values, eigen_vectors, (MatrixDenseSymmetricDouble)A.Clone(), -1.0, 1.0,1.0e-6);
  ts = DateTime.Now - dt;
  T.WL("計算時間=" + ts);
  int n = egv.Count;
  T.WL("計算された固有値数=" + n);
  double s = 0;
  for (int j = 1; j <= n; ++j)
  {
     VectorDenseDouble x = eigen_vectors.Columns[j]; // 固有ベクトル
     // Ax = λx の確認（λは固有値==eigen_values[j]、x は固有ベクトル。）
     _Blas.LeMv(x, A, (VectorDenseDouble)x.Clone(), 1.0, -eigen_values[j]);
     s += x.NormE();
  }
  T.WL("残差=" + s);

密な行列の結果は以下の通りです。密な対称行列ですので対角要素を含む 上三角部分の全要素数 4128501（=2873(2873+1)/2）が「非ゼロ要素数」（確保された領域数）になります。

　読み込み時間=00:00:00.0468001 行列サイズ=2873 非ゼロ要素数=4128501
　計算時間=00:00:05.6784100
　計算された固有値数=2857
　残差=0.000627841236062497

　読み込み時間=00:00:00.0468001 行列サイズ=2873 非ゼロ要素数=3530
　計算時間=00:01:17.2513357
　計算された固有値数=2857
　残差=3.47111366778501E-13