11from time import time
22import numpy as np
33import scipy
4+ import scipy .sparse as sp
45from mpi4py import MPI
56import sys
67import numpy as np
78import argparse
8- from pyclassify .utils import read_config , poisson_2d_structure
9+ from pyclassify .utils import read_config , poisson_2d_structure , make_symmetric
910from pyclassify .eigenvalues import Lanczos_PRO
1011
1112
@@ -58,27 +59,26 @@ def compute_eigvals(A, n_procs):
5859density = kwargs ["density" ]
5960n_procs = kwargs ["n_processes" ]
6061
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()
6465
6566# 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 )
7170
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
7574
7675print ("---------------\n Calling Lanczos a first time to compile it..." )
7776Q , diag , off_diag = Lanczos_PRO (A_np , np .ones_like (np .diag (A_np )) * 1.0 )
7877print ("Done! Now we compute the eigenvalues.\n ---------------" )
7978
8079eigvals , 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 )
8282print ("---------------" )
8383sorted_indices = np .argsort (exact_eigvals )
8484exact_eigvals = exact_eigvals [sorted_indices ]
@@ -87,5 +87,6 @@ def compute_eigvals(A, n_procs):
8787max_error = np .max (np .abs (exact_eigvals - eigvals ))
8888print (f"The maximum error between real and computed eigenvalues is { max_error } " )
8989
90+
9091if max_error < 1e-8 :
9192 print ("Pretty small, huh?" )
0 commit comments