-- True
--
qr :: forall m n a. (Arity m, Arity n, Eq a, Algebraic.C a, Ring.C a)
- => Mat m n a -> (Mat m m a, Mat m n a)
+ => Mat (S m) (S n) a
+ -> (Mat (S m) (S m) a, Mat (S m) (S n) a)
qr matrix =
ifoldl col_function initial_qr columns
where
ifoldl (f col_idx) (q,r) col
-- | Process the entries in a column, doing basically the same
- -- thing as col_dunction does. It updates the QR factorization,
+ -- thing as col_function does. It updates the QR factorization,
-- maybe, and returns the current one.
f col_idx (q,r) idx _ -- ignore the current element
| idx <= col_idx = (q,r) -- leave it alone
| otherwise = (q*rotator, (transpose rotator)*r)
where
y = r !!! (idx, col_idx)
- rotator :: Mat m m a
+ rotator :: Mat (S m) (S m) a
rotator = givens_rotator col_idx idx (r !!! (col_idx, col_idx)) y
-- iterated QR algorithm (see Golub and Van Loan, \"Matrix
-- Computations\").
--
+-- Warning: this may not converge if there are repeated eigenvalues
+-- (in magnitude).
+--
-- Examples:
--
-- >>> import Linear.Matrix ( Col2, Col3, Mat2, Mat3 )
=> Int
-> Mat (S m) (S m) a
-> Col (S m) a
-eigenvalues iterations matrix =
- diagonal (ut_approximation iterations)
- where
- ut_approximation :: Int -> Mat (S m) (S m) a
- ut_approximation 0 = matrix
- ut_approximation k = rk*qk where (qk,rk) = qr (ut_approximation (k-1))
+eigenvalues iterations matrix
+ | iterations < 0 = error "negative iterations requested"
+ | iterations == 0 = diagonal matrix
+ | otherwise =
+ diagonal (ut_approximation (iterations - 1))
+ where
+ ut_approximation :: Int -> Mat (S m) (S m) a
+ ut_approximation 0 = matrix
+ ut_approximation k = ut_next
+ where
+ ut_prev = ut_approximation (k-1)
+ (qk,rk) = qr ut_prev
+ ut_next = rk*qk
+
-- | Compute the eigenvalues and eigenvectors of a symmetric matrix
-- references see Goluv and Van Loan, \"Matrix Computations\", or
-- \"Calculation of Gauss Quadrature Rules\" by Golub and Welsch.
--
+-- Warning: this may not converge if there are repeated eigenvalues
+-- (in magnitude).
+--
-- Examples:
--
-- >>> import Linear.Matrix ( Col2, Col3, Mat2, Mat3 )
--- >>> import Linear.Matrix ( column', frobenius_norm, fromList )
+-- >>> import Linear.Matrix ( column, frobenius_norm, fromList )
-- >>> import Linear.Matrix ( identity_matrix, vec3d )
-- >>> import Normed ( Normed(..) )
--
-- >>> let v1 = (1 / (norm v1') :: Double) *> v1'
-- >>> let v2' = vec3d (-4, -2, 5) :: Col3 Double
-- >>> let v2 = (1 / (norm v2') :: Double) *> v2'
--- >>> frobenius_norm ((column' vecs 0) - v0) < 1e-12
+-- >>> frobenius_norm ((column vecs 0) - v0) < 1e-12
-- True
--- >>> frobenius_norm ((column' vecs 1) - v1) < 1e-12
+-- >>> frobenius_norm ((column vecs 1) - v1) < 1e-12
-- True
--- >>> frobenius_norm ((column' vecs 2) - v2) < 1e-12
+-- >>> frobenius_norm ((column vecs 2) - v2) < 1e-12
-- True
--
eigenvectors_symmetric :: forall m a. (Arity m, Algebraic.C a, Eq a)
where
(t_prev, p_prev) = tp_pair (k-1)
(qk,rk) = qr t_prev
- tk = rk*qk
pk = p_prev*qk
-
+ tk = rk*qk
(values, vectors) = (first diagonal) (tp_pair iterations)