帯行列の計算
2013年9月19日
後 保範
1
帯行列の三角分解
0
0
0
・
=
0
0
0
0
0
A
規則的疎行列
UT
U
半帯の中が非ゼロになる
2
スカイライン行列との関係
帯幅の位置
まで広げる
有限要
素法
不規則
疎行列
不規則
疎行列
帯行列
3
対称正定値行列の発生
 div(k  grad(u))  f
有限要素法で離散化
対称正定置帯行列
軸交換なし、上半分で分解
4
対称正定値帯行列
消去対象
この部分
だけ計算
5
対称帯行列の消去
j
k
Wにコピー
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
6
対称帯行列の消去計算
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]; }
}
}
7
対称帯行列の前進後退代入
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];
}
8
非対称行列の発生
 div(k  grad(u))  v  grad(u)  f
有限要素法で離散化
非対称帯行列
軸交換有り、半幅拡大して分解
9
非対称帯行列
軸交換により
この部分が
増える
消去対象
10
非対称帯行列の軸交換
k
IPK
|A[ip][k-IPK+m]|
が最大値
A[k][j-k+m]
交換
A[ip][j-IPK+m]
コピー
W[j]
11
非対称帯行列の消去
A[i][j-i+m] -=T*W[j]
k
j
i
T=A[i][k-i+m]/A[k][m]
消去計算
12
非対称帯行列の消去計算１（全体）
for (k=0; k<n; k++)
{ 軸交換処理(最大値探査と軸交換)
if (最大値 > ε)
{ 消去計算 }
else
{ 特異行列と判定 }
}
13
消去計算２（軸交換）
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]; }
14
消去計算３（消去計算）
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]; }
}
15
非対称帯行列の前進代入
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]; }
}
16
非対称帯行列の後退代入
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];
}
17

NO.280411 潜水作業船・ヤンマーDE43C 3級