Untitled

def ensure_nonzero_diagonal(matrix, size, index):
    """
    Ensures that the diagonal element matrix[index, index] is non-zero. If zero, 
    swaps rows and columns with another non-zero element in the same row.
    """
    # If the diagonal element is already non-zero, return the matrix unchanged.
    if matrix[index][index] != 0:
        return matrix, np.eye(size)

    col = index + 1
    # Search for a non-zero element in the row.
    while col < size and matrix[index][col] == 0:
        col += 1

    # Construct an identity matrix to perform the swap.
    transform = np.eye(size)
    transform[col][index] = 1

    # Apply row and column swaps.
    matrix = transform.T @ matrix @ transform

    return matrix, transform

def create_transformation(matrix, size, index):
    """
    Creates a transformation matrix to eliminate off-diagonal elements in the row.
    """
    transform = np.eye(size)
    diag_element = matrix[index][index]

    # Compute transformation values for off-diagonal elements in the row.
    for col in range(index + 1, size):
        transform[index][col] = -matrix[index][col] / diag_element

    return transform

def symmetric_diagonalization(matrix):
    """
    Diagonalizes a symmetric matrix by sequentially eliminating off-diagonal entries.
    """
    # Validate symmetry of the matrix.
    if not np.array_equal(matrix, matrix.T):
        return "Error: Input matrix is not symmetric."

    size = len(matrix)
    combined_transform = np.eye(size)

    # Iterate through rows to diagonalize the matrix.
    for row in range(size - 1):
        matrix, temp_transform = ensure_nonzero_diagonal(matrix, size, row)
        combined_transform = combined_transform @ temp_transform

        row_transform = create_transformation(matrix, size, row)
        combined_transform = combined_transform @ row_transform

        matrix = row_transform.T @ matrix @ row_transform

    return matrix, combined_transform

# Example input matrix and function call
input_matrix = np.array([
    [1, 0, 0, 0],
    [0, -4, 0, 0],
    [0, 0, 0, -1/2],
    [0, 0, -0.5, 16]
])

diagonalized, transformation_matrix = symmetric_diagonalization(input_matrix)
print("Original matrix: \n", input_matrix)
print("Diagonalized matrix: \n", diagonalized)
print("Transformation matrix: \n", transformation_matrix)