﻿ Sequential Block-Based Gauss-Jordan Algorithm – WNJXYKのBlog
WNJXYK
Thanks to the cruel world.
Sequential Block-Based Gauss-Jordan Algorithm

A implement of sequential block-based Gauss-Jordan algorithm by Python which was mentioned in the Papler 《Parallel processing on block-based Gauss-Jordan algorithm for desktop grid》. This algorithm can be adapted for calculating matrix inversion with parallel processes.

import numpy, math
import sys, time

# Divided Matrix into Matrix Blocks
def block_divide(matrix, blk, blc):
return [[numpy.matrix([[matrix[i*blk+ii, j*blk+jj] for jj in range(blk)] for ii in range(blk)]) for j in range(blc)] for i in range(blc)]

# Merge Matrix Blocks to  Matrix
def block_merge(matrix, blk, blc):
return numpy.matrix([[matrix[i][j][ii, jj] for j in range(blc) for jj in range(blk) ] for i in range(blc) for ii in range(blk) ])

# Check A == B
def check_matrix(A, B):
EPS = 1e-6
A, B = numpy.array(A).flatten(), numpy.array(B).flatten()
if len(A) != len(B): return False
for i in range(len(A)):
if math.fabs(A[i]-B[i])>EPS: return False
return True

# Sequential Block-based Gauss-Jordan Algorithm
def BlockBased_GaussJordan(A, B, blk, blc):
lastA, lastB = A, B
for k in range(blc):
A = [[numpy.matrix(numpy.eye(blk))]*blc for i in range(blc)]
B = [[numpy.matrix(numpy.eye(blk))]*blc for i in range(blc)]
A[k][k] = lastA[k][k].I
B[k][k] = A[k][k]
for j in range(k+1, blc): A[k][j] = A[k][k]*lastA[k][j]
for j in range(k): B[k][j] = A[k][k]*lastB[k][j]
for j in range(k+1, blc):
for i in range(blc):
if i==k: continue
A[i][j] = lastA[i][j]-lastA[i][k]*A[k][j]
for j in range(k):
for i in range(blc):
if i==k: continue
B[i][j] = lastB[i][j]-lastA[i][k]*B[k][j]
for i in range(blc):
if i==k: continue
B[i][k] = -lastA[i][k]*A[k][k]
lastA, lastB = A, B
return B

def main(siz, blk):
# Assert siz%blk == 0
if siz%blk !=0:
print("Make Sure Block_Size(%d) | Matrix_Size(%d)!" % (blk, siz))
return

# Calculte Matrix Block Size
blc = int(siz/blk)

# Get A Random Matrix A
A = numpy.matrix(numpy.random.rand(siz, siz))
Ab = block_divide(A, blk, blc)

# Get An Eye Matrix B
B = numpy.matrix(numpy.eye(siz, siz))
Bb = block_divide(B, blk, blc)

# Calculate A.Inverse By Block-Based Gauss-Jordan Algorithm
Ib = BlockBased_GaussJordan(Ab, Bb, blk, blc)
I = block_merge(Ib, blk, blc)

print("Correct:", check_matrix(I * A, numpy.matrix(numpy.eye(siz))))

if __name__=='__main__': main(128, 2)


### 推荐文章

WNJXYKのBlog

Sequential Block-Based Gauss-Jordan Algorithm
A implement of sequential block-based Gauss-Jordan algorithm by Python which was mentioned in the Papler 《Parallel processing on block-based Gauss-Jordan algorithm for deskt…

2018-11-03
<--! http2https -->