「FORTRAN77数値計算プログラミング」のプログラムをCとPythonに移植してみる（その3）

前回の記事
「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回目、数式書くのが結構めんどくさかったです。