3

MatlabからC+GSLに移行していますが、行列Bを計算するための最も効率的な方法を知りたいです。

B[i][j] = exp(A[i][j])

ここで、iは[0、Ny]に、jは[0、Nx]にあります。

これは行列指数とは異なることに注意してください。

B = exp(A)

これは、GSL（linalg.h）の不安定な/サポートされていないコードで実現できます。

強引な解決策（「for」ループのカップル）を見つけましたが、それを行うためのよりスマートな方法はありますか？

編集

ドリューホールのソリューションポストからの結果

すべての結果は1024x1024for(for)ループからのものであり、各反復で2つのdouble値（複素数）が割り当てられます。時間は、100回の実行の平均時間です。

  • {Row、Column}を考慮した場合の結果-行列を格納するためのメジャーモード：
    • Row-Majorモード（ケース1）で、内側のループの行をループする場合、226.56ミリ秒。
    • 行メジャーモード（ケース2）で内部ループの列をループする場合は223.22ミリ秒。
    • GSLが提供する機能を使用する場合は224.60ミリ秒gsl_matrix_complex_set（ケース3）。

ケース1のソースコード：

for(i=0; i<Nx; i++)
{
    for(j=0; j<Ny; j++)
    {
        /* Operations to obtain c_value (including exponentiation) */
        matrix[2*(i*s_tda + j)] = GSL_REAL(c_value);
        matrix[2*(i*s_tda + j)+1] = GSL_IMAG(c_value);
    }
}

ケース2のソースコード：

for(i=0; i<Nx; i++)
{
    for(j=0; j<Ny; j++)
    {
        /* Operations to obtain c_value (including exponentiation) */
        matrix->data[2*(j*s_tda + i)] = GSL_REAL(c_value);
        matrix->data[2*(j*s_tda + i)+1] = GSL_IMAG(c_value);
    }
}

ケース3のソースコード：

for(i=0; i<Nx; i++)
{
    for(j=0; j<Ny; j++)
    {
        /* Operations to obtain c_value (including exponentiation) */
        gsl_matrix_complex_set(matrix, i, j, c_value);
    }
}
4

4 に答える 4

5

すべての要素を繰り返し処理し、各要素を呼び出すexp()か同等のものにすることを回避する方法はありません。しかし、反復する方法は速くも遅くもあります。

特に、目標はキャッシュミスを最小限に抑えることです。データが行優先または列優先のどちらの順序で保存されているかを確認し、内側のループがメモリに連続して保存されている要素を繰り返し処理し、外側のループが次の行に大きく進むようにループを配置してください（行メジャーの場合）または列（列メジャーの場合）。これは些細なことのように見えますが、パフォーマンスに大きな違いをもたらす可能性があります（マトリックスのサイズによって異なります）。

キャッシュを処理したら、次の目標はループのオーバーヘッドを取り除くことです。最初のステップ（マトリックスAPIがサポートしている場合）は、ネストされたループ（M＆N境界）から、基になるデータを反復処理する単一のループ（M N境界）に移行することです。これを行うには、基になるメモリブロックへの生のポインタ（つまり、doubleではなくdouble **）を取得する必要があります。

最後に、ループの展開をスローして（つまり、ループの反復ごとに8または16の要素を実行し）、ループのオーバーヘッドをさらに削減します。これは、おそらく可能な限り高速です。残りの要素をクリーンアップするには、フォールスルーを含む最後のswitchステートメントが必要になる可能性があります（配列サイズ％ブロックサイズ！= 0の場合）。

于 2010-07-24T02:48:27.980 に答える
3

いいえ、聞いたことのない奇妙な数学的癖がない限り、2つのforループで要素をループする必要があります。

于 2010-07-23T21:38:58.460 に答える
2

数値の配列に適用したいだけの場合exp、実際には近道はありません。あなたはそれを（Nx * Ny）回呼ばなければなりません。行列要素の一部が0のように単純な場合、または要素が繰り返されている場合は、メモ化が役立つ可能性があります。

ただし、本当に必要なのが行列指数関数（これは非常に便利です）である場合、依存するアルゴリズムはDGPADMです。これはFortranにありますが、f2cを使用してCに変換できます。これがその論文です。

于 2010-07-24T02:04:45.763 に答える
0

ループの内容が表示されていないため、c_valueを計算するビットは、コードのパフォーマンスがメモリ帯域幅によって制限されているのか、CPUによって制限されているのかわかりません。確実に知る唯一の方法は、プロファイラーと洗練されたプロファイラーを使用することです。メモリレイテンシ、つまりCPUがRAMからデータが到着するのを待ってアイドル状態になっている時間を測定できる必要があります。

メモリ帯域幅によって制限されている場合、メモリに順番にアクセスすると、できることは多くありません。CPUとメモリは、データが順番にフェッチされるときに最適に機能します。データをRAMからキャッシュにフェッチする必要がある可能性が高いため、ランダムアクセスはスループットに影響を与えます。あなたはいつでもより速いRAMを手に入れようとすることができます。

CPUによって制限されている場合は、さらにいくつかのオプションを利用できます。浮動小数点コードを手動でコーディングする場合と同様に、SIMDを使用することも1つのオプションです（C / C ++コンパイラは多くの理由でFPUコードに優れていません）。これが私であり、内側のループのコードで許可されている場合、配列への2つのポインターがあります。1つは先頭に、もう1つは配列の4/5にあります。反復ごとに、最初のポインターを使用してSIMD操作が実行され、2番目のポインターを使用してスカラーFPU操作が実行されるため、ループの各反復は5つの値を実行します。次に、レイテンシコストを軽減するために、SIMD命令とFPU命令をインターリーブします。（少なくともPentiumでは）MMUは最大4つのデータストリームを同時にストリーミングできるため、これはキャッシュに影響を与えることはありません（つまり、

于 2010-07-25T23:12:54.723 に答える