前回の記事
「FORTRAN77数値計算プログラミング」のプログラムをCとPythonに移植してみる（その2）
はじめに
その1、その2とFortran77数値計算プログラミングの1章のコードを確認してきました。
本では、1章にこのあと3つプログラムが載っているのですが、これ、全部誤差がどうなるかを確認しているものです。前回@cure_honey さんに教えていただいたとおり、16進数を基数とした浮動小数点数を使った計算結果が今のIEEE754方式の浮動小数点数での結果とは合わず、本での検証とはずれていくばかりなので、このあたりのコードはもう飛ばすことにします。
誤差について飛ばすと次は乱数のお話で、乱数生成ってCにもPythonにも関数があるのでそれも飛ばすと、お次は5章、「連立1次方程式とLU分解」です。
LU分解
本の5章、P.52〜の式を引きつつ、LU分解についてまとめます。
連立一次方程式の解き方
$n\times n$行列の行列$A$
A=\begin{pmatrix}
a_{11} & a_{12} & \cdots & a_{1n} \\
a_{21} & a_{22} & \cdots & a_{2n} \\
\vdots & \vdots & & \vdots \\
a_{n1} & a_{n2} & \cdots & a_{nn}
\end{pmatrix}\tag{1}
があって、$b$および$x$はそれぞれ次のようなベクトルのとき、
b = \begin{pmatrix}
b_1\\
b_2\\
\vdots\\
b_n
\end{pmatrix}\tag{2}
\\
x = \begin{pmatrix}
x_1\\
x_2\\
\vdots\\
x_n
\end{pmatrix}\tag{3}
\\
この$A$と$b$に対して
Ax = b\tag{4}
となる$x$を求める連立一次方程式を考えます。
LU分解とは
LU分解というのは、与えられた行列$A$を左下三角行列$L$と右上三角行列$U$との積で表すことです。つまり、
A=\begin{pmatrix}
■&0&0&0&0\\
■&■&0&0&0\\
■&■&■&0&0\\
■&■&■&■&0\\
■&■&■&■&■
\end{pmatrix}
\begin{pmatrix}
■&■&■&■&■\\
0&■&■&■&■\\
0&0&■&■&■\\
0&0&0&■&■\\
0&0&0&0&■
\end{pmatrix}
\equiv
LU
になるように■を定めることです。ただし実際は、
A=\begin{pmatrix}
1&0&0&0&0\\
■&1&0&0&0\\
■&■&1&0&0\\
■&■&■&1&0\\
■&■&■&■&1
\end{pmatrix}
\begin{pmatrix}
■&■&■&■&■\\
0&■&■&■&■\\
0&0&■&■&■\\
0&0&0&■&■\\
0&0&0&0&■
\end{pmatrix}
\equiv
LU\tag{5}
という形に限定します。
こうすることで、方程式$(4)$が
LUx=b\tag{6}
となって、これを解くためには、まず
Ly=b\tag{7}
の解$y$を求め、つぎにこの解$y$を右辺に持つ方程式
Ux=y\tag{8}
を解いて$x$を求めればよくなります。式が2つに増えてしまって、めんどくさそうですが、行列が三角行列になると、ずっと簡単に解くことが出来ます。
Fortranのコード
sdecom.f
PROGRAM SDECOM
* *
* Sample program of DECOMP and SOLVE *
* *
PARAMETER (NDIM = 20)
*
REAL A(NDIM,NDIM),B(NDIM),X(NDIM)
REAL W(NDIM),V(NDIM)
INTEGER IP(NDIM)
*
DO 10 NN = 1, 2
*
N = 10 * NN
*
CALL LUDATA (NDIM, N, A, B)
*
CALL ANCOMP (NDIM, N, A, ANORM)
*
WRITE (*,2000) ANORM
2000 FORMAT (/' ---- DECOMP and SOLVE ----'/
$ ' 1-norm of A = ',1PE13.6)
CALL DECOMP (NDIM, N, A, W, IP)
CALL SOLVE (NDIM, N, A, B, X, IP)
CALL ESTCND (NDIM, N, A, V, IP, W,
$ ANORM, COND)
*
WRITE (*,2001) (X(I), I=1,N)
2001 FORMAT (' solution X(I)'/(1X,5F9.5))
*
WRITE (*,2002) COND
2002 FORMAT (' condition number = ',1PE13.6)
*
10 CONTINUE
*
END
ludata.f
SUBROUTINE LUDATA (NDIM, N, A, B)
* *
* Sample data of DECOMP and SOLVE *
* *
REAL A(NDIM,N),B(N)
DO 510 I =1, N
DO 520 J = 1, N
A(I,J) = 0.0
520 CONTINUE
510 CONTINUE
*
A(1,1) = 4.0
A(1,2) = 1.0
DO 530 I = 2, N - 1
A(I,I-1) = 1.0
A(I,I) = 4.0
A(I,I+1) = 1.0
530 CONTINUE
A(N,N-1) = 1.0
A(N,N) = 4.0
DO 540 I = 1, N
B(I) = 0.0
DO 550 J = 1, N
B(I) = B(I) + A(I,J) * J
550 CONTINUE
540 CONTINUE
*
RETURN
END
ancomp.f
SUBROUTINE ANCOMP (NDIM, N, A, ANORM)
* *
* Compute 1-norm of A toestimate *
* condition number of A *
* *
REAL A(NDIM,N)
*
ANORM = 0.0
DO 10 K = 1, N
S = 0.0
DO 520 I = 1, N
S = S + ABS(A(I,K))
520 CONTINUE
IF (S .GT. ANORM) ANORM = S
10 CONTINUE
*
RETURN
END
decomp.f
SUBROUTINE DECOMP (NDIM, N, A, W, IP)
* *
* LU decomposition of N * N matrix A *
* *
REAL A(NDIM,N),W(N)
INTEGER IP(N)
*
DATA EPS / 1.0E-75 /
*
DO 510 K = 1, N
IP(K) = K
510 CONTINUE
*
DO 10 K = 1, N
L = K
AL = ABS(A(IP(L),K))
DO 520 I = K+1, N
IF (ABS(A(IP(I),K)) .GT. AL) THEN
L = I
AL = ABS(A(IP(L),K))
END IF
520 CONTINUE
*
IF (L .NE. K) THEN
*
* ---- row exchange ----
*
LV = IP(K)
IP(K) = IP(L)
IP(L) = LV
END IF
*
IF (ABS(A(IP(K),K)) .LE. EPS) GO TO 900
*
* ---- Gauss elimination ----
*
A(IP(K),K) = 1.0 / A(IP(K),K)
DO 30 I = K+1, N
A(IP(I),K) = A(IP(I),K) * A(IP(K),K)
DO 540 J = K+1, N
W(J) = A(IP(I),J) - A(IP(I),K) * A(IP(K),J)
540 CONTINUE
DO 550 J = K+1, N
A(IP(I),J) = W(J)
550 CONTINUE
30 CONTINUE
*
10 CONTINUE
*
RETURN
*
* ---- matrix singular ----
*
900 CONTINUE
WRITE (*,2001) K
2001 FORMAT (' (DECOMP) matrix singular ( at ',I4,
$ '-th pivot)')
N = - K
RETURN
*
END
solve.f
SUBROUTINE SOLVE (NDIM, N, A, B, X, IP)
* *
* Solve system of linear equations *
* A * X = B *
* *
REAL A(NDIM,N),B(N),X(N)
INTEGER IP(N)
*
REAL*8 T
*
* ---- forward substitution ----
*
DO 10 I = 1,N
T = B(IP(I))
DO 520 J = 1, I-1
T = T - A(IP(I),J) * DBLE(X(J))
520 CONTINUE
X(I) = T
10 CONTINUE
*
* ---- backward substitution ----
*
DO 30 I = N, 1, -1
T = X(I)
DO 540 J = I+1, N
T = T - A(IP(I),J) * DBLE(X(J))
540 CONTINUE
X(I) = T * A(IP(I),I)
30 CONTINUE
*
RETURN
END
estcnd.f
SUBROUTINE ESTCND(NDIM, N, A, V, IP, Y,
$ ANORM, COND)
* *
* Estimate condition number of A *
* *
REAL A(NDIM,N),V(N),Y(N)
INTEGER IP(N)
*
REAL*8 T
*
DO 10 K = 1,N
T = 0.0
DO 520 I = 1, k-1
T = T - A(IP(I),K) * DBLE(V(I))
520 CONTINUE
S = 1.0
IF (T .LT. 0.0) S = -1.0
V(K) = (S + T) * A(IP(K),K)
10 CONTINUE
*
DO 30 K = N, 1, -1
T = V(K)
DO 530 I = K+1, N
T = T - A(IP(I),K) * DBLE(Y(IP(I)))
530 CONTINUE
Y(IP(K)) = T
30 CONTINUE
*
YMAX = 0.0
DO 540 I = 1, N
IF (ABS(Y(I)) .GT. YMAX) YMAX = ABS(Y(I))
540 CONTINUE
*
COND = ANORM * YMAX
*
RETURN
END
コードの理解
sdecom
メインプログラム。例題のための配列の大きさを与えるNDIMをPARAMETER文で定義しています。
ludata
ludataというサブルーチンはまず、
N=10
のときに
A=\begin{pmatrix}
4.0 & 1.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\
1.0 & 4.0 & 1.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\
0.0 & 1.0 & 4.0 & 1.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\
0.0 & 0.0 & 1.0 & 4.0 & 1.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\
0.0 & 0.0 & 0.0 & 1.0 & 4.0 & 1.0 & 0.0 & 0.0 & 0.0 & 0.0 \\
0.0 & 0.0 & 0.0 & 0.0 & 1.0 & 4.0 & 1.0 & 0.0 & 0.0 & 0.0 \\
0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 1.0 & 4.0 & 1.0 & 0.0 & 0.0 \\
0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 1.0 & 4.0 & 1.0 & 0.0 \\
0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 4.0 & 1.0 \\
0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 1.0 & 4.0 \\
\end{pmatrix}
b=\begin{pmatrix}
6.00000 & 12.00000 & 18.00000 & 24.00000 & 30.00000 & 36.00000 & 42.00000 & 48.00000 & 54.00000 & 49.00000
\end{pmatrix}
という行列を作っています。
このとき、行列$b$は、
b_i = \sum_{j=1}^na_{ij}j
で求められています。
やりたいことは、この$A$と$b$に対して
Ax = b
となる$x$を求めることでした。
ancomp
次に、ancompというサブルーチンは、
||A||_1 = \max_{1\leq k\leq n} \sum_{i=1}^n|a_{ik}|\tag{9}
という行列のノルムに基づいて条件数を設定しています。これが$ANORM$という変数に設定されています。
ここで作業用変数$S$に倍精度実数型の宣言を行っていないのは、このノルム自身はそれほど制度を必要としないからだそう。
decomp
decompが、本題のLU分解を行うサブルーチン。
引数は、行列$A$に対応する2次元配列$A$とインデックスベクトル$p$に対応する1次元の整数型配列$IP$、行列の次元$N$です。
まず配列$IP$を、
IP(k) = k, \quad k=1,2,\cdots, n
で初期化してから、ピボットの部分選択を行いながらガウス消去法によってLU分解を実行します。
solve
LU分解した結果を使って方程式$(4)$を解くサブルーチンがsolveです。
estcnd
行列$A$の条件数の推定値を計算するサブルーチンがestcnd。ただし、行列$A$があらかじめdecompによってLU分解されていることが前提。
Cのコード
sdecom.h
#ifndef NDIM
#define NDIM 20
#endif
void ludata(int n, float a[NDIM][NDIM], float b[NDIM]);
float ancomp(int n, float a[NDIM][NDIM]);
int decomp(int n, float a[NDIM][NDIM], float w[NDIM], int ip[NDIM]);
void solve(int n, float a[NDIM][NDIM], float b[NDIM], float x[NDIM], int ip[NDIM]);
float estcnd(int n, float a[NDIM][NDIM], float v[NDIM], int ip[NDIM], float y[NDIM], float anorm);
sdecom.c
#include <stdio.h>
#include<stdlib.h>
#include "sdecom.h"
int main(void)
{
float a[NDIM][NDIM], b[NDIM], x[NDIM];
float w[NDIM], v[NDIM];
int ip[NDIM];
int nn;
int n;
float anorm, cond;
int i,j;
for(nn = 1; nn < 3; nn++){
n = 10 * nn;
ludata(n, a, b);
anorm = ancomp(n, a);
printf(" ---- DECOMP and SOLVE ----\n 1-norm of A = %13.6E\n"
, anorm);
n = decomp(n, a, w, ip);
solve(n, a, b, x, ip);
cond = estcnd(n, a, v, ip, w, anorm);
printf(" solution X(I)\n");
for(i=0; i<n; i++){
if(i != 0 && i % 5 == 0) {
printf("\n");
}
printf(" %.5f", x[i]);
}
printf("\n condition number = %13.6E\n", cond);
}
return 0;
}
ludata.c
#include <stdio.h>
#include "sdecom.h"
void ludata(int n, float a[NDIM][NDIM], float b[NDIM])
{
int i, j;
for(i=0; i<n; i++){
for(j=0; j<n; j++){
a[i][j] = 0.0f;
}
}
a[0][0] = 4.0f;
a[0][1] = 1.0f;
for(i=1; i<n-1; i++){
a[i][i-1] = 1.0f;
a[i][i] = 4.0f;
a[i][i+1] = 1.0f;
}
a[n-1][n-2] = 1.0f;
a[n-1][n-1] = 4.0f;
for(i=0; i<n; i++){
b[i] = 0.0f;
for(j=0; j<n; j++) {
b[i] += (a[i][j] * (j+1));
}
}
}
ancomp.c
#include <stdio.h>
#include <math.h>
#include "sdecom.h"
float ancomp(int n, float a[NDIM][NDIM])
{
int i,k;
float s;
float anorm = 0.0f;
for(k=0; k<n; k++){
s = 0.0f;
for(i=0; i<n; i++) {
s = s + fabsf(a[i][k]);
}
if (s > anorm){
anorm = s;
}
}
return anorm;
}
decomp.c
#include <stdio.h>
#include <math.h>
#include <float.h>
#include "sdecom.h"
int decomp(int n, float a[NDIM][NDIM], float w[NDIM], int ip[NDIM])
{
//double eps = 1.0e-75;
int i,j,k,l;
float al;
int lv;
for(k=0; k<n; k++){
ip[k] = k;
}
for(k=0; k<n; k++){
l = k;
al = fabsf(a[ip[l]][k]);
for(i=k+1; i<n; i++){
if(fabsf(a[ip[i]][k]) > al) {
l = i;
al = fabsf(a[ip[l]][k]);
}
}
if(l != k) {
lv = ip[k];
ip[k] = ip[l];
ip[l] = lv;
}
if((fabsf(a[ip[k]][k]) <= DBL_EPSILON)) {
printf(" (DECOMP) matrix singular ( at %4d-th pivot)\n", k);
n = -k;
break;
} else {
a[ip[k]][k] = 1.0f / a[ip[k]][k];
for(i=k+1; i<n; i++) {
a[ip[i]][k] = a[ip[i]][k] * a[ip[k]][k];
for(j=k+1; j<n; j++) {
w[j] = a[ip[i]][j] - a[ip[i]][k] * a[ip[k]][j];
}
for(j=k+1; j<n; j++) {
a[ip[i]][j] = w[j];
}
}
}
}
return n;
}
solve.c
#include <stdio.h>
#include "sdecom.h"
void solve(int n, float a[NDIM][NDIM], float b[NDIM], float x[NDIM], int ip[NDIM])
{
double t;
int i,j;
for(i=0; i<n; i++){
x[i] = 0.0f;
}
for(i=0; i<n; i++){
t = b[ip[i]];
for(j=0; j<i; j++){
t = t - a[ip[i]][j] * (double)x[j];
}
x[i] = t;
}
for(i=n-1; i>-1; i--){
t = x[i];
for(j=i+1; j <n; j++){
t = t - a[ip[i]][j] * (double)x[j];
}
x[i] = t * a[ip[i]][i];
}
}
estcnd.c
#include <stdio.h>
#include <math.h>
#include "sdecom.h"
float estcnd(int n, float a[NDIM][NDIM], float v[NDIM], int ip[NDIM], float y[NDIM], float anorm)
{
float cond;
double t;
int i,k;
float s;
float ymax;
for(k=0; k<n; k++){
t = 0.0;
for(i=0; i<k; i++){
t = t - a[ip[i]][k] * (double)v[i];
}
s = 1.0f;
if( t < 0.0f) {
s = -1.0f;
}
v[k] = (s + t) * a[ip[k]][k];
}
for(k=n-1; k>-1; k--){
t = v[k];
for(i=k+1; i<n; i++){
t = t - a[ip[i]][k] * (double)y[ip[i]];
}
y[ip[k]] = t;
}
ymax = 0.0f;
for(i=0; i<n; i++){
if(fabsf(y[i]) > ymax){
ymax = fabsf(y[i]);
}
}
cond = anorm * ymax;
return cond;
}
ただ愚直にFortranのコードを置き換え。
変数をいっそのことグローバルで定義しちゃおうかな、とか思ったものの、個人的に無理でやめました。
Pythonのコード
import scipy.linalg as linalg
import numpy as np
def ludata(n):
a = [[0.0] * n for i in range(n) ]
a[0][0] = 4.0
a[0][1] = 1.0
for i in range(1,n-1):
a[i][i-1] = 1.0
a[i][i] = 4.0
a[i][i+1] = 1.0
a[n-1][n-2] = 1.0
a[n-1][n-1] = 4.0
b = [0.0] * n
for i in range(n):
for j in range(n):
b[i] += (a[i][j] * (j+1))
return a,b
for nn in range(1,3):
n = 10 * nn
a,b = ludata(n)
LU = linalg.lu_factor(a)
x = linalg.lu_solve(LU, b)
print(x)
Pythonには、SciPyという強力な数値計算ライブラリがあり、当然のようにLU分解の関数もあるわけですね…。
そんなわけでゴリゴリ書くのはデータ作ってるところだけでした。
まとめ
探せばC言語にもなんか良さげなライブラリ公開されているんじゃないかなーとか思ったり。
とりあえず、だいぶ寝かせてしまった3回目、数式書くのが結構めんどくさかったです。