1 /* linalg/hermtd.c
2 *
3 * Copyright (C) 2001, 2007, 2009 Brian Gough
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 3 of the License, or (at
8 * your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18 */
19
20 /* Factorise a hermitian matrix A into
21 *
22 * A = U T U'
23 *
24 * where U is unitary and T is real symmetric tridiagonal. Only the
25 * diagonal and lower triangular part of A is referenced and modified.
26 *
27 * On exit, T is stored in the diagonal and first subdiagonal of
28 * A. Since T is symmetric the upper diagonal is not stored.
29 *
30 * U is stored as a packed set of Householder transformations in the
31 * lower triangular part of the input matrix below the first subdiagonal.
32 *
33 * The full matrix for U can be obtained as the product
34 *
35 * U = U_N ... U_2 U_1
36 *
37 * where
38 *
39 * U_i = (I - tau_i * v_i * v_i')
40 *
41 * and where v_i is a Householder vector
42 *
43 * v_i = [0, ..., 0, 1, A(i+2,i), A(i+3,i), ... , A(N,i)]
44 *
45 * This storage scheme is the same as in LAPACK. See LAPACK's
46 * chetd2.f for details.
47 *
48 * See Golub & Van Loan, "Matrix Computations" (3rd ed), Section 8.3 */
49
50 #include <config.h>
51 #include <stdlib.h>
52 #include <gsl/gsl_math.h>
53 #include <gsl/gsl_vector.h>
54 #include <gsl/gsl_matrix.h>
55 #include <gsl/gsl_blas.h>
56 #include <gsl/gsl_complex_math.h>
57
58 #include <gsl/gsl_linalg.h>
59
60 int
gsl_linalg_hermtd_decomp(gsl_matrix_complex * A,gsl_vector_complex * tau)61 gsl_linalg_hermtd_decomp (gsl_matrix_complex * A, gsl_vector_complex * tau)
62 {
63 if (A->size1 != A->size2)
64 {
65 GSL_ERROR ("hermitian tridiagonal decomposition requires square matrix",
66 GSL_ENOTSQR);
67 }
68 else if (tau->size + 1 != A->size1)
69 {
70 GSL_ERROR ("size of tau must be (matrix size - 1)", GSL_EBADLEN);
71 }
72 else
73 {
74 const size_t N = A->size1;
75 size_t i;
76
77 const gsl_complex zero = gsl_complex_rect (0.0, 0.0);
78 const gsl_complex one = gsl_complex_rect (1.0, 0.0);
79 const gsl_complex neg_one = gsl_complex_rect (-1.0, 0.0);
80
81 for (i = 0 ; i < N - 1; i++)
82 {
83 gsl_vector_complex_view c = gsl_matrix_complex_column (A, i);
84 gsl_vector_complex_view v = gsl_vector_complex_subvector (&c.vector, i + 1, N - (i + 1));
85 gsl_complex tau_i = gsl_linalg_complex_householder_transform (&v.vector);
86
87 /* Apply the transformation H^T A H to the remaining columns */
88
89 if ((i + 1) < (N - 1)
90 && !(GSL_REAL(tau_i) == 0.0 && GSL_IMAG(tau_i) == 0.0))
91 {
92 gsl_matrix_complex_view m =
93 gsl_matrix_complex_submatrix (A, i + 1, i + 1,
94 N - (i+1), N - (i+1));
95 gsl_complex ei = gsl_vector_complex_get(&v.vector, 0);
96 gsl_vector_complex_view x = gsl_vector_complex_subvector (tau, i, N-(i+1));
97 gsl_vector_complex_set (&v.vector, 0, one);
98
99 /* x = tau * A * v */
100 gsl_blas_zhemv (CblasLower, tau_i, &m.matrix, &v.vector, zero, &x.vector);
101
102 /* w = x - (1/2) tau * (x' * v) * v */
103 {
104 gsl_complex xv, txv, alpha;
105 gsl_blas_zdotc(&x.vector, &v.vector, &xv);
106 txv = gsl_complex_mul(tau_i, xv);
107 alpha = gsl_complex_mul_real(txv, -0.5);
108 gsl_blas_zaxpy(alpha, &v.vector, &x.vector);
109 }
110
111 /* apply the transformation A = A - v w' - w v' */
112 gsl_blas_zher2(CblasLower, neg_one, &v.vector, &x.vector, &m.matrix);
113
114 gsl_vector_complex_set (&v.vector, 0, ei);
115 }
116
117 gsl_vector_complex_set (tau, i, tau_i);
118 }
119
120 return GSL_SUCCESS;
121 }
122 }
123
124
125 /* Form the orthogonal matrix U from the packed QR matrix */
126
127 int
gsl_linalg_hermtd_unpack(const gsl_matrix_complex * A,const gsl_vector_complex * tau,gsl_matrix_complex * U,gsl_vector * diag,gsl_vector * sdiag)128 gsl_linalg_hermtd_unpack (const gsl_matrix_complex * A,
129 const gsl_vector_complex * tau,
130 gsl_matrix_complex * U,
131 gsl_vector * diag,
132 gsl_vector * sdiag)
133 {
134 if (A->size1 != A->size2)
135 {
136 GSL_ERROR ("matrix A must be sqaure", GSL_ENOTSQR);
137 }
138 else if (tau->size + 1 != A->size1)
139 {
140 GSL_ERROR ("size of tau must be (matrix size - 1)", GSL_EBADLEN);
141 }
142 else if (U->size1 != A->size1 || U->size2 != A->size1)
143 {
144 GSL_ERROR ("size of U must match size of A", GSL_EBADLEN);
145 }
146 else if (diag->size != A->size1)
147 {
148 GSL_ERROR ("size of diagonal must match size of A", GSL_EBADLEN);
149 }
150 else if (sdiag->size + 1 != A->size1)
151 {
152 GSL_ERROR ("size of subdiagonal must be (matrix size - 1)", GSL_EBADLEN);
153 }
154 else
155 {
156 const size_t N = A->size1;
157 gsl_vector_complex_const_view zd = gsl_matrix_complex_const_diagonal(A);
158 gsl_vector_complex_const_view zsd = gsl_matrix_complex_const_subdiagonal(A, 1);
159 gsl_vector_const_view d = gsl_vector_complex_const_real(&zd.vector);
160 gsl_vector_const_view sd = gsl_vector_complex_const_real(&zsd.vector);
161 gsl_vector_complex * work = gsl_vector_complex_alloc(N);
162 size_t i;
163
164 /* initialize U to the identity */
165 gsl_matrix_complex_set_identity (U);
166
167 for (i = N - 1; i-- > 0;)
168 {
169 gsl_complex ti = gsl_vector_complex_get (tau, i);
170 gsl_vector_complex_const_view h =
171 gsl_matrix_complex_const_subcolumn (A, i, i + 1, N - i - 1);
172 gsl_matrix_complex_view m =
173 gsl_matrix_complex_submatrix (U, i + 1, i + 1, N - i - 1, N - i - 1);
174 gsl_vector_complex_view w = gsl_vector_complex_subvector(work, 0, N - i - 1);
175
176 gsl_linalg_complex_householder_left (ti, &h.vector, &m.matrix, &w.vector);
177 }
178
179 /* copy diagonal into diag */
180 gsl_vector_memcpy(diag, &d.vector);
181
182 /* copy subdiagonal into sdiag */
183 gsl_vector_memcpy(sdiag, &sd.vector);
184
185 gsl_vector_complex_free(work);
186
187 return GSL_SUCCESS;
188 }
189 }
190
191 int
gsl_linalg_hermtd_unpack_T(const gsl_matrix_complex * A,gsl_vector * diag,gsl_vector * sdiag)192 gsl_linalg_hermtd_unpack_T (const gsl_matrix_complex * A,
193 gsl_vector * diag,
194 gsl_vector * sdiag)
195 {
196 if (A->size1 != A->size2)
197 {
198 GSL_ERROR ("matrix A must be sqaure", GSL_ENOTSQR);
199 }
200 else if (diag->size != A->size1)
201 {
202 GSL_ERROR ("size of diagonal must match size of A", GSL_EBADLEN);
203 }
204 else if (sdiag->size + 1 != A->size1)
205 {
206 GSL_ERROR ("size of subdiagonal must be (matrix size - 1)", GSL_EBADLEN);
207 }
208 else
209 {
210 gsl_vector_complex_const_view zd = gsl_matrix_complex_const_diagonal(A);
211 gsl_vector_complex_const_view zsd = gsl_matrix_complex_const_subdiagonal(A, 1);
212 gsl_vector_const_view d = gsl_vector_complex_const_real(&zd.vector);
213 gsl_vector_const_view sd = gsl_vector_complex_const_real(&zsd.vector);
214
215 /* copy diagonal into diag */
216 gsl_vector_memcpy(diag, &d.vector);
217
218 /* copy subdiagonal into sdiag */
219 gsl_vector_memcpy(sdiag, &sd.vector);
220
221 return GSL_SUCCESS;
222 }
223 }
224