帯行列の計算

帯行列の計算 2013年9月19日後　保範

帯行列の三角分解 0 0 0 0 = ・ 0 0 0 0 A UT U 規則的疎行列半帯の中が非ゼロになる

スカイライン行列との関係 帯幅の位置まで広げる有限要素法不規則疎行列不規則疎行列帯行列

対称正定値行列の発生 有限要素法で離散化対称正定置帯行列軸交換なし、上半分で分解

対称正定値帯行列 消去対象この部分だけ計算

対称帯行列の消去 Wにコピー j k A[i][j-i] -= T*W[j] i T = A[i][i-k]/A[k][0] 消去計算 0 <= k <n k< i <=k+m, i<= j <i+m

対称帯行列の消去計算 for (k=0; k<n; k++) { bend = Math.min(k+m, n-1); for (j=k+1; j<=bend; j++) { W[j] = A[k][j-k]; } for (i=k+1; i<=bend; i++) { T = A[k][i-k] / A[k][0]; for (j=i; j<=bend; j++) { A[i][j-i] -= T*W[j]; } } }

対称帯行列の前進後退代入 for (k=0; k<n; k++) { BK = B[k][0]; bend = Math.min(k+m,n-1) for (i=k+1; i<=bend; i++) { B[i] -= BK*A[k][i-k]; } } for (k=n-1; k>=0; k--) { S = B[k]; bend = Math.min(k+m,n-1); for (j=k+1; j<=bend; j++) { S -= B[j]*A[k][j-k]; } B[k] = S / A[k][0]; }

非対称行列の発生 有限要素法で離散化非対称帯行列軸交換有り、半幅拡大して分解

非対称帯行列 軸交換によりこの部分が増える消去対象

非対称帯行列の軸交換 A[k][j-k+m] k 交換 A[ip][j-IPK+m] IPK コピー W[j] |A[ip][k-IPK+m]| 　が最大値

非対称帯行列の消去 A[i][j-i+m] -=T*W[j] j k i 消去計算 T=A[i][k-i+m]/A[k][m]

非対称帯行列の消去計算１（全体） for (k=0; k<n; k++) { 軸交換処理(最大値探査と軸交換) if (最大値 > ε) { 消去計算 } else { 特異行列と判定 } }

消去計算２（軸交換） AMAX = Math.abs(A[k][m]); IPK = Math.abs(A[k][m]); for (i=k+1; i<= Math.min(k+m,n-1); i++) { AIK = Math.abs(A[i][k-i+m]); if ( AIK > AMAX) { IPk = i; AMAX = AIK; } } IP[k] = IPK; for (j=k; j <= Math,min(k+2*m,n-1); j++) { W[j] = A[IPK][j-IPK+m]; A[IPK][j-IPK+m] = A[k][j-k+m]; A[k][j-k+m] = W[j]; }

消去計算３（消去計算） for (i=k+1; i <= Math.min(k+m,n-1); i++) { T = A[i][k-i+m] / A[k][m]; A[i][k-i+m] = T; for (j=k+1; j <= Math.min(k+2*m,n-1; j++) { A[i][j-i+m] -= T*W[j]; } }

非対称帯行列の前進代入 for (k=0; k<n; k++) { IPK = IP[k]; if ( IPK != k) { WK = B[IPK]; B[IPK]= B[k]; B[k]=W; } for (i=k+1; Math.min(k+m,n-1); i++) { B[i] -= B[k]*A[i][k-i+m]; } }

非対称帯行列の後退代入 for (k=n-1; k >= 0; k--) { S = B[k]; for (j=k+1; j<= Math.min(k+2*m,n-1); j++) { S -= B[j]*A[k][j-k+m]; } B[k] = S / A[k][m]; }

帯行列の計算

帯行列の計算

Presentation Transcript