Previous article Porting the "FORTRAN77 Numerical Computing Programming" program to C and Python (Part 2)
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".
The LU decomposition is summarized by drawing the formulas from Chapter 5, P.52 of the book.
$ 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 $.
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.
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 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.
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.
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.
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