import numpy as np

def gaussian_elimination_with_pivoting(A, b):
   
    n = len(b)
    # Augmented matrix
    augmented_matrix = np.hstack((A, b.reshape(-1, 1)))
    
    # Forward elimination
    for i in range(n):
        # Partial pivoting: Find the maximum element in the current column
        max_row = np.argmax(np.abs(augmented_matrix[i:, i])) + i
        if augmented_matrix[max_row, i] == 0:
            raise ValueError("Matrix is singular or nearly singular.")
        
        # Swap rows
        augmented_matrix[[i, max_row]] = augmented_matrix[[max_row, i]]
        
        # Eliminate entries below the pivot
        for j in range(i + 1, n):
            factor = augmented_matrix[j, i] / augmented_matrix[i, i]
            augmented_matrix[j, i:] -= factor * augmented_matrix[i, i:]
    
    # Back substitution
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (augmented_matrix[i, -1] - np.dot(augmented_matrix[i, i + 1:n], x[i + 1:])) / augmented_matrix[i, i]
    
    return x

# Example usage
A = np.array([[2, -1, 1],
              [1, 3, 2],
              [1, -1, 2]], dtype=float)
b = np.array([2, 9, -3], dtype=float)

x = gaussian_elimination_with_pivoting(A, b)
print("Solution vector x:", x)
