1 /* linalg/householder.c
2 *
3 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004, 2007, 2010 Gerard Jungman, 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 #include <config.h>
21 #include <stdlib.h>
22 #include <gsl/gsl_math.h>
23 #include <gsl/gsl_vector.h>
24 #include <gsl/gsl_matrix.h>
25 #include <gsl/gsl_blas.h>
26
27 #include <gsl/gsl_linalg.h>
28
29 /*
30 gsl_linalg_householder_transform()
31 Compute a householder transformation (tau,v) of a vector
32 x so that P x = [ I - tau*v*v' ] x annihilates x(1:n-1)
33
34 Inputs: v - on input, x vector
35 on output, householder vector v
36
37 Notes:
38 1) on output, v is normalized so that v[0] = 1. The 1 is
39 not actually stored; instead v[0] = -sign(x[0])*||x|| so
40 that:
41
42 P x = v[0] * e_1
43
44 Therefore external routines should take care when applying
45 the projection matrix P to vectors, taking into account
46 that v[0] should be 1 when doing so.
47 */
48
49 double
gsl_linalg_householder_transform(gsl_vector * v)50 gsl_linalg_householder_transform (gsl_vector * v)
51 {
52 /* replace v[0:n-1] with a householder vector (v[0:n-1]) and
53 coefficient tau that annihilate v[1:n-1] */
54
55 const size_t n = v->size ;
56
57 if (n == 1)
58 {
59 return 0.0; /* tau = 0 */
60 }
61 else
62 {
63 double alpha, beta, tau ;
64
65 gsl_vector_view x = gsl_vector_subvector (v, 1, n - 1) ;
66
67 double xnorm = gsl_blas_dnrm2 (&x.vector);
68
69 if (xnorm == 0)
70 {
71 return 0.0; /* tau = 0 */
72 }
73
74 alpha = gsl_vector_get (v, 0) ;
75 beta = - GSL_SIGN(alpha) * hypot(alpha, xnorm);
76 tau = (beta - alpha) / beta ;
77
78 {
79 double s = (alpha - beta);
80
81 if (fabs(s) > GSL_DBL_MIN)
82 {
83 gsl_blas_dscal (1.0 / s, &x.vector);
84 gsl_vector_set (v, 0, beta) ;
85 }
86 else
87 {
88 gsl_blas_dscal (GSL_DBL_EPSILON / s, &x.vector);
89 gsl_blas_dscal (1.0 / GSL_DBL_EPSILON, &x.vector);
90 gsl_vector_set (v, 0, beta) ;
91 }
92 }
93
94 return tau;
95 }
96 }
97
98 /*
99 gsl_linalg_householder_transform2()
100 Compute a householder transformation P so that
101
102 P [ alpha ] = [ beta ]
103 [ x(1:n-1) ] [ 0 ]
104
105 where
106
107 P = I - tau [ 1 ] [ 1 v' ]
108 [ v ]
109
110 Inputs: alpha - on input, alpha scalar
111 on output, beta scalar
112 v - length n vector
113 on input, v(1:n-1) contains x vector
114 on output, v(1:n-1) householder vector v
115 v(n) is not modified
116 */
117
118 double
gsl_linalg_householder_transform2(double * alpha,gsl_vector * v)119 gsl_linalg_householder_transform2 (double * alpha, gsl_vector * v)
120 {
121 const size_t n = v->size;
122
123 if (n == 1)
124 {
125 return 0.0; /* tau = 0 */
126 }
127 else
128 {
129 double beta, tau;
130 gsl_vector_view x = gsl_vector_subvector (v, 0, n - 1);
131 double xnorm = gsl_blas_dnrm2 (&x.vector);
132
133 if (xnorm == 0)
134 {
135 return 0.0; /* tau = 0 */
136 }
137
138 beta = - GSL_SIGN(*alpha) * hypot(*alpha, xnorm);
139 tau = (beta - *alpha) / beta;
140
141 {
142 double s = (*alpha - beta);
143
144 if (fabs(s) > GSL_DBL_MIN)
145 {
146 gsl_blas_dscal (1.0 / s, &x.vector);
147 *alpha = beta;
148 }
149 else
150 {
151 gsl_blas_dscal (GSL_DBL_EPSILON / s, &x.vector);
152 gsl_blas_dscal (1.0 / GSL_DBL_EPSILON, &x.vector);
153 *alpha = beta;
154 }
155 }
156
157 return tau;
158 }
159 }
160
161 int
gsl_linalg_householder_hm(double tau,const gsl_vector * v,gsl_matrix * A)162 gsl_linalg_householder_hm (double tau, const gsl_vector * v, gsl_matrix * A)
163 {
164 /* applies a householder transformation v,tau to matrix m */
165
166 if (tau == 0.0)
167 {
168 return GSL_SUCCESS;
169 }
170
171 #ifdef USE_BLAS
172 {
173 gsl_vector_const_view v1 = gsl_vector_const_subvector (v, 1, v->size - 1);
174 gsl_matrix_view A1 = gsl_matrix_submatrix (A, 1, 0, A->size1 - 1, A->size2);
175 size_t j;
176
177 for (j = 0; j < A->size2; j++)
178 {
179 double wj = 0.0;
180 gsl_vector_view A1j = gsl_matrix_column(&A1.matrix, j);
181 gsl_blas_ddot (&A1j.vector, &v1.vector, &wj);
182 wj += gsl_matrix_get(A,0,j);
183
184 {
185 double A0j = gsl_matrix_get (A, 0, j);
186 gsl_matrix_set (A, 0, j, A0j - tau * wj);
187 }
188
189 gsl_blas_daxpy (-tau * wj, &v1.vector, &A1j.vector);
190 }
191 }
192 #else
193 {
194 size_t i, j;
195
196 for (j = 0; j < A->size2; j++)
197 {
198 /* Compute wj = Akj vk */
199
200 double wj = gsl_matrix_get(A,0,j);
201
202 for (i = 1; i < A->size1; i++) /* note, computed for v(0) = 1 above */
203 {
204 wj += gsl_matrix_get(A,i,j) * gsl_vector_get(v,i);
205 }
206
207 /* Aij = Aij - tau vi wj */
208
209 /* i = 0 */
210 {
211 double A0j = gsl_matrix_get (A, 0, j);
212 gsl_matrix_set (A, 0, j, A0j - tau * wj);
213 }
214
215 /* i = 1 .. M-1 */
216
217 for (i = 1; i < A->size1; i++)
218 {
219 double Aij = gsl_matrix_get (A, i, j);
220 double vi = gsl_vector_get (v, i);
221 gsl_matrix_set (A, i, j, Aij - tau * vi * wj);
222 }
223 }
224 }
225 #endif
226
227 return GSL_SUCCESS;
228 }
229
230 int
gsl_linalg_householder_mh(double tau,const gsl_vector * v,gsl_matrix * A)231 gsl_linalg_householder_mh (double tau, const gsl_vector * v, gsl_matrix * A)
232 {
233 /* applies a householder transformation v,tau to matrix m from the
234 right hand side in order to zero out rows */
235
236 if (tau == 0)
237 return GSL_SUCCESS;
238
239 /* A = A - tau w v' */
240
241 #ifdef USE_BLAS
242 {
243 gsl_vector_const_view v1 = gsl_vector_const_subvector (v, 1, v->size - 1);
244 gsl_matrix_view A1 = gsl_matrix_submatrix (A, 0, 1, A->size1, A->size2-1);
245 size_t i;
246
247 for (i = 0; i < A->size1; i++)
248 {
249 double wi = 0.0;
250 gsl_vector_view A1i = gsl_matrix_row(&A1.matrix, i);
251 gsl_blas_ddot (&A1i.vector, &v1.vector, &wi);
252 wi += gsl_matrix_get(A,i,0);
253
254 {
255 double Ai0 = gsl_matrix_get (A, i, 0);
256 gsl_matrix_set (A, i, 0, Ai0 - tau * wi);
257 }
258
259 gsl_blas_daxpy(-tau * wi, &v1.vector, &A1i.vector);
260 }
261 }
262 #else
263 {
264 size_t i, j;
265
266 for (i = 0; i < A->size1; i++)
267 {
268 double wi = gsl_matrix_get(A,i,0);
269
270 for (j = 1; j < A->size2; j++) /* note, computed for v(0) = 1 above */
271 {
272 wi += gsl_matrix_get(A,i,j) * gsl_vector_get(v,j);
273 }
274
275 /* j = 0 */
276
277 {
278 double Ai0 = gsl_matrix_get (A, i, 0);
279 gsl_matrix_set (A, i, 0, Ai0 - tau * wi);
280 }
281
282 /* j = 1 .. N-1 */
283
284 for (j = 1; j < A->size2; j++)
285 {
286 double vj = gsl_vector_get (v, j);
287 double Aij = gsl_matrix_get (A, i, j);
288 gsl_matrix_set (A, i, j, Aij - tau * wi * vj);
289 }
290 }
291 }
292 #endif
293
294 return GSL_SUCCESS;
295 }
296
297 int
gsl_linalg_householder_hv(double tau,const gsl_vector * v,gsl_vector * w)298 gsl_linalg_householder_hv (double tau, const gsl_vector * v, gsl_vector * w)
299 {
300 /* applies a householder transformation v to vector w */
301 const size_t N = v->size;
302
303 if (tau == 0)
304 return GSL_SUCCESS ;
305
306 {
307 /* compute d = v'w */
308
309 double w0 = gsl_vector_get(w,0);
310 double d1, d;
311
312 gsl_vector_const_view v1 = gsl_vector_const_subvector(v, 1, N-1);
313 gsl_vector_view w1 = gsl_vector_subvector(w, 1, N-1);
314
315 /* compute d1 = v(2:n)'w(2:n) */
316 gsl_blas_ddot (&v1.vector, &w1.vector, &d1);
317
318 /* compute d = v'w = w(1) + d1 since v(1) = 1 */
319 d = w0 + d1;
320
321 /* compute w = w - tau (v) (v'w) */
322
323 gsl_vector_set (w, 0, w0 - tau * d);
324 gsl_blas_daxpy (-tau * d, &v1.vector, &w1.vector);
325 }
326
327 return GSL_SUCCESS;
328 }
329
330 /*
331 gsl_linalg_householder_left()
332 Apply a Householder reflector
333
334 H = I - tau v v^T
335
336 to a M-by-N matrix A from the left
337
338 Inputs: tau - Householder coefficient
339 v - Householder vector, length M
340 A - (input/output) M-by-N matrix on input; on output, H*A
341 work - workspace, length N
342
343 Notes:
344 1) This routine replaces gsl_linalg_householder_hm
345 */
346
347 int
gsl_linalg_householder_left(const double tau,const gsl_vector * v,gsl_matrix * A,gsl_vector * work)348 gsl_linalg_householder_left(const double tau, const gsl_vector * v, gsl_matrix * A, gsl_vector * work)
349 {
350 const size_t M = A->size1;
351 const size_t N = A->size2;
352
353 if (v->size != M)
354 {
355 GSL_ERROR ("matrix must match Householder vector dimensions", GSL_EBADLEN);
356 }
357 else if (work->size != N)
358 {
359 GSL_ERROR ("workspace must match matrix", GSL_EBADLEN);
360 }
361 else
362 {
363 /* quick return */
364 if (tau == 0.0)
365 return GSL_SUCCESS;
366
367 /* work := A^T v */
368 gsl_blas_dgemv(CblasTrans, 1.0, A, v, 0.0, work);
369
370 /* A := A - tau v work^T */
371 gsl_blas_dger(-tau, v, work, A);
372
373 return GSL_SUCCESS;
374 }
375 }
376
377 /*
378 gsl_linalg_householder_right()
379 Apply a Householder reflector
380
381 H = I - tau v v^T
382
383 to a M-by-N matrix A from the right
384
385 Inputs: tau - Householder coefficient
386 v - Householder vector, length N
387 A - (input/output) M-by-N matrix on input; on output, A*H
388 work - workspace, length M
389
390 Notes:
391 1) v(1) is modified but is restored on output
392
393 2) This routine replaces gsl_linalg_householder_mh
394 */
395
396 int
gsl_linalg_householder_right(const double tau,const gsl_vector * v,gsl_matrix * A,gsl_vector * work)397 gsl_linalg_householder_right(const double tau, const gsl_vector * v, gsl_matrix * A, gsl_vector * work)
398 {
399 const size_t M = A->size1;
400 const size_t N = A->size2;
401
402 if (v->size != N)
403 {
404 GSL_ERROR ("matrix must match Householder vector dimensions", GSL_EBADLEN);
405 }
406 else if (work->size != M)
407 {
408 GSL_ERROR ("workspace must match matrix", GSL_EBADLEN);
409 }
410 else
411 {
412 double v0;
413
414 /* quick return */
415 if (tau == 0.0)
416 return GSL_SUCCESS;
417
418 v0 = gsl_vector_get(v, 0);
419 v->data[0] = 1.0;
420
421 /* work := A v */
422 gsl_blas_dgemv(CblasNoTrans, 1.0, A, v, 0.0, work);
423
424 /* A := A - tau work v^T */
425 gsl_blas_dger(-tau, work, v, A);
426
427 v->data[0] = v0;
428
429 return GSL_SUCCESS;
430 }
431 }
432
433 int
gsl_linalg_householder_hm1(double tau,gsl_matrix * A)434 gsl_linalg_householder_hm1 (double tau, gsl_matrix * A)
435 {
436 /* applies a householder transformation v,tau to a matrix being
437 build up from the identity matrix, using the first column of A as
438 a householder vector */
439
440 if (tau == 0)
441 {
442 size_t i,j;
443
444 gsl_matrix_set (A, 0, 0, 1.0);
445
446 for (j = 1; j < A->size2; j++)
447 {
448 gsl_matrix_set (A, 0, j, 0.0);
449 }
450
451 for (i = 1; i < A->size1; i++)
452 {
453 gsl_matrix_set (A, i, 0, 0.0);
454 }
455
456 return GSL_SUCCESS;
457 }
458
459 /* w = A' v */
460
461 #ifdef USE_BLAS
462 {
463 gsl_matrix_view A1 = gsl_matrix_submatrix (A, 1, 0, A->size1 - 1, A->size2);
464 gsl_vector_view v1 = gsl_matrix_column (&A1.matrix, 0);
465 size_t j;
466
467 for (j = 1; j < A->size2; j++)
468 {
469 double wj = 0.0; /* A0j * v0 */
470
471 gsl_vector_view A1j = gsl_matrix_column(&A1.matrix, j);
472 gsl_blas_ddot (&A1j.vector, &v1.vector, &wj);
473
474 /* A = A - tau v w' */
475
476 gsl_matrix_set (A, 0, j, - tau * wj);
477
478 gsl_blas_daxpy(-tau*wj, &v1.vector, &A1j.vector);
479 }
480
481 gsl_blas_dscal(-tau, &v1.vector);
482
483 gsl_matrix_set (A, 0, 0, 1.0 - tau);
484 }
485 #else
486 {
487 size_t i, j;
488
489 for (j = 1; j < A->size2; j++)
490 {
491 double wj = 0.0; /* A0j * v0 */
492
493 for (i = 1; i < A->size1; i++)
494 {
495 double vi = gsl_matrix_get(A, i, 0);
496 wj += gsl_matrix_get(A,i,j) * vi;
497 }
498
499 /* A = A - tau v w' */
500
501 gsl_matrix_set (A, 0, j, - tau * wj);
502
503 for (i = 1; i < A->size1; i++)
504 {
505 double vi = gsl_matrix_get (A, i, 0);
506 double Aij = gsl_matrix_get (A, i, j);
507 gsl_matrix_set (A, i, j, Aij - tau * vi * wj);
508 }
509 }
510
511 for (i = 1; i < A->size1; i++)
512 {
513 double vi = gsl_matrix_get(A, i, 0);
514 gsl_matrix_set(A, i, 0, -tau * vi);
515 }
516
517 gsl_matrix_set (A, 0, 0, 1.0 - tau);
518 }
519 #endif
520
521 return GSL_SUCCESS;
522 }
523