Functions | |
| template<typename value_type > | |
| void | get_rotation (value_type const &a, value_type const &b, value_type &c, value_type &s) |
| template<typename value_type > | |
| void | set_rotation (value_type &a, value_type &b, value_type const &c, value_type const &s) |
| template<typename size_type , typename matrix_type , typename vector_type > | |
| void | hessenberg_matrix_transform (size_type j, matrix_type &H, vector_type &g, vector_type &c, vector_type &s) |
| template<typename size_type , typename vector_type > | |
| bool | is_orthonormal (size_type m, std::vector< vector_type > const &v) |
| template<typename size_type , typename matrix_type > | |
| bool | is_upper_triangular (size_type m, matrix_type const &H) |
| template<typename size_type , typename vector_type , typename matrix_type > | |
| void | update (vector_type &x, size_type m, matrix_type &H, vector_type const &g, std::vector< vector_type > const &v) |
| void OpenTissue::math::big::detail::get_rotation | ( | value_type const & | a, | |
| value_type const & | b, | |||
| value_type & | c, | |||
| value_type & | s | |||
| ) | [inline] |
This method determines the rotation around the z-axis, which rotates the vector [a,b]^T into a vector with only a positive x-component.
| a | The x-component of the vector to be rorated. | |
| b | The y-component of the vector to be rotated. | |
| c | Upon return this argument holds the value of cos(theta), where theta is the angle the vector needs to be rotated around the z-axis. | |
| s | Upon return this argument holds the value of sin(theta), where theta is the angle the vector needs to be rotated around the z-axis. |
| void OpenTissue::math::big::detail::hessenberg_matrix_transform | ( | size_type | j, | |
| matrix_type & | H, | |||
| vector_type & | g, | |||
| vector_type & | c, | |||
| vector_type & | s | |||
| ) | [inline] |
This method transforms least squares subproblem into upper triangular form. Below follows an outline of how this works.
This method is invoked after the Arnoldi method has generated the next vector for the span of the Krylow subspace.
The Arnoldi method results in the factorization
A V_M = V_{M+1} H_{M}
Where H_{M} is a Hessenberg matrix of dimension (M+1 M). The columns of V_M+1 is an orthonormal basis for the Krylow subspace K_m = span{ r_0, A r_0, ..., A^{M-1} r_0}, where the initial residual r_0 = b - A x_0 is given from an initial solution guess x_0.
GMRES minimizes the square of the residual norm restricted to the subspace K_m.
So, the residual norm is given by
{b-Ax}
For a projection method we have x = x_0 + V_m y, so
b - A x = b - A(x_0 + V_m y) = r_0 - A V_M y
Use the ``factorization'' we get from the Arnoldi Method
b - A x = r_0 - V_{M+1} H_{M} y
Define beta as {r_0}, note that the first column, v_1, of V_{m+1} is a unit vector in the direction of r_0 (because the Arnoldi method is basically a modified Graham-Schmidt ortonormalization of K_M).
b - A x = beta v_1 - V_{M+1} H_{M} y = V_{M+1} ( beta e_1 - H_{M} y )
Since the columns of V_{M+1} are orthonormal we have
| r | = | b - A x | = | beta e_1 - H_{M} y |
It can be shown that the solution for the minimizer of the square of the residual norm is given by
y = R_M^-1 g_m
where R_M is an M M matrix obtained from transforming H_{M} into upper triangular form and deleting the last row. g_m is the resulting vector from applying the same transformations to the right-hand side beta e_1, and deleting the last entry.
The transformation to upper triangular form can be done iteratively in tandem with the Arnoldi method.
The idea is to apply Givens rotations and store these. Whenever H_M is expanded with a new row and column. The neat thing about this is that after having transformed the system into upper triangular form the last entry of the right hand side vector contains the value of the current residual iterate.
| j | The index of the last added column to the matrix H. The new row will have index j+1. | |
| H | The hessenberg matrix. | |
| g | The right hand side vector. | |
| c | A vector storing the value of cos(theta) of all the previously applied Given rotations. | |
| s | A vector storing the value of sin(theta) of all the previously applied Given rotations. |
| bool OpenTissue::math::big::detail::is_orthonormal | ( | size_type | m, | |
| std::vector< vector_type > const & | v | |||
| ) | [inline] |
Test if basis is orthonormal. This method is only used for internal testing purpose.
| m | The number of vectors in the basis for the preconditioned Krylov subspace: K_m = span{r_0, M^{-1} A r_0,...,( M^{-1} A )^{m-1} r_0, } . That is the indices of the basis vectors range from v_0 to v_{m-1}. | |
| v | The vectors for the basis of the Krylow subspace. |
| bool OpenTissue::math::big::detail::is_upper_triangular | ( | size_type | m, | |
| matrix_type const & | H | |||
| ) | [inline] |
Test if matrix is upper triangular. This method is only used for internal testing purpose.
| m | The number of vectors in the basis for the preconditioned Krylov subspace: K_m = span{r_0, M^{-1} A r_0,...,( M^{-1} A )^{m-1} r_0, } . That is the indices of the basis vectors range from v_0 to v_{m-1}. | |
| H | The matrix that should be tested. |
| void OpenTissue::math::big::detail::set_rotation | ( | value_type & | a, | |
| value_type & | b, | |||
| value_type const & | c, | |||
| value_type const & | s | |||
| ) | [inline] |
This method applies a rotation, theta, around the z-axis, to the vector [a,b]^T.
| a | Upon call this argument holds the intial value of the x-component of the vector to be rorated. Upon return this argument holds the result of the rotation. | |
| b | Upon call this argument holds the intial value of the y-component of the vector to be rorated. Upon return this argument holds the result of the rotation. | |
| c | This argument holds the value of cos(theta), where theta is the angle the vector needs to be rotated around the z-axis. | |
| s | This argument holds the value of sin(theta), where theta is the angle the vector needs to be rotated around the z-axis. |
| void OpenTissue::math::big::detail::update | ( | vector_type & | x, | |
| size_type | m, | |||
| matrix_type & | H, | |||
| vector_type const & | g, | |||
| std::vector< vector_type > const & | v | |||
| ) | [inline] |
Update the current iterate. This method computes
x = x_0 + V_M y
where
y = arg min { b - A x} y
| x | The current value of the iterate (intially x_0). Upon return this argument holds the updated value of the iterate x. | |
| m | The number of vectors in the basis for the preconditioned Krylov subspace: K_m = span{r_0, M^{-1} A r_0,...,( M^{-1} A )^{m-1} r_0, } . That is the indices of the basis vectors range from v_0 to v_{m-1}. | |
| H | The Hessenberg matrix transformed into triangular form. The triangular form consist of the m m submatrix of H. | |
| g | The right hand side vector for the least square sub problem transformed into triangular form. The transformed vector is the m-dimensional subpart of g. | |
| v | The vectors for the basis of the Krylow subspace. |
1.7.1