1 /* linalg/qr_band.c
2 *
3 * Copyright (C) 2020 Patrick Alken
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 #include <config.h>
21 #include <stdlib.h>
22 #include <string.h>
23 #include <gsl/gsl_linalg.h>
24 #include <gsl/gsl_math.h>
25 #include <gsl/gsl_vector.h>
26 #include <gsl/gsl_matrix.h>
27 #include <gsl/gsl_blas.h>
28
29 /* Factorise a (p,q) banded M x N matrix A into
30 *
31 * A = Q R
32 *
33 * where Q is orthogonal (M x M) and R is upper triangular (M x N).
34 *
35 * Example with M = 7, N = 6, (p,q) = (2,1)
36 *
37 * A = [ A11 A12 0 0 0 0 ]
38 * [ A21 A22 A23 0 0 0 ]
39 * [ A31 A32 A33 A34 0 0 ]
40 * [ 0 A42 A43 A44 A45 0 ]
41 * [ 0 0 A53 A54 A55 A56 ]
42 * [ 0 0 0 A64 A65 A66 ]
43 * [ 0 0 0 0 A75 A76 ]
44 *
45 * AB has dimensions N-by-(2p + q + 1)
46 *
47 * INPUT: OUTPUT:
48 *
49 * AB = [ * * * A11 A21 A31 ] AB = [ * * * R11 V21 V31 ]
50 * [ * * A12 A22 A32 A42 ] [ * * R12 R22 V32 V42 ]
51 * [ * 0 A23 A33 A43 A53 ] [ * R13 R23 R33 V43 V53 ]
52 * [ 0 0 A34 A44 A54 A64 ] [ R14 R24 R34 R44 V54 V64 ]
53 * [ 0 0 A45 A55 A65 A75 ] [ R25 R35 R45 R55 V65 V75 ]
54 * [ 0 0 A56 A66 A76 * ] [ R36 R46 R56 R66 V76 * ]
55 * -p-- -q- -1- ---p--- ---p--- -q- -1- ---p---
56 *
57 * The full matrix for Q can be obtained as the product
58 *
59 * Q = Q_k .. Q_2 Q_1
60 *
61 * where k = MIN(M,N) and
62 *
63 * Q_i = (I - tau_i * v_i * v_i')
64 *
65 * and where v_i is a Householder vector
66 */
67
68 int
gsl_linalg_QR_band_decomp_L2(const size_t M,const size_t p,const size_t q,gsl_matrix * AB,gsl_vector * tau)69 gsl_linalg_QR_band_decomp_L2 (const size_t M, const size_t p, const size_t q, gsl_matrix * AB, gsl_vector * tau)
70 {
71 const size_t N = AB->size1;
72
73 if (tau->size != N)
74 {
75 GSL_ERROR ("tau must have length N", GSL_EBADLEN);
76 }
77 else if (AB->size2 != 2*p + q + 1)
78 {
79 GSL_ERROR ("dimensions of AB are inconsistent with (p,q)", GSL_EBADLEN);
80 }
81 else
82 {
83 const size_t minMN = GSL_MIN(M, N);
84 size_t j;
85
86 /* set AB(:,1:p) to zero */
87 if (p > 0)
88 {
89 gsl_matrix_view m = gsl_matrix_submatrix(AB, 0, 0, N, p);
90 gsl_matrix_set_zero(&m.matrix);
91 }
92
93 for (j = 0; j < minMN; ++j)
94 {
95 /* Compute the Householder transformation to reduce the j-th
96 column of the matrix to a multiple of the j-th unit vector */
97
98 size_t k1 = GSL_MIN(p + 1, M - j); /* number of non-zero elements of this column, including diagonal element */
99 size_t k2 = GSL_MIN(p + q, N - j - 1); /* number of columns to update */
100 gsl_vector_view c = gsl_matrix_subrow(AB, j, p + q, k1);
101 double tau_j = gsl_linalg_householder_transform (&(c.vector));
102 double * ptr = gsl_vector_ptr(&(c.vector), 0);
103
104 gsl_vector_set (tau, j, tau_j);
105
106 /* apply the transformation to the remaining columns */
107 if (k2 > 0)
108 {
109 gsl_matrix_view m = gsl_matrix_submatrix (AB, j + 1, p + q - 1, k2, k1);
110 gsl_vector_view work = gsl_vector_subvector(tau, j + 1, k2);
111 double tmp = *ptr;
112
113 m.matrix.tda -= 1; /* unskew matrix */
114
115 /* we want to compute H*A(j:j+k1-1,j+1:j+k2), but due to our using row-major order, the
116 * matrix m contains A(j:j+k1-1,j+1:j+k2)^T. So therefore we apply H from the right,
117 *
118 * [H*A]^T = A^T H^T = A^T H
119 */
120
121 *ptr = 1.0;
122 gsl_linalg_householder_right(tau_j, &(c.vector), &(m.matrix), &(work.vector));
123 *ptr = tmp;
124 }
125 }
126
127 return GSL_SUCCESS;
128 }
129 }
130
131 int
gsl_linalg_QR_band_unpack_L2(const size_t p,const size_t q,const gsl_matrix * QRB,const gsl_vector * tau,gsl_matrix * Q,gsl_matrix * R)132 gsl_linalg_QR_band_unpack_L2 (const size_t p, const size_t q, const gsl_matrix * QRB, const gsl_vector * tau,
133 gsl_matrix * Q, gsl_matrix * R)
134 {
135 const size_t M = Q->size1;
136 const size_t N = QRB->size1;
137
138 if (Q->size2 != M)
139 {
140 GSL_ERROR ("Q matrix must be square", GSL_ENOTSQR);
141 }
142 else if (R->size1 != M || R->size2 != N)
143 {
144 GSL_ERROR ("R matrix must be M x N", GSL_ENOTSQR);
145 }
146 else if (tau->size < GSL_MIN (M, N))
147 {
148 GSL_ERROR ("size of tau must be at least MIN(M,N)", GSL_EBADLEN);
149 }
150 else if (QRB->size2 != 2*p + q + 1)
151 {
152 GSL_ERROR ("dimensions of QRB are inconsistent with (p,q)", GSL_EBADLEN);
153 }
154 else
155 {
156 size_t i;
157
158 /* form matrix Q */
159 gsl_matrix_set_identity (Q);
160
161 for (i = GSL_MIN (M, N); i-- > 0;)
162 {
163 size_t k1 = GSL_MIN(p + 1, M - i); /* number of non-zero elements of this column, including diagonal element */
164 gsl_vector_const_view h = gsl_matrix_const_subrow(QRB, i, p + q, k1);
165 gsl_matrix_view m = gsl_matrix_submatrix (Q, i, i, k1, M - i);
166 double ti = gsl_vector_get (tau, i);
167 gsl_vector_view work = gsl_matrix_subcolumn(R, 0, 0, M - i);
168 double * ptr = gsl_vector_ptr((gsl_vector *) &h.vector, 0);
169 double tmp = *ptr;
170
171 *ptr = 1.0;
172 gsl_linalg_householder_left (ti, &h.vector, &m.matrix, &work.vector);
173 *ptr = tmp;
174 }
175
176 /* form matrix R */
177 gsl_matrix_set_zero(R);
178
179 for (i = 0; i <= GSL_MIN(p + q, N - 1); ++i)
180 {
181 gsl_vector_const_view src = gsl_matrix_const_subcolumn(QRB, p + q - i, i, GSL_MIN(M, N - i));
182 gsl_vector_view dest = gsl_matrix_superdiagonal(R, i);
183 gsl_vector_memcpy(&dest.vector, &src.vector);
184 }
185
186 return GSL_SUCCESS;
187 }
188 }
189