matrix_conv.py

def load_element_matrices():
    A = ()
    P = ()
    singulars = 0
    # Read matrix
    for file_index in range(0, nproc):
        fn = f"{prefix}0.{str(file_index).rjust(5, '0')}.matrix.bin"
        with open(fn, "rb") as f:
            N, size = struct.unpack('QQ', f.read(16))
            if A == ():
                A = lil_matrix((N, N), dtype=np.float64)
                P = lil_matrix((N, N), dtype=np.float64)
            # Iterate over local matrices
            for data_index in range(0, size):
                local_size = struct.unpack('Q', f.read(8))[0]
                local_data = np.fromfile(f, np.float64, local_size * local_size)
                local_data = np.reshape(local_data, (local_size, local_size))
                local_idxs = np.fromfile(f, np.uint32, local_size)
                for i in range(0, local_size):
                    for j in range(0, local_size):
                        xc = local_idxs[i]
                        yc = local_idxs[j]
                        numb = local_data[i, j]
                        if abs(numb) > 1e-10:
                            A[xc, yc] += numb
                chunksize = 8
                for i in range(0, local_size//chunksize):
                    #j = i
                    #local_piece = np.zeros(chunksize, chunksize)
                    local_piece = local_data[i*chunksize:(i+1)*chunksize,
                                      i * chunksize:(i + 1) * chunksize]
                    inverted_piece = local_piece
                    piece_idxs = list(map(lambda x: local_idxs[x], range(i*chunksize, i*chunksize+chunksize)))
                    for i in range(0, chunksize):
                        for j in range(0, chunksize):
                            numb = inverted_piece[i, j]
                            xc = piece_idxs[i]
                            yc = piece_idxs[j]
                            if abs(numb) > 1e-10:
                                P[xc, yc] += numb
                if data_index % 100 == 0:
                    print(f"Loaded chunk {data_index}, singulars: {singulars}")
        print(f"Loaded file {fn}")
    return csr_matrix(A), csr_matrix(P)