Try porting the "FORTRAN77 Numerical Computing Programming" program to C and Python (Part 3)

Previous article Porting the "FORTRAN77 Numerical Computing Programming" program to C and Python (Part 2)

Introduction

I have checked the code of Chapter 1 of Part 1, Part 2 and Fortran77 Numerical Computing Programming. In the book, there are three programs after this in Chapter 1, but all of them are checking what the error will be. As @cure_honey told me last time, the calculation result using floating point numbers based on hexadecimal numbers does not match the result with the current IEEE754 floating point number, and it deviates from the verification in the book. I'm just skipping the code around here. If you skip the error, next is the story of random numbers. Since there are functions in both C and Python for random number generation, if you skip that as well, the next is Chapter 5, "Simultaneous linear equations and LU decomposition".

LU decomposition

The LU decomposition is summarized by drawing the formulas from Chapter 5, P.52 of the book.

How to solve simultaneous linear equations

$ n \ times n $ matrix matrix $ 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}

And when $ b $ and $ x $ are the following vectors, respectively,

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}
\\

For this $ A $ and $ b $

Ax = b\tag{4}

Consider a system of linear equations for $ x $.

What is LU decomposition?

LU decomposition means that the given matrix $ A $ is expressed as the product of the lower left triangular matrix $ L $ and the upper right triangular matrix $ U $. In other words

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

It is to determine ■ so that it becomes. However, in reality

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}

It is limited to the form. By doing this, the equation $ (4) $

LUx=b\tag{6}

To solve this, first

Ly=b\tag{7}

Find the solution $ y $ of, then the equation with this solution $ y $ on the right side

Ux=y\tag{8}

All you have to do is solve to find $ x $. It seems to be annoying because the number of equations has increased to two, but when the matrix becomes a triangular matrix, it is much easier to solve.

Fortran code

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 

Understanding the code

sdecom Main program. The PARAMETER statement defines NDIM, which gives the size of the array for the example.

ludata First of all, the subroutine called ludata

N=10

When

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}

I am making a procession called. At this time, the matrix $ b $ is

b_i = \sum_{j=1}^na_{ij}j

Is sought after. What I want to do is for this $ A $ and $ b $

Ax = b

Was to find $ x $.

ancomp Next, the subroutine called ancomp

||A||_1 = \max_{1\leq k\leq n} \sum_{i=1}^n|a_{ik}|\tag{9}

The conditional number is set based on the norm of the matrix. This is set in a variable called $ ANORM $. The reason why the double-precision real type is not declared for the working variable $ S $ here is that this norm itself does not require much institution.

decomp A subroutine in which decomp performs the LU decomposition of the main subject. The arguments are a two-dimensional array $ A $ corresponding to the matrix $ A $, a one-dimensional integer array $ IP $ corresponding to the index vector $ p $, and a matrix dimension $ N $. First, the array $ IP $,

IP(k) = k, \quad k=1,2,\cdots, n

After initializing with, LU decomposition is performed by Gaussian elimination while performing partial selection of the pivot.

solve The subroutine that solves the equation $ (4) $ using the result of LU decomposition is solve.

estcnd The subroutine that calculates the estimated number of conditions in the matrix $ A $ is estcnd. However, it is assumed that the matrix $ A $ has been LU-decomposed by decomp in advance.

C code

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;
}

Just honestly replace the Fortran code. I wondered if I could define variables globally, but I personally couldn't do it.

Python code

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 has a powerful numerical calculation library called SciPy, and of course there is also an LU decomposition function ... That's why I only wrote data when I was making data.

Summary

If you look for it, I wonder if there is a nice library published in C language as well. For the time being, it was quite annoying to write the formula for the third time, which I had put to sleep for a long time.

Recommended Posts

Try porting the "FORTRAN77 Numerical Computing Programming" program to C and Python (Part 1)
Try porting the "FORTRAN77 Numerical Computing Programming" program to C and Python (Part 3)
Try porting the "FORTRAN77 Numerical Computing Programming" program to C and Python (Part 2)
Try to solve the programming challenge book with python3
[Introduction to Python3 Day 1] Programming and Python
[Introduction to Python] I compared the naming conventions of C # and Python.
Try to write a program that abuses the program and sends 100 emails
Try to implement and understand the segment tree step by step (python)
Put Cabocha 0.68 on Windows and try to analyze the dependency with Python
C language to see and remember Part 2 Call C language from Python (argument) string
C language to see and remember Part 1 Call C language from Python (hello world)
Porting and modifying doublet-solver from python2 to python3.
Try to write python code to generate go code --Try porting JSON-to-Go and so on
Python amateurs try to summarize the list ①
C language to see and remember Part 4 Call C language from Python (argument) double
C language to see and remember Part 5 Call C language from Python (argument) Array
C language to see and remember Part 3 Call C language from Python (argument) c = a + b
How to generate permutations in Python and C ++
Try to solve the Python class inheritance problem
Try using the Python web framework Tornado Part 1
Try using the Python web framework Tornado Part 2
"Cython" tutorial to make Python explosive: Pass the C ++ class object to the class object on the Python side. Part 2
The story of porting code from C to Go and getting hooked (and to the language spec)
Think about how to program Python on the iPad
[Python] Try to read the cool answer to the FizzBuzz problem
Try to visualize the room with Raspberry Pi, part 1
Try to solve the internship assignment problem with Python
Try to operate DB with Python and visualize with d3
Try embedding Python in a C ++ program with pybind11
Tips and precautions when porting MATLAB programs to Python
Specifies the function to execute when the python program ends
I felt that I ported the Python code to C ++ 98.
How to enjoy Python on Android !! Programming on the go !!
Go language to see and remember Part 7 C language in GO language
[C / C ++] Pass the value calculated in C / C ++ to a python function to execute the process, and use that value in C / C ++.
"Cython" tutorial to make Python explosive: Pass a C ++ class object to a class object on the Python side. Part ①
How to start the PC at a fixed time every morning and execute the python program
Try to import to the database by manipulating ShapeFile of national land numerical information with Python
A story about porting the code of "Try and understand how Linux works" to Rust
[Python] A program to find the number of apples and oranges that can be harvested
[Introduction to cx_Oracle] (Part 6) DB and Python data type mapping
Try to make it using GUI and PyQt in Python
Try hitting the Twitter API quickly and easily with Python
Try to get the function list of Python> os package
I wanted to solve the Panasonic Programming Contest 2020 with Python
Try connecting to Supervisord via XMLRPC to start / stop the program
Measure the execution result of the program in C ++, Java, Python.
Write a python program to find the editing distance [python] [Levenshtein distance]
Consider If Programming Was An Anime in Python and C
[Part.2] Crawling with Python! Click the web page to move!
I tried to illustrate the time and time in C language
Try to bring up a subwindow with PyQt5 and Python
Build a Python environment and transfer data to the server
I tried programming the chi-square test in Python and Java.
Try to automate the operation of network devices with Python
Solving with Ruby and Python AtCoder ABC011 C Dynamic programming
I want to know the features of Python and pip
I tried to enumerate the differences between java and python
Object-oriented in C: Refactored "○ ✕ game" and ported it to Python
Program to determine leap year from the Christian era [Python]
Just try to receive a webhook in ngrok and python