インテル® MKL クックブック

定常非線形熱伝導方程式の近似解を求める

目的

熱方程式の境界値問題の解と、その解に依存する熱係数を得る。

ソリューション

固定小数点の反復アプローチ [ Amos10] を使用し、インテル® MKL PARDISO によって各外部反復の線形問題を解きます。

  1. CSR (Compressed Sparse Rows: 圧縮スパース行) 形式で行列構造を設定します。

  2. 残差のノルムが許容差未満になるまで固定小数点の反復を行います。

    1. pardiso ルーチンを使用して現在の反復の線形化されたシステムを解きます。

    2. dcopy ルーチンを使用して、システムの解を主方程式の次の近似に設定します。

    3. 新しい近似に基づいて、行列の新しい要素を計算します。

    4. mkl_cspblas_dcsrgemv ルーチンを使用して現在の解の残差を計算します。

    5. dnrm2 ルーチンを使用して残差のノルムを計算し、許容差と比較します。

  3. ソルバーの内部メモリーを解放します。

ソースコード: サンプル (http://software.intel.com/en-us/mkl_cookbook_samples) の sparse フォルダーを参照してください。

インテル® MKL PARDISO、スパース BLAS、BLAS を使用して近似解を求める

      CONSTRUCT_MATRIX_STRUCTURE (nx, ny, nz, &ia, &ja, &a, &error);
      CONSTRUCT_MATRIX_VALUES (nx, ny, nz, us, ia, ja, a, &error);
      DO WHILE res > tolerance
        phase = 13;
        PARDISO (pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, &idum, 
          &nrhs, iparm, &msglvl, f, u, &error );      
        DCOPY (&n, u, &one, u_next, &one);
        CONSTRUCT_MATRIX_VALUES (nx, ny, nz, u_next, ia, ja, a, &error);
        MKL_CSPBLAS_DCSRGEMV (uplo, &n, a, ia, ja, u, temp );
        DAXPY (&n, &minus_one, f, &one, temp, &one);
        res = DNRM2 (&n, temp, &one);
      END DO
      phase = -1;
      PARDISO ( pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, &idum,
        &nrhs, iparm, &msglvl, f, u, &error );

使用するルーチン

タスク

ルーチン

説明

現在の反復の線形化されたシステムを解き、ソルバーの内部メモリーを解放する

PARDISO

複数の右辺を備えた 1 セットのスパース線形方程式の解を計算します。

見つかった解を主方程式の次の近似として設定する

DCOPY

ベクトルを別のベクトルにコピーします。

現在の非線形反復の残差を計算する

MKL_CSPBLAS_DCSRGEMV

ゼロベースでインデックスされ、CSR 形式で格納されたスパース一般行列の行列-ベクトル積を計算します。

DAXPY

ベクトル-スカラー積を計算して結果をベクトルに追加します。

残差のノルムを計算し、停止条件と比較する

DNRM2

ベクトルのユークリッド・ノルムを計算します。

説明

定常非線形熱伝導方程式は、非線形偏微分方程式の境界値問題として説明できます。

ここで、領域 D は立方体であると仮定します。D = (0, 1)3 および v(x, y, z) は温度の未知関数です。

デモ目的のため、問題は解の熱係数の線形従属性に制限されています。

数の解を得るため、領域 D でグリッドステップ h の等距離のグリッドが選択され、偏微分方程式は差分を使用して近似されます。 このプロシージャーから [Smith86] 非線形代数方程式のシステムが得られます。

各方程式は、7 つのグリッドポイントで、未知グリッド関数 u の値とそれぞれの右辺の値を関連付けます。 方程式の左辺は、解に依存するグリッド関数値と係数の線形の組み合わせとして表せます。これらの係数から構成される行列を利用すると、方程式をベクトル-行列形式で書き直すことができます。

係数行列 A は疎である (各行に非ゼロ要素が 7 つしかない) ため、この反復アルゴリズムで解くために、行列を CSR 形式の配列に格納して (『インテル® MKL リファレンス・マニュアル』の「スパース行列格納形式」を参照)、PARDISO* ソルバーを使用することは適切です。

  1. u を初期値 u0 に設定します。

  2. 残差 r = A(u)u - g を計算します。

  3. ||r|| < 許容差の場合:

    1. 方程式 A(u)w = g を w について解きます。

    2. u = w を設定します。

    3. 残差 r = A(u)u - g を計算します。

親トピック: インテル® MKL クックブック