Let's 'C#で数値計算' with Mol & Dsl

行列・ベクトルのデータ構造と配列演算

Molがサポートするベクトルには「密(dense)・疎(sparse)」、それと Complex・double を組み合わせた 4 種類、それに「密」な int 型、計 5 種類があります。
同様に行列には「密・疎」、「一般(general)・対称(symmetric)・三角(triangle)・帯(band)」、さらに double・Complex 等が組み合わさって様々な種類があります。

以下は、全てのベクトルと行列のクラス階層です。左、または左上、が各クラスの親(=スーパークラス)になります。


_Array _Vector                                     VectorDenseInt
               _VectorDouble                       VectorDenseDouble
                                                   VectorSparseDouble
               _VectorComplex                      VectorDenseComplex
                                                   VectorSparseComplex
       _Matrix _MatrixDouble _MatrixGeneralDouble  MatrixDenseGeneralDouble
                                                   MatrixDenseBandDouble
                                                   MatrixSparseGeneralDouble
                             _MatrixSquareDouble   MatrixDenseUpperTriangleDouble
                                                   MatrixDenseLowerTriangleDouble
                                                   MatrixDenseSymmetricDouble
                                                   MatrixSparseUpperTriangleDouble
                                                   MatrixSparseLowerTriangleDouble
                                                   MatrixSparseSymmetricDouble
               _MatriComplex _MatrixGeneralComplex MatrixDenseGeneralComplex
                                                   MatrixDenseBandComplex
                                                   MatrixSparseGeneralComplex
                             _MatrixSquareComplex  MatrixDenseUpperTriangleComplex
                                                   MatrixDenseLowerTriangleComplex
                                                   MatrixDenseSymmetricComplex
                                                   MatrixDenseHermite
                                                   MatrixSparseUpperTriangleComplex
                                                   MatrixSparseLowerTriangleComplex
                                                   MatrixSparseSymmetricComplex
                                                   MatrixSparseHermite


名前の先頭にアンダースコア '_' の付くクラスはスーパークラスになること(だけ)を想定されています。 独自のクラスをそれらから直接継承しないでください。

内部データ構造は、まず「密」な行列では

等の構造を持っています。詳細はMKLのドキュメント(Matrix Arguments)を参照するか、「mkl matrix packed storage format」等で検索すれば多くの情報が得られます。

同様に「疎」な行列は MKLのドキュメント(Sparse Matrix Storage Formats)を参照するか、「mkl matrix sparse storage format」等で検索すれば多くの情報が得られます。
「疎」な対称・三角行列(MatrixSparseSymmtricDouble等)は 「疎」な構造を持つとともに対角要素を含めた三角部分のみを保持します。

Mol で採用している内部データ構造は MKL の構造と完全に一致しています。
値は最終的に、全て一次元の配列(全てのベクトル・行列の親クラス _Array が管理します)に格納されます。

例えば A を密一般行列(サイズは M x N)、とすると _Array が管理する一次元配列の配列サイズは M x N となり、要素 A[i,j] は一次元配列[M*(j-1)+i] の位置に格納されます。
※疎な行列では最初の代入で、代入された行列要素のメモリー領域が確保されます(配列サイズは動的に増加します)。

行列のデータ構造が異なっていても、 Mol では宣言時の宣言文が異なるだけで、使用方法に違いはありません。 例えば、行列とベクトルの(A が行列、e、x、それに b はベクトル)演算 e = A*x - b は行列 A が何であれ、このように 書くことができるし、配列要素 A[i,j] にアクセスするには v = A[i,j] と書くだけです。

Mol ではこのようなデータ構造の違いをプログラマーから隠蔽、全ての行列やベクトル同士の演算(四則演算 + - * / )に対応しています。

さらに、構造保持のため、以下のように制約が強制されます。

メモリー管理

Mol は .Net 用のクラスライブラリー(Mol.Net.Dll)です。一方、MKL(その他の機能も含みます)は C や FORTRAN で作成された従来通りのDLL(Mol.C++.Dll)です。 C# 等でプログラミングする場合は、当然 Mol.Net.Dll が全ての管理を行います。 Mol.Net.Dll は連立方程式の解等、必要に応じて Mol.C++.Dll のエントリーを呼び出します。

さらに、ベクトルや行列の内部構造(配列要素値を格納する一次元配列等)は Mol.C++.Dll 側で確保・管理されます。
従って、例えば C# 等で配列要素 A[i,j] にアクセスする場合は .Net(Mol.Net.Dll)と従来のDLL(Mol.C++.Dll)間の データ交換操作が実行されますが、比較的コストの高い処理です。 以下のような頻繁に行列やベクトル要素にアクセスすることは避けて、簡単に C = A*B; 等と記述して、 計算は Mol.C++.Dll 内(.Net の外側)で行うようにしてください。

 int N = 1000;
 MatrixDenseGeneralDouble   A = new MatrixDenseGeneralDouble(N,N);
 MatrixDenseSymmetricDouble B = new MatrixDenseSymmetricDouble(N);
 ...........
 MatrixDenseGeneralDouble C = A*B; // 計算は Mol.C++.Dll 内(.Net の外側)で実行される。
 // 結果は同じでも、以下のような書き方は避けること。
 // for(int i=1;i<=N;++i)
 // {
 //    for(int j=0;j<=N;++j)
 //    {
 //        double s = 0;
 //        for(int k=1;k<=N;++k) s += A[i,k]*B[k,j]; // 掛け算は .Net 環境で実行される。
 //        C[i,j] = s;
 //    }
 // }
この様な四則演算子 + - * / 等を使う方法よりさらに高速・効率的な _Blas や _VML のメソッド群も 用意されています。

_Array クラス

_Array クラスには幾つかのプロパティやメソッドが用意されています。 _Array は行列やベクトルの要素を一次元の配列で管理しますが、行列の行や列の数等の情報は存在しません。 以下の表記で A は _Array クラスのインスタンス(ベクトルか行列)であるとします。 static なメソッドは _Array.SetConjugated(A) のように表記します。 更に複素数(Complex)の値を使用する場合、配列タイプが実数(または整数)の場合、虚数部分は使用されません。
呼出し形式説明備考
Complex A[int I]一次元配列 I 番目の要素にアクセスします。子クラスで再定義されるので ((_Array)A)[I] とキャストしてください。
int A.SizeLimit一次元配列の最大値。密な配列では以下の Size や Count は SizeLimit と同じ値になります。読み込みのみ
int A.Size現時点で確保されている一次元配列のサイズ(Size ≦ SizeLimit)。 読み込みのみ。
int A.Count現時点で値の設定されている要素の数(Count ≦ Size)。 読み込みのみ。
A.Negate()
void _Array.Negate(A)
A の全要素の符号を反転します。 
int A.Fill(Complex v)
int _Array.Fill(A,Complex v)
A の全要素に値 v を代入します。戻り値は A.Count
int _Array.CopyN(_Array target, int ix_t, _Array source, int ix_s, int N)source[ix_s]から N 個の要素を target[ix_t] の位置からコピーします。どちらかの配列サイズを超えない範囲で実行されます。戻り値はコピーした要素数
void A.SetConjugated()
void _Array.SetConjugated(A)
A の全要素を複素共役で置き換えます。 複素数配列のみ有効
double A.MaxAbs()
double _Array.MaxAbs(A)
Aの最大絶体値を返します。
double A.SumSquare()
double _Array.SumSquare(A)
Aの全要素の二乗和を返します。
double A.SumAbs()
double _Array.SumAbs(A)
Aの全要素の絶対値の和を返します。
ARRAY_TYPE ArrayType()
ARRAY_TYPE Complex ArrayType(A)
Aの配列タイプを返します。

以下は引数に行列やベクトルとスカラーを設定する _Array の(全て staticな)メソッド です。 c と a は _VectorDouble,_VectorComplex,_MatrixDouble,または _MatrixComplex。b は a と同じタイプ(double か Complex)の スカラー。 第一引数 c に null を指定すると a と同じサイズの c が新たに作成されてリターンされます。計算は c または a の配列サイズ範囲内で実行されます。 スカラー b は a の全ての(疎な行列・ベクトルの場合は非ゼロ)要素に作用します。結果は c に格納され、戻り値は c です。
呼出し形式説明
c = _Array.Add(c,a,b)c[i]=a[i]+b
c = _Array.Add(c,b,a)c[i]=b+a[i]
c = _Array.Sub(c,a,b)c[i]=a[i]-b
c = _Array.Sub(c,b,a)c[i]=b-a[i]
c = _Array.Mul(c,a,b)c[i]=a[i]*b
c = _Array.Mul(c,b,a)c[i]=b*a[i]
c = _Array.Div(c,a,b)c[i]=a[i]/b
c = _Array.Div(c,b,a)c[i]=b/a[i]
※ _Array クラスの static メソッドは最下位に位置するメソッドですが、要素全体を一気に操作するような場合に便利です。
※ 行列やベクトルにも同様の Add/Sub/Mul/Div が多数定義されています(従って、キャストが必要な場合があります。詳細はリファレンスを参照してください。)。
※ a がエルミート行列の場合、全要素に同じ値(特に虚数部がゼロでない複素数)を足したり引いたりするとエルミート性はなくなります。 ただし、エルミート行列は対角要素を含む上半分しか格納されていません。また、Mol の演算では対角要素の実数部分は無視(ゼロと扱われる)されるので見かけ上エルミート性は保存されます。 とはいえ、対角要素の虚数部値は、そのまま保存されるので _Blas 等 MKL 固有の処理部分に影響を与える可能性があるので注意してください。