WNJXYK
Thanks to the cruel world.
WNJXYKのBlog
Sequential Block-Based Gauss-Jordan Algorithm
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.

使用Python语言对序列的分块高斯算法的一个实现,可以用来做矩阵求逆的并行计算。

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)

    # Output & Check Answer
    print("Correct:", check_matrix(I * A, numpy.matrix(numpy.eye(siz))))

if __name__=='__main__': main(128, 2)
赞赏
https://secure.gravatar.com/avatar/f83b57c055136369e9feba5d6671d6b5?s=256&r=g

WNJXYK

文章作者

一个蒟蒻

发表评论

textsms
account_circle
email

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