#!/usr/bin/env python
# -*- coding: utf-8 -*-

import numpy as np

matrix_2_2 = True
matrix_3_3 = False

if matrix_2_2 == True:
    comp_x = 9999.0
    matrix_a = []
    vec_solu = np.array([0, 0], 'f4')
    eq_1 = np.array([3, 2, 7], 'f4')
    eq_2 = np.array([4, 6, 16], 'f4')

    if eq_1[0] > eq_1[1] and eq_2[1] > eq_2[0]:
        pass

    else:
        print 'El sistema de ecuaciones no es una matriz diagonalmente dominante, no habrá convergencia con el método Jacobi'
        quit()

    eq_1_desp = [eq_1[2], -eq_1[1]]
    for i in range(2):
        eq_1_desp[i] = eq_1_desp[i]/eq_1[0]
    matrix_a.append(eq_1_desp)

    eq_2_desp = [eq_2[2], -eq_2[0]]
    for i in range(2):
        eq_2_desp[i] = eq_2_desp[i]/eq_2[1]
    matrix_a.append(eq_2_desp)

    matrix_a = np.asarray(matrix_a, 'f4')

    sol_x = matrix_a[0][0] + (matrix_a[0][1] * vec_solu[1])
    sol_y = matrix_a[1][0] + (matrix_a[1][1] * vec_solu[0])
    print sol_x
    print sol_y
    print
    print

    vec_solu[0] = sol_x
    vec_solu[1] = sol_y


    while 0.001 < abs(comp_x - sol_x):
        comp_x = vec_solu[0]

        sol_x = matrix_a[0][0] + (matrix_a[0][1] * vec_solu[1])
        sol_y = matrix_a[1][0] + (matrix_a[1][1] * vec_solu[0])
        print sol_x
        print sol_y
        print
        print
        vec_solu[0] = sol_x
        vec_solu[1] = sol_y

    print vec_solu

elif matrix_3_3 == True:
    comp_x = 9999.0
    matrix_a = []
    vec_solu = np.array([0, 0, 0], 'f4')
    eq_1 = np.array([4, -1, -1, 3], 'f4')
    eq_2 = np.array([-2, 6, 1, 9], 'f4')
    eq_3 = np.array([-1, 1, 7, -6], 'f4')

    if eq_1[0] > eq_1[1] + eq_1[2] and eq_2[1] > eq_2[0] + eq_2[2] and eq_3[2] > eq_3[0] + eq_3[1]:
        pass

    else:
        print 'El sistema de ecuaciones no es una matriz diagonalmente dominante, no habrá convergencia con el método Jacobi'
        quit()

    eq_1_desp = [eq_1[3], -eq_1[1], -eq_1[2]]
    for i in range(3):
        eq_1_desp[i] = eq_1_desp[i]/eq_1[0]
    matrix_a.append(eq_1_desp)

    eq_2_desp = [eq_2[3], -eq_2[0], -eq_2[2]]
    for i in range(3):
        eq_2_desp[i] = eq_2_desp[i]/eq_2[1]
    matrix_a.append(eq_2_desp)

    eq_3_desp = [eq_3[3], -eq_3[0], -eq_3[1]]
    for i in range(3):
        eq_3_desp[i] = eq_3_desp[i]/eq_3[2]
    matrix_a.append(eq_3_desp)

    matrix_a = np.asarray(matrix_a, 'f4')

    sol_x = matrix_a[0][0] + (matrix_a[0][1] * vec_solu[1]) + (matrix_a[0][2] * vec_solu[2])
    sol_y = matrix_a[1][0] + (matrix_a[1][1] * vec_solu[0]) + (matrix_a[1][2] * vec_solu[2])
    sol_z = matrix_a[2][0] + (matrix_a[2][1] * vec_solu[0]) + (matrix_a[2][2] * vec_solu[1])

    vec_solu[0] = sol_x
    vec_solu[1] = sol_y
    vec_solu[2] = sol_z


    while 0.00001 < abs(comp_x - sol_x):
        comp_x = vec_solu[0]

        sol_x = matrix_a[0][0] + (matrix_a[0][1] * vec_solu[1]) + (matrix_a[0][2] * vec_solu[2])
        sol_y = matrix_a[1][0] + (matrix_a[1][1] * vec_solu[0]) + (matrix_a[1][2] * vec_solu[2])
        sol_z = matrix_a[2][0] + (matrix_a[2][1] * vec_solu[0]) + (matrix_a[2][2] * vec_solu[1])

        vec_solu[0] = sol_x
        vec_solu[1] = sol_y
        vec_solu[2] = sol_z

    print vec_solu