1
1
from time import time
2
2
import numpy as np
3
3
import scipy
4
+ import scipy .sparse as sp
4
5
from mpi4py import MPI
5
6
import sys
6
7
import numpy as np
7
8
import argparse
8
- from pyclassify .utils import read_config , poisson_2d_structure
9
+ from pyclassify .utils import read_config , poisson_2d_structure , make_symmetric
9
10
from pyclassify .eigenvalues import Lanczos_PRO
10
11
11
12
@@ -58,27 +59,26 @@ def compute_eigvals(A, n_procs):
58
59
density = kwargs ["density" ]
59
60
n_procs = kwargs ["n_processes" ]
60
61
61
-
62
- A = poisson_2d_structure (dim )
63
- A_np = A .toarray ()
62
+ # You could use (for low values of dim, else accuracy suffers):
63
+ # A = poisson_2d_structure(dim)
64
+ # A_np = A.toarray()
64
65
65
66
# Alternatively, consider for instance:
66
- # A = sp.random(dim, dim, density=density, format="csr") # uncomment if you want to use a random matrix instead
67
- # or
68
- # eig = np.arange(1, dim + 1)
69
- # A = np.diag(eig)
70
- # U = scipy.stats.ortho_group.rvs(dim)
67
+ eig = np .arange (1 , dim + 1 )
68
+ A = np .diag (eig )
69
+ U = scipy .stats .ortho_group .rvs (dim )
71
70
72
- # A = U @ A @ U.T
73
- # A = make_symmetric(A)
74
- # A_sp = sp.csr_matrix(A)
71
+ A = U @ A @ U .T
72
+ A = make_symmetric (A )
73
+ A_np = A
75
74
76
75
print ("---------------\n Calling Lanczos a first time to compile it..." )
77
76
Q , diag , off_diag = Lanczos_PRO (A_np , np .ones_like (np .diag (A_np )) * 1.0 )
78
77
print ("Done! Now we compute the eigenvalues.\n ---------------" )
79
78
80
79
eigvals , eigvecs , delta_t , total_mem_children = compute_eigvals (A_np , n_procs )
81
- exact_eigvals , exact_eigvecs = np .linalg .eig (A .toarray ())
80
+ print ("Now computing eigenvalues using np" )
81
+ exact_eigvals , exact_eigvecs = np .linalg .eig (A_np )
82
82
print ("---------------" )
83
83
sorted_indices = np .argsort (exact_eigvals )
84
84
exact_eigvals = exact_eigvals [sorted_indices ]
@@ -87,5 +87,6 @@ def compute_eigvals(A, n_procs):
87
87
max_error = np .max (np .abs (exact_eigvals - eigvals ))
88
88
print (f"The maximum error between real and computed eigenvalues is { max_error } " )
89
89
90
+
90
91
if max_error < 1e-8 :
91
92
print ("Pretty small, huh?" )
0 commit comments