248 lines
6.8 KiB
Python
248 lines
6.8 KiB
Python
import numpy as np
|
|
import multiprocessing as mp
|
|
import time
|
|
|
|
def sequential_determinant(matrix):
|
|
n = matrix.shape[0]
|
|
|
|
A = matrix.astype(np.float64).copy()
|
|
|
|
sign = 1
|
|
|
|
for i in range(n):
|
|
max_row = i
|
|
max_val = abs(A[i, i])
|
|
|
|
for k in range(i + 1, n):
|
|
if abs(A[k, i]) > max_val:
|
|
max_row = k
|
|
max_val = abs(A[k, i])
|
|
|
|
if max_row != i:
|
|
A[[i, max_row]] = A[[max_row, i]]
|
|
sign *= -1
|
|
|
|
if abs(A[i, i]) < 1e-15:
|
|
return 0.0, -np.inf
|
|
|
|
for k in range(i + 1, n):
|
|
factor = A[k, i] / A[i, i]
|
|
A[k, i:] -= factor * A[i, i:]
|
|
|
|
log_abs_det = 0.0
|
|
for i in range(n):
|
|
if A[i, i] < 0:
|
|
sign *= -1
|
|
log_abs_det += np.log(abs(A[i, i]))
|
|
|
|
return sign, log_abs_det
|
|
|
|
def to_float(sign, log_abs_det):
|
|
if log_abs_det == -np.inf:
|
|
return 0.0
|
|
|
|
if log_abs_det > 700:
|
|
return sign * np.inf
|
|
|
|
return sign * np.exp(log_abs_det)
|
|
|
|
|
|
def worker_lu_decomposition(args):
|
|
matrix, start, end = args
|
|
n = matrix.shape[0]
|
|
A = matrix.copy()
|
|
perm_count = 0
|
|
|
|
for i in range(start, min(end, n)):
|
|
max_row = i
|
|
for k in range(i + 1, n):
|
|
if abs(A[k, i]) > abs(A[max_row, i]):
|
|
max_row = k
|
|
|
|
if max_row != i:
|
|
A[[i, max_row]] = A[[max_row, i]]
|
|
perm_count += 1
|
|
|
|
if abs(A[i, i]) > 1e-15:
|
|
for k in range(i + 1, n):
|
|
factor = A[k, i] / A[i, i]
|
|
A[k, i:] -= factor * A[i, i:]
|
|
|
|
return A, perm_count
|
|
|
|
|
|
def parallel_determinant(matrix, num_processes = None):
|
|
if num_processes is None:
|
|
num_processes = mp.cpu_count()
|
|
|
|
n = matrix.shape[0]
|
|
|
|
if n < 100 or num_processes == 1:
|
|
return sequential_determinant(matrix)
|
|
|
|
A = matrix.astype(np.float64).copy()
|
|
sign = 1
|
|
|
|
with mp.Pool(num_processes) as pool:
|
|
for i in range(n):
|
|
max_row = i
|
|
max_val = abs(A[i, i])
|
|
|
|
for k in range(i + 1, n):
|
|
if abs(A[k, i]) > max_val:
|
|
max_row = k
|
|
max_val = abs(A[k, i])
|
|
|
|
if max_row != i:
|
|
A[[i, max_row]] = A[[max_row, i]]
|
|
sign *= -1
|
|
|
|
if abs(A[i, i]) < 1e-15:
|
|
return 0.0, -np.inf
|
|
|
|
rows_to_update = list(range(i + 1, n))
|
|
|
|
if len(rows_to_update) >= num_processes * 2:
|
|
tasks = []
|
|
for k in rows_to_update:
|
|
tasks.append((A[k].copy(), A[i].copy(), A[k, i], A[i, i], i))
|
|
|
|
updated_rows = pool.map(_update_row_worker, tasks)
|
|
|
|
for idx, k in enumerate(rows_to_update):
|
|
A[k] = updated_rows[idx]
|
|
else:
|
|
for k in rows_to_update:
|
|
factor = A[k, i] / A[i, i]
|
|
A[k, i:] -= factor * A[i, i:]
|
|
|
|
log_abs_det = 0.0
|
|
for i in range(n):
|
|
if A[i, i] < 0:
|
|
sign *= -1
|
|
log_abs_det += np.log(abs(A[i, i]))
|
|
|
|
return sign, log_abs_det
|
|
|
|
|
|
def _update_row_worker(args):
|
|
row_k, row_i, a_ki, a_ii, i = args
|
|
|
|
import os
|
|
os.environ['OMP_NUM_THREADS'] = '1'
|
|
os.environ['MKL_NUM_THREADS'] = '1'
|
|
|
|
factor = a_ki / a_ii
|
|
|
|
row_k[i:] -= factor * row_i[i:]
|
|
|
|
return row_k
|
|
|
|
|
|
def generate_stable_matrix(size, scale = 1.0) -> np.ndarray:
|
|
np.random.seed(13)
|
|
|
|
random_matrix = np.random.randn(size, size)
|
|
Q, R = np.linalg.qr(random_matrix)
|
|
|
|
eigenvalues = np.random.uniform(0.5, 2.0, size) * scale
|
|
D = np.diag(eigenvalues)
|
|
|
|
A = Q @ D @ Q.T
|
|
|
|
return A
|
|
|
|
|
|
def benchmark(sizes, thread_counts):
|
|
results = []
|
|
|
|
print("=" * 80)
|
|
print("БЕНЧМАРК ВЫЧИСЛЕНИЯ ДЕТЕРМИНАНТА")
|
|
print("=" * 80)
|
|
print()
|
|
|
|
for size in sizes:
|
|
print(f"\n{'=' * 80}")
|
|
print(f"Размер матрицы: {size}x{size}")
|
|
print(f"{'=' * 80}")
|
|
|
|
matrix = generate_stable_matrix(size)
|
|
|
|
try:
|
|
sign_numpy, logdet_numpy = np.linalg.slogdet(matrix)
|
|
det_numpy = sign_numpy * np.exp(logdet_numpy) if logdet_numpy < 700 else sign_numpy * np.inf
|
|
print(f"Детерминант (numpy): {det_numpy:.6e}")
|
|
print(f" Sign: {sign_numpy}, Log|det|: {logdet_numpy:.6f}")
|
|
except Exception as e:
|
|
sign_numpy, logdet_numpy = None, None
|
|
det_numpy = None
|
|
print(f"Детерминант (numpy): Ошибка: {e}")
|
|
|
|
print()
|
|
|
|
start_time = time.time()
|
|
sign_seq, logdet_seq = sequential_determinant(matrix)
|
|
seq_time = time.time() - start_time
|
|
det_seq = to_float(sign_seq, logdet_seq)
|
|
|
|
print(f"Последовательный алгоритм:")
|
|
print(f" Время: {seq_time:.4f} сек")
|
|
print(f" Детерминант: {det_seq:.6e}")
|
|
print(f" Sign: {sign_seq}, Log|det|: {logdet_seq:.6f}")
|
|
|
|
results.append({
|
|
'size': size,
|
|
'threads': 1,
|
|
'time': seq_time,
|
|
'determinant': det_seq,
|
|
'speedup': 1.0
|
|
})
|
|
|
|
print()
|
|
|
|
for num_threads in thread_counts:
|
|
if num_threads == 1:
|
|
continue
|
|
|
|
start_time = time.time()
|
|
sign_par, logdet_par = parallel_determinant(matrix, num_threads)
|
|
par_time = time.time() - start_time
|
|
det_par = to_float(sign_par, logdet_par)
|
|
|
|
speedup = seq_time / par_time if par_time > 0 else 0
|
|
|
|
print(f"Параллельный алгоритм ({num_threads} потоков):")
|
|
print(f" Время: {par_time:.4f} сек")
|
|
print(f" Детерминант: {det_par:.6e}")
|
|
print(f" Ускорение: {speedup:.2f}x")
|
|
|
|
results.append({
|
|
'size': size,
|
|
'threads': num_threads,
|
|
'time': par_time,
|
|
'determinant': det_par,
|
|
'speedup': speedup
|
|
})
|
|
|
|
print()
|
|
|
|
return results
|
|
|
|
if __name__ == "__main__":
|
|
sizes = [100, 300, 500]
|
|
thread_counts = [1, 2, 4, mp.cpu_count()]
|
|
|
|
thread_counts = sorted(list(set(thread_counts)))
|
|
|
|
results = benchmark(sizes, thread_counts)
|
|
|
|
print("\n" + "=" * 80)
|
|
print("СВОДНАЯ ТАБЛИЦА РЕЗУЛЬТАТОВ")
|
|
print("=" * 80)
|
|
print()
|
|
print(f"{'Размер':<10} {'Потоки':<10} {'Время (сек)':<15} {'Ускорение':<12}")
|
|
print("-" * 80)
|
|
|
|
for r in results:
|
|
print(f"{r['size']:<10} {r['threads']:<10} {r['time']:<15.4f} {r['speedup']:<12.2f}")
|