#!/usr/bin/env python
# coding: utf-8

# # Table of Contents
#  <p><div class="lev1 toc-item"><a href="#3x3行列のLU分解と解" data-toc-modified-id="3x3行列のLU分解と解-1"><span class="toc-item-num">1&nbsp;&nbsp;</span>3x3行列のLU分解と解</a></div><div class="lev1 toc-item"><a href="#pivot付きのLU分解" data-toc-modified-id="pivot付きのLU分解-2"><span class="toc-item-num">2&nbsp;&nbsp;</span>pivot付きのLU分解</a></div><div class="lev1 toc-item"><a href="#Jacobi-vs-Gauss-Seidel" data-toc-modified-id="Jacobi-vs-Gauss-Seidel-3"><span class="toc-item-num">3&nbsp;&nbsp;</span>Jacobi vs Gauss-Seidel</a></div>

# # 3x3行列のLU分解と解

# In[1]:


from pprint import pprint
import scipy.linalg as linalg   # SciPy Linear Algebra Library
import numpy as np

# A = np.array([[1,1,-2,-4],[1,-2,1,5],[2,-2,-1,2]])
A = np.array([[1.0,4,3],[1,-2,1],[2,-2,-1]])
b = np.array([[11.0],[11],[11]])
pprint(A)
pprint(b)


# In[2]:


inv_A = linalg.inv(A)
print(inv_A)
np.dot(inv_A,b)


# In[3]:


e_A = np.hstack((A, b))
pprint(e_A)
pprint(np.column_stack((A,b)))


# In[4]:


P, L, U = linalg.lu(e_A)

print("P:")
pprint(P)

print("L:")
pprint(L)

print("U:")
pprint(U)


# In[5]:


np.dot(L,U)


# In[1]:


from pprint import pprint
import scipy.linalg as linalg   # SciPy Linear Algebra Library
import numpy as np


A = np.array([[1.0,4,3],[1,-2,1],[2,-2,-1]])
A0 = np.array([[1.0,4,3],[1,-2,1],[2,-2,-1]])
b = np.array([[11.0],[11],[11]])


n = 3
L = np.identity(n)
T = []
for i in range(n):        #i行目
    T.append(np.identity(n))
    for j in range(i+1, n):
        am = A[j,i]/A[i,i]    #i行の要素を使って，i+1行目の先頭を消す係数を求める
        T[i][j,i]=-am         #i番目の消去行列に要素を入れる
        L[j,i]=am             #LTMの要素
        for k in range(n):
            A[j,k] -= am*A[i,k] #もとの行列をUTMにしていく
        b[j] -= b[i]*am   #数値ベクトルも操作
pprint(A)
pprint(b)


# In[2]:


np.dot(L,A)


# In[3]:


pprint(L)


# # pivot付きのLU分解

# In[5]:


from pprint import pprint
import scipy.linalg as linalg   # SciPy Linear Algebra Library
import numpy as np

A=np.array([[3,2,2,1],[3,2,3,1],[1,-2,-3,1],[5,3,-2,5]])
b=np.array([-6,2,-9,2])
ex_A = np.column_stack((A,b))
P,L,U = linalg.lu(ex_A)
pprint(P)
pprint(L)
pprint(U)

pprint(np.dot(P,A))
pprint(np.dot(L,U))


# # Jacobi vs Gauss-Seidel
# 
# ```
# [-1.2       0.666667  0.777778  0.6     ]
# [-1.608889  0.607407  0.562963  0.751111]
# [-1.584296  0.764938  0.54749   0.782519]
# [-1.618989  0.751429  0.518705  0.676892]
# [-1.589405  0.807797  0.506116  0.680422]
# ```

# In[6]:


# Jacobi iterative method for solving Ax=b by bob
import numpy as np
np.set_printoptions(precision=6, suppress=True)

A=np.array([[5,1,1,1],[1,3,1,1],[1,-2,-9,1],[1,3,-2,5]])
b=np.array([-6,2,-7,3])
n=4
x0=np.zeros(n)
x1=np.zeros(n)

for iter in range(0, 5):
    for i in range(0, n):
        x1[i]=b[i]
        for j in range(0, n):
            x1[i]=x1[i]-A[i][j]*x0[j]
        x1[i]=x1[i]+A[i][i]*x0[i]
        x1[i]=x1[i]/A[i][i]
        x0[i]=x1[i]
    print(x0)


# In[15]:


# Jacobi iterative method for solving Ax=b by bob
import numpy as np
np.set_printoptions(precision=6, suppress=True)

A=np.array([[5,1,1,1],[1,3,1,1],[1,-2,-9,1],[1,3,-2,5]])
b=np.array([-6,2,-7,3])
n=4
x0=np.zeros(n)

for iter in range(0, 5):
    for i in range(0, n):
        x1 = b[i]
        for j in range(0, n):
            x1 -= A[i][j]*x0[j]
        x1 += A[i][i]*x0[i]
        x1 /= A[i][i]
        x0[i] = x1
    print(x0)


# In[ ]: