# leto/math--matrixreal

Fetching contributors…
Cannot retrieve contributors at this time
130 lines (118 sloc) 3.86 KB
 use Test::More tests => 13; use File::Spec; use lib File::Spec->catfile("..","lib"); use Math::MatrixReal; do 'funcs.pl'; my \$DEBUG2 = 0; my \$string = "[ 1 2 3 ]\n[ 2 2 -1 ]\n[ 1 1 1 ]\n"; my \$matrix33 = Math::MatrixReal->new_from_string(\$string); print "\$matrix33" if \$DEBUG; # # Tests eigenvalues & eigenvectors computation # # # First the tridiagonal reduction (Householder) on the 3x3 # my \$symm = \$matrix33 + ~\$matrix33; ok( not(\$matrix33->is_symmetric())); ok( \$symm->is_symmetric(), 'A + ~A is symmetric'); print "Matrix 3x3 for eigenvalues & eigenvectors computation:\n\$symm" if \$DEBUG; print "Householder reduction...\n" if \$DEBUG; my (\$T, \$Q) = \$symm->householder(); print "T=\n\$T Q=\n\$Q" if \$DEBUG2; print "Is Q orthogonal?\n" if \$DEBUG; print (\$Q * ~\$Q) if \$DEBUG2; ok_matrix_orthogonal(\$Q); ok_matrix( \$symm, \$Q * \$T * ~\$Q, 'symmetric householder reduction works'); print "Diagonalization of tridiagonal...\n" if \$DEBUG; my (\$L1, \$V1) = \$T->tri_diagonalize(\$Q); print "eigenvalues L:\n\$L1 eigenvectors:\n\$V1" if \$DEBUG2; ok_eigenvectors(\$symm, \$L1, \$V1); # Get first eigenvector my \$aev1 = \$V1->column(1); my \$al1 = \$L1->element(1,1); my \$ap1_1 = \$symm * \$aev1; # A * x my \$ap1_2 = \$al1 * \$aev1; # lambda *x print "Original computation of A*ev1:\n\$ap1_1 Scaled eigenvector:\n\$ap1_2" if \$DEBUG2; ok_matrix( \$ap1_1, \$ap1_2, 'eigenvectors match'); print "Direct diagonalization...\n" if \$DEBUG; my (\$L12, \$V12) = \$symm->sym_diagonalize(); print "eigenvalues L:\n\$L12 eigenvectors:\n\$V12" if \$DEBUG2; ok_eigenvectors(\$symm, \$L12, \$V12); ok_matrix_orthogonal(\$V12); # Double check the equality ok_matrix( \$L12, \$L1); ok_matrix( \$V12, \$V1); # # Now test the eigenvalues only computations... # print "Recomputing: Eigenvalues only.\n 3x3\n" if \$DEBUG; my \$altT = \$symm->householder_tridiagonal(); ok_matrix( \$altT, \$T,'householder_tridiagonal works'); my \$altL1 = \$altT->tri_eigenvalues(); ok_matrix( \$altL1, \$L1,'tri_eigenvalues works'); my \$altL12 = \$symm->sym_eigenvalues(); ok_matrix( \$altL12, \$L12, 'sym_eigenvalues works'); __END__ # Attempt: # Obtain the eigenvectors when eigenvalues are known # using inverse iteration. # We solve (M - lambda * I) * b(k+1) = b(k) # with b(0) a random unit vector and b(k) is # normalized at each step. # This should converge towards the eigenvector, # but there are problems: # - for some value, there is not convergence # (the above system is rather singular, so...) # - the solution can oscillate between v and -v (?) # Rodolphe Ortalo, 99/06/14 sub obtain_eigenvector (\$\$) { my (\$M, \$eigenvalue) = @_; # Form the linear system A - lamda1 * I my \$inv_it = \$M->shadow(); \$inv_it->one(); \$inv_it->multiply_scalar(\$inv_it, (-1.0 * \$eigenvalue)); \$inv_it->add(\$M, \$inv_it); print "Linear system matrix:\n \$inv_it" if \$DEBUG2; # Creates a random vector my (\$rows, \$cols) = \$inv_it->dim(); my \$b = Math::MatrixReal->new(\$rows, 1); for (my \$i = 1; \$i <= \$rows; \$i++) { \$b->assign(\$i, 1, rand()); } # Normalize it my \$l = \$b->length(); \$b->multiply_scalar(\$b, (1.0 / \$l)); # Now do LR decomposition for linear system my \$inv_it_LR = \$inv_it->decompose_LR(); # Check iterations my \$iter = 0; my \$delta; do { my (\$dim, \$b_base, \$base) = \$inv_it_LR->solve_LR(\$b); # Normalize my \$l = \$b_base->length(); \$b_base->multiply_scalar(\$b_base, (-1.0 / \$l)); # print "b_base=\n\$b_base"; \$b->subtract(\$b_base,\$b); \$delta = \$b->norm_one(); print "delta=\$delta\n"; \$b = \$b_base; } while ((\$delta >= 1e-10) && (\$iter++ <= 10)); return \$b; } # # Now, try to find one eigenvector again... # (Using Steffen's functions...:-) # my \$ev = obtain_eigenvector(\$symm, \$al1); print "Real ev:\n \$aev1 Found ev:\n \$ev" if \$DEBUG; ok_matrix(15, \$ev, \$aev1); ok_matrix(16, obtain_eigenvector(\$symm, \$L1->element(2,1)), \$V1->column(2)); ok_matrix(17, obtain_eigenvector(\$symm, \$L1->element(3,1)), \$V1->column(3));