```import random
from math import sqrt
from utilities import mat_mult, mat_dot, mat_minus
from copy import deepcopy
import sys
# sys.path.append('/Users/jim/Python_3/math_libraries')
sys.path.append('/home/pi/python3/math_libraries')
from number_theory import DataError

def positive_definite_system(verbose, a, rhs, scaling_limit, lam_iters, mat_iters):
n = len(a)
rng_n = range(n)
rescale, scales = do_rescale(verbose, a, rhs, rng_n, n, scaling_limit)
lam = largest_eigenvalue(a, rng_n, lam_iters)
if verbose:
print('\nAfter {0:d} iterations:\n   Largest eigenvalue <= {1:.1f}'.format(lam_iters, lam))
normalize_max_lamda(a, rhs, lam, rng_n)
k = build_k_matrix(a, rhs, n, rng_n)
pm, pmm1 = initialize_solution_vectors(rng_n, n, rhs)
mag_rhs = magnitude_of_rhs(rhs)
sol = iteration(verbose, a, rhs, rescale, scales, mat_iters, rng_n, k, pm, pmm1, mag_rhs)

return sol

def scale_iteration(verbose, it_num, a, rhs, rng_n, rescale, scales, pmp1, mag_rhs):
sol = []
for i in rng_n:
sol.append(pmp1[i][0] / pmp1[-2][0])
if rescale:
for i in rng_n:
if scales[i] != 1.0:
sol[i] /= scales[i]

if verbose:
print('\nIterate {0:d}:'.format(it_num))
for i in rng_n:
print('   Variable {0:d} = {1:7.5f}'.format((i + 1), sol[i]))
print('\n   Relative error within 5% if matrix\'s condition number is < {0:.1f}'.format(((it_num + 2) / 2.56)**2))
print actual_relative_error(a, rhs, rng_n, sol, mag_rhs)

return sol

# Rescale columns (and rows, for symmetry), if necessary
def do_rescale(verbose, a, rhs, rng_n, n, scaling_limit):
column_sq_total, column_sums = 0.0, []
for j in rng_n:
sum = 0.0
for i in rng_n:
sum += a[i][j]**2
column_sq_total += sum
column_sums.append(sum)

column_sq_average = column_sq_total / n

rescale, scales = False, []
for i in rng_n:
temp = column_sums[i] / column_sq_average
if temp > scaling_limit or temp < 1.0 / scaling_limit:
if verbose:
print('Column {0:d} needs to be rescaled'.format(i))
# Use square root of temp so that squared sum of column
# coefficient is rescaled to column_sq_average
scales.append(sqrt(temp))
rescale = True
else:
scales.append(1.0) # sentinel

if rescale:
if verbose:
print('Rescaling')
for k in rng_n:
if scales[k] != 1.0:
s_k = scales[k]
# rows first
for j in rng_n:
a[k][j] /= s_k
# right hand side
rhs[k] /= s_k
# columns next
for i in rng_n:
a[i][k] /= s_k
else:
if verbose:
print('Rescaling unnecessary with scaling limit = {0:d}'.format(scaling_limit))

return rescale, scales

def largest_eigenvalue(a, rng_n, lam_iters):
# Initial vector of random values for starting power method
b0 = []
for _ in rng_n:
b0.append([random.random()])

try:
b = []
last_b = mat_mult(a, b0)
for _ in range(lam_iters - 1):
b.append(last_b)
last_b = mat_mult(a, last_b)

b.append(last_b)
bmp2, bmp1, bm = b[-1], b[-2], b[-3]
c0 = mat_dot(bm, bmp1)
c1 = mat_dot(bmp1, bmp1)
c2 = mat_dot(bmp1, bmp2)
c3 = mat_dot(bmp2, bmp2)

r1, r2, r3 = c1 / c0, c2 / c1, c3 / c2
lam = r3
num, den = 2 * (r3 - r2)**2, 2 * r2 - r1 - r3
if num != 0.0 and den != 0.0:
lam += num / den
except DataError as e:
print('DataError exception occurred: {0:s}'.format(e.value))

return lam

def normalize_max_lamda(a, rhs, lam, rng_n):
# Normalize a so maximum eigenvalue is < 1
for row in a:
for j in rng_n:
row[j] /= lam

# Normalize rhs too
for i in rng_n:
rhs[i] /= lam

def build_k_matrix(a, rhs, n, rng_n):
k = []
for i in rng_n:
row = []
for j in range(n + 2):
if j < n:
row.append(-4.0 * a[i][j])
if i == j:
row[-1] += 2.0
elif j == n:
row.append(4.0 * rhs[i])
else:
row.append(0.0)
k.append(row)

row = [0.0] * n
row.append(2.0)
row.append(2.0)
k.append(row)

next_row = deepcopy(row)
next_row[n] = 0.0
k.append(next_row)

return k

def initialize_solution_vectors(rng_n, n, rhs):
pmm1 = [[0.0]] * n
pmm1.append([1.0])
pmm1.append([1.0])

pm = []
for i in rng_n:
pm.append([4.0 * rhs[i]])

pm.append([4.0])
pm.append([1.0])

return pm, pmm1

def magnitude_of_rhs(rhs):
sum = 0.0
for el in rhs:
sum += el**2

return sqrt(sum)

def iteration(verbose, a, rhs, rescale, scales, mat_iters, rng_n, k, pm, pmm1, mag_rhs):
for i in range(mat_iters):
temp_pm = mat_mult(k, pm)
pmp1 = mat_minus(temp_pm, pmm1)
pmm1, pm = pm, pmp1
# pmm1 = pm
# pm = pmp1
sol = scale_iteration(verbose, i + 1, a, rhs, rng_n, rescale, scales, pmp1, mag_rhs)

return sol

def actual_relative_error(a, rhs, rng_n, sol, mag_rhs):
vec_sol = [[el] for el in sol]
approx_rhs = mat_mult(a, vec_sol)
delt = 0.0
for i in rng_n:
delt += (rhs[i] - approx_rhs[i][0])**2

return '   Actual relative error = {0:.3f}%'.format(100 * sqrt(delt) / mag_rhs)

if __name__ == "__main__":
import sys
a, rhs = [[3.0, 3.0, 5.0], [3.0, 5.0, 9.0], [5.0, 9.0, 17.0]], [24.0, 40.0, 74.0]
x = [1.0, 2.0, 3.0]
lam_iters, mat_iters, scaling_limit = 30, 15, 10
rng_n, verbose = range(len(a)), True
sol = positive_definite_system(verbose, a, rhs, scaling_limit, lam_iters, mat_iters)
print('\n')
for i in rng_n:
print('x{0:d} = {1:7.5f} (exact {2:7.5f})'.format((i + 1), sol[i], x[i]))

# Typical results:
# ...
# Iterate 15:
#    Variable 1 = 1.03584
#    Variable 2 = 1.81987
#    Variable 3 = 3.07039
#
#    Relative error within 5% if matrix's condition number is < 44.1
#    Actual relative error = 0.347%

# x1 = 1.03584 (exact 1.00000)
# x2 = 1.81987 (exact 2.00000)
# x3 = 3.07039 (exact 3.00000)
```