1 /* linalg/cod.c
2 *
3 * Copyright (C) 2016, 2017 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_math.h>
24 #include <gsl/gsl_vector.h>
25 #include <gsl/gsl_matrix.h>
26 #include <gsl/gsl_permute_vector.h>
27 #include <gsl/gsl_blas.h>
28 #include <gsl/gsl_linalg.h>
29
30 /*
31 * This module contains routines for factoring an M-by-N matrix A as:
32 *
33 * A P = Q R Z^T
34 *
35 * known as the Complete Orthogonal Decomposition, where:
36 *
37 * P is a N-by-N permutation matrix
38 * Q is M-by-M orthogonal
39 * R has an r-by-r upper triangular block
40 * Z is N-by-N orthogonal
41 *
42 * When A is full rank, Z = I and this becomes the QR decomposition
43 * with column pivoting. When A is rank deficient, then
44 *
45 * R = [ R11 0 ] where R11 is r-by-r and r = rank(A)
46 * [ 0 0 ]
47 */
48
49 static int cod_RZ(gsl_matrix * A, gsl_vector * tau);
50 static double cod_householder_transform(double *alpha, gsl_vector * v);
51 static int cod_householder_mh(const double tau, const gsl_vector * v,
52 gsl_matrix * A, gsl_vector * work);
53 static int cod_householder_hv(const double tau, const gsl_vector * v, gsl_vector * w);
54 static int cod_householder_Zvec(const gsl_matrix * QRZT, const gsl_vector * tau_Z, const size_t rank,
55 gsl_vector * v);
56 static int cod_trireg_solve(const gsl_matrix * R, const double lambda, const gsl_vector * b,
57 gsl_matrix * S, gsl_vector * x, gsl_vector * work);
58
59 int
gsl_linalg_COD_decomp_e(gsl_matrix * A,gsl_vector * tau_Q,gsl_vector * tau_Z,gsl_permutation * p,double tol,size_t * rank,gsl_vector * work)60 gsl_linalg_COD_decomp_e(gsl_matrix * A, gsl_vector * tau_Q, gsl_vector * tau_Z,
61 gsl_permutation * p, double tol, size_t * rank, gsl_vector * work)
62 {
63 const size_t M = A->size1;
64 const size_t N = A->size2;
65
66 if (tau_Q->size != GSL_MIN (M, N))
67 {
68 GSL_ERROR ("size of tau_Q must be MIN(M,N)", GSL_EBADLEN);
69 }
70 else if (tau_Z->size != GSL_MIN (M, N))
71 {
72 GSL_ERROR ("size of tau_Z must be MIN(M,N)", GSL_EBADLEN);
73 }
74 else if (p->size != N)
75 {
76 GSL_ERROR ("permutation size must be N", GSL_EBADLEN);
77 }
78 else if (work->size != N)
79 {
80 GSL_ERROR ("work size must be N", GSL_EBADLEN);
81 }
82 else
83 {
84 int status, signum;
85 size_t r;
86
87 /* decompose: A P = Q R */
88 status = gsl_linalg_QRPT_decomp(A, tau_Q, p, &signum, work);
89 if (status)
90 return status;
91
92 /* estimate rank of A */
93 r = gsl_linalg_QRPT_rank(A, tol);
94
95 if (r < N)
96 {
97 /*
98 * matrix is rank-deficient, so that the R factor is
99 *
100 * R = [ R11 R12 ] =~ [ R11 R12 ]
101 * [ 0 R22 ] [ 0 0 ]
102 *
103 * compute RZ decomposition of upper trapezoidal matrix
104 * [ R11 R12 ] = [ R11~ 0 ] Z
105 */
106 gsl_matrix_view R_upper = gsl_matrix_submatrix(A, 0, 0, r, N);
107 gsl_vector_view t = gsl_vector_subvector(tau_Z, 0, r);
108
109 cod_RZ(&R_upper.matrix, &t.vector);
110 }
111
112 *rank = r;
113
114 return GSL_SUCCESS;
115 }
116 }
117
118 int
gsl_linalg_COD_decomp(gsl_matrix * A,gsl_vector * tau_Q,gsl_vector * tau_Z,gsl_permutation * p,size_t * rank,gsl_vector * work)119 gsl_linalg_COD_decomp(gsl_matrix * A, gsl_vector * tau_Q, gsl_vector * tau_Z,
120 gsl_permutation * p, size_t * rank, gsl_vector * work)
121 {
122 return gsl_linalg_COD_decomp_e(A, tau_Q, tau_Z, p, -1.0, rank, work);
123 }
124
125 /*
126 gsl_linalg_COD_lssolve()
127 Find the least squares solution to the overdetermined system
128
129 min ||b - A x||^2
130
131 for M >= N using the COD factorization A P = Q R Z
132
133 Inputs: QRZT - matrix A, in COD compressed format, M-by-N
134 tau_Q - Householder scalars for Q, length min(M,N)
135 tau_Z - Householder scalars for Z, length min(M,N)
136 perm - permutation matrix
137 rank - rank of A
138 b - rhs vector, length M
139 x - (output) solution vector, length N
140 residual - (output) residual vector, b - A x, length M
141 */
142
143 int
gsl_linalg_COD_lssolve(const gsl_matrix * QRZT,const gsl_vector * tau_Q,const gsl_vector * tau_Z,const gsl_permutation * perm,const size_t rank,const gsl_vector * b,gsl_vector * x,gsl_vector * residual)144 gsl_linalg_COD_lssolve (const gsl_matrix * QRZT, const gsl_vector * tau_Q, const gsl_vector * tau_Z,
145 const gsl_permutation * perm, const size_t rank, const gsl_vector * b,
146 gsl_vector * x, gsl_vector * residual)
147 {
148 const size_t M = QRZT->size1;
149 const size_t N = QRZT->size2;
150
151 if (M < N)
152 {
153 GSL_ERROR ("QRZT matrix must have M>=N", GSL_EBADLEN);
154 }
155 else if (M != b->size)
156 {
157 GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
158 }
159 else if (rank > GSL_MIN (M, N))
160 {
161 GSL_ERROR ("rank must be <= MIN(M,N)", GSL_EBADLEN);
162 }
163 else if (N != x->size)
164 {
165 GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
166 }
167 else if (M != residual->size)
168 {
169 GSL_ERROR ("matrix size must match residual size", GSL_EBADLEN);
170 }
171 else
172 {
173 gsl_matrix_const_view R11 = gsl_matrix_const_submatrix (QRZT, 0, 0, rank, rank);
174 gsl_vector_view QTb1 = gsl_vector_subvector(residual, 0, rank);
175 gsl_vector_view x1 = gsl_vector_subvector(x, 0, rank);
176
177 gsl_vector_set_zero(x);
178
179 /* compute residual = Q^T b = [ c1 ; c2 ] */
180 gsl_vector_memcpy(residual, b);
181 gsl_linalg_QR_QTvec (QRZT, tau_Q, residual);
182
183 /* solve x1 := R11^{-1} (Q^T b)(1:r) */
184 gsl_vector_memcpy(&(x1.vector), &(QTb1.vector));
185 gsl_blas_dtrsv(CblasUpper, CblasNoTrans, CblasNonUnit, &(R11.matrix), &(x1.vector));
186
187 /* compute Z ( R11^{-1} x1; 0 ) */
188 cod_householder_Zvec(QRZT, tau_Z, rank, x);
189
190 /* compute x = P Z^T ( R11^{-1} x1; 0 ) */
191 gsl_permute_vector_inverse(perm, x);
192
193 /* compute residual = b - A x = Q (Q^T b - R [ R11^{-1} x1; 0 ]) = Q [ 0 ; c2 ] */
194 gsl_vector_set_zero(&(QTb1.vector));
195 gsl_linalg_QR_Qvec(QRZT, tau_Q, residual);
196
197 return GSL_SUCCESS;
198 }
199 }
200
201 /*
202 gsl_linalg_COD_lssolve2()
203 Find the least squares solution to the Tikhonov regularized
204 system in standard form:
205
206 min ||b - A x||^2 + lambda^2 ||x||^2
207
208 for M >= N using the COD factorization A P = Q R Z
209
210 Inputs: lambda - parameter
211 QRZT - matrix A, in COD compressed format, M-by-N
212 tau_Q - Householder scalars for Q, length min(M,N)
213 tau_Z - Householder scalars for Z, length min(M,N)
214 perm - permutation matrix
215 rank - rank of A
216 b - rhs vector, length M
217 x - (output) solution vector, length N
218 residual - (output) residual vector, b - A x, length M
219 S - workspace, rank-by-rank
220 work - workspace, length rank
221 */
222
223 int
gsl_linalg_COD_lssolve2(const double lambda,const gsl_matrix * QRZT,const gsl_vector * tau_Q,const gsl_vector * tau_Z,const gsl_permutation * perm,const size_t rank,const gsl_vector * b,gsl_vector * x,gsl_vector * residual,gsl_matrix * S,gsl_vector * work)224 gsl_linalg_COD_lssolve2 (const double lambda, const gsl_matrix * QRZT, const gsl_vector * tau_Q, const gsl_vector * tau_Z,
225 const gsl_permutation * perm, const size_t rank, const gsl_vector * b,
226 gsl_vector * x, gsl_vector * residual, gsl_matrix * S, gsl_vector * work)
227 {
228 const size_t M = QRZT->size1;
229 const size_t N = QRZT->size2;
230
231 if (M < N)
232 {
233 GSL_ERROR ("QRZT matrix must have M>=N", GSL_EBADLEN);
234 }
235 else if (M != b->size)
236 {
237 GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
238 }
239 else if (rank > GSL_MIN (M, N))
240 {
241 GSL_ERROR ("rank must be <= MIN(M,N)", GSL_EBADLEN);
242 }
243 else if (N != x->size)
244 {
245 GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
246 }
247 else if (M != residual->size)
248 {
249 GSL_ERROR ("matrix size must match residual size", GSL_EBADLEN);
250 }
251 else if (S->size1 != rank || S->size2 != rank)
252 {
253 GSL_ERROR ("S must be rank-by-rank", GSL_EBADLEN);
254 }
255 else if (work->size != rank)
256 {
257 GSL_ERROR ("work must be length rank", GSL_EBADLEN);
258 }
259 else
260 {
261 gsl_matrix_const_view R11 = gsl_matrix_const_submatrix (QRZT, 0, 0, rank, rank);
262 gsl_vector_view c1 = gsl_vector_subvector(residual, 0, rank);
263 gsl_vector_view y1 = gsl_vector_subvector(x, 0, rank);
264
265 gsl_vector_set_zero(x);
266
267 /* compute residual = Q^T b = [ c1 ; c2 ]*/
268 gsl_vector_memcpy(residual, b);
269 gsl_linalg_QR_QTvec (QRZT, tau_Q, residual);
270
271 /* solve [ R11 ; lambda*I ] y1 = [ (Q^T b)(1:r) ; 0 ] */
272 cod_trireg_solve(&(R11.matrix), lambda, &(c1.vector), S, &(y1.vector), work);
273
274 /* save y1 for later residual calculation */
275 gsl_vector_memcpy(work, &(y1.vector));
276
277 /* compute Z [ y1; 0 ] */
278 cod_householder_Zvec(QRZT, tau_Z, rank, x);
279
280 /* compute x = P Z^T ( y1; 0 ) */
281 gsl_permute_vector_inverse(perm, x);
282
283 /* compute residual = b - A x = Q (Q^T b - [ R11 y1; 0 ]) = Q [ c1 - R11*y1 ; c2 ] */
284
285 /* work = R11*y1 */
286 gsl_blas_dtrmv(CblasUpper, CblasNoTrans, CblasNonUnit, &(R11.matrix), work);
287
288 gsl_vector_sub(&(c1.vector), work);
289 gsl_linalg_QR_Qvec(QRZT, tau_Q, residual);
290
291 return GSL_SUCCESS;
292 }
293 }
294
295 /*
296 gsl_linalg_COD_unpack()
297 Unpack encoded COD decomposition into the matrices Q,R,Z,P
298
299 Inputs: QRZT - encoded COD decomposition
300 tau_Q - Householder scalars for Q
301 tau_Z - Householder scalars for Z
302 rank - rank of matrix (as determined from gsl_linalg_COD_decomp)
303 Q - (output) M-by-M matrix Q
304 R - (output) M-by-N matrix R
305 Z - (output) N-by-N matrix Z
306 */
307
308 int
gsl_linalg_COD_unpack(const gsl_matrix * QRZT,const gsl_vector * tau_Q,const gsl_vector * tau_Z,const size_t rank,gsl_matrix * Q,gsl_matrix * R,gsl_matrix * Z)309 gsl_linalg_COD_unpack(const gsl_matrix * QRZT, const gsl_vector * tau_Q,
310 const gsl_vector * tau_Z, const size_t rank, gsl_matrix * Q,
311 gsl_matrix * R, gsl_matrix * Z)
312 {
313 const size_t M = QRZT->size1;
314 const size_t N = QRZT->size2;
315
316 if (tau_Q->size != GSL_MIN (M, N))
317 {
318 GSL_ERROR ("size of tau_Q must be MIN(M,N)", GSL_EBADLEN);
319 }
320 else if (tau_Z->size != GSL_MIN (M, N))
321 {
322 GSL_ERROR ("size of tau_Z must be MIN(M,N)", GSL_EBADLEN);
323 }
324 else if (rank > GSL_MIN (M, N))
325 {
326 GSL_ERROR ("rank must be <= MIN(M,N)", GSL_EBADLEN);
327 }
328 else if (Q->size1 != M || Q->size2 != M)
329 {
330 GSL_ERROR ("Q must by M-by-M", GSL_EBADLEN);
331 }
332 else if (R->size1 != M || R->size2 != N)
333 {
334 GSL_ERROR ("R must by M-by-N", GSL_EBADLEN);
335 }
336 else if (Z->size1 != N || Z->size2 != N)
337 {
338 GSL_ERROR ("Z must by N-by-N", GSL_EBADLEN);
339 }
340 else
341 {
342 size_t i;
343 gsl_matrix_view R11 = gsl_matrix_submatrix(R, 0, 0, rank, rank);
344 gsl_matrix_const_view QRZT11 = gsl_matrix_const_submatrix(QRZT, 0, 0, rank, rank);
345
346 /* form Q matrix */
347
348 gsl_matrix_set_identity(Q);
349
350 for (i = GSL_MIN (M, N); i-- > 0;)
351 {
352 gsl_vector_const_view h = gsl_matrix_const_subcolumn (QRZT, i, i, M - i);
353 gsl_matrix_view m = gsl_matrix_submatrix (Q, i, i, M - i, M - i);
354 gsl_vector_view work = gsl_matrix_subcolumn (R, 0, 0, M - i);
355 double ti = gsl_vector_get (tau_Q, i);
356 double * ptr = gsl_vector_ptr((gsl_vector *) &h.vector, 0);
357 double tmp = *ptr;
358
359 *ptr = 1.0;
360 gsl_linalg_householder_left (ti, &h.vector, &m.matrix, &work.vector);
361 *ptr = tmp;
362 }
363
364 /* form Z matrix */
365 gsl_matrix_set_identity(Z);
366
367 if (rank < N)
368 {
369 gsl_vector_view work = gsl_matrix_row(R, 0); /* temporary workspace, size N */
370
371 /* multiply I by Z from the right */
372 gsl_linalg_COD_matZ(QRZT, tau_Z, rank, Z, &work.vector);
373 }
374
375 /* copy rank-by-rank upper triangle of QRZT into R and zero the rest */
376 gsl_matrix_set_zero(R);
377 gsl_matrix_tricpy(CblasUpper, CblasNonUnit, &R11.matrix, &QRZT11.matrix);
378
379 return GSL_SUCCESS;
380 }
381 }
382
383 /*
384 gsl_linalg_COD_matZ
385 Multiply an M-by-N matrix A on the right by Z (N-by-N)
386
387 Inputs: QRZT - encoded COD matrix
388 tau_Z - Householder scalars for Z
389 rank - matrix rank
390 A - on input, M-by-N matrix
391 on output, A * Z
392 work - workspace of length M
393 */
394
395 int
gsl_linalg_COD_matZ(const gsl_matrix * QRZT,const gsl_vector * tau_Z,const size_t rank,gsl_matrix * A,gsl_vector * work)396 gsl_linalg_COD_matZ(const gsl_matrix * QRZT, const gsl_vector * tau_Z, const size_t rank,
397 gsl_matrix * A, gsl_vector * work)
398 {
399 const size_t M = A->size1;
400 const size_t N = A->size2;
401
402 if (tau_Z->size != GSL_MIN (QRZT->size1, QRZT->size2))
403 {
404 GSL_ERROR("tau_Z must be GSL_MIN(M,N)", GSL_EBADLEN);
405 }
406 else if (QRZT->size2 != N)
407 {
408 GSL_ERROR("QRZT must have N columns", GSL_EBADLEN);
409 }
410 else if (work->size != M)
411 {
412 GSL_ERROR("workspace must be length M", GSL_EBADLEN);
413 }
414 else
415 {
416 /* if rank == N, then Z = I and there is nothing to do */
417 if (rank < N)
418 {
419 size_t i;
420
421 for (i = rank; i > 0 && i--; )
422 {
423 gsl_vector_const_view h = gsl_matrix_const_subrow (QRZT, i, rank, N - rank);
424 gsl_matrix_view m = gsl_matrix_submatrix (A, 0, i, M, N - i);
425 double ti = gsl_vector_get (tau_Z, i);
426 cod_householder_mh (ti, &h.vector, &m.matrix, work);
427 }
428 }
429
430 return GSL_SUCCESS;
431 }
432 }
433
434
435 /*********************************************
436 * INTERNAL ROUTINES *
437 *********************************************/
438
439 /*
440 cod_RZ()
441 Perform RZ decomposition of an upper trapezoidal matrix,
442
443 A = [ A11 A12 ] = [ R 0 ] Z
444
445 where A is M-by-N with N >= M, A11 is M-by-M upper triangular,
446 and A12 is M-by-(N-M). On output, Z is stored as Householder
447 reflectors in the A12 portion of A,
448
449 Z = Z(1) Z(2) ... Z(M)
450
451 Inputs: A - M-by-N matrix with N >= M
452 On input, upper trapezoidal matrix [ A11 A12 ]
453 On output, A11 is overwritten by R (subdiagonal elements
454 are not touched), and A12 is overwritten by Z in packed storage
455 tau - (output) Householder scalars, size M
456 */
457
458 static int
cod_RZ(gsl_matrix * A,gsl_vector * tau)459 cod_RZ(gsl_matrix * A, gsl_vector * tau)
460 {
461 const size_t M = A->size1;
462 const size_t N = A->size2;
463
464 if (tau->size != M)
465 {
466 GSL_ERROR("tau has wrong size", GSL_EBADLEN);
467 }
468 else if (N < M)
469 {
470 GSL_ERROR("N must be >= M", GSL_EINVAL);
471 }
472 else if (M == N)
473 {
474 /* quick return */
475 gsl_vector_set_all(tau, 0.0);
476 return GSL_SUCCESS;
477 }
478 else
479 {
480 size_t k;
481
482 for (k = M; k > 0 && k--; )
483 {
484 double *alpha = gsl_matrix_ptr(A, k, k);
485 gsl_vector_view z = gsl_matrix_subrow(A, k, M, N - M);
486 double tauk;
487
488 /* compute Householder reflection to zero [ A(k,k) A(k,M+1:N) ] */
489 tauk = cod_householder_transform(alpha, &z.vector);
490 gsl_vector_set(tau, k, tauk);
491
492 if ((tauk != 0) && (k > 0))
493 {
494 gsl_vector_view w = gsl_vector_subvector(tau, 0, k);
495 gsl_matrix_view B = gsl_matrix_submatrix(A, 0, k, k, N - k);
496
497 cod_householder_mh(tauk, &z.vector, &B.matrix, &w.vector);
498 }
499 }
500
501 return GSL_SUCCESS;
502 }
503 }
504
505 static double
cod_householder_transform(double * alpha,gsl_vector * v)506 cod_householder_transform(double *alpha, gsl_vector * v)
507 {
508 double beta, tau;
509 double xnorm = gsl_blas_dnrm2(v);
510
511 if (xnorm == 0)
512 {
513 return 0.0; /* tau = 0 */
514 }
515
516 beta = - (*alpha >= 0.0 ? +1.0 : -1.0) * gsl_hypot(*alpha, xnorm);
517 tau = (beta - *alpha) / beta;
518
519 {
520 double s = (*alpha - beta);
521
522 if (fabs(s) > GSL_DBL_MIN)
523 {
524 gsl_blas_dscal (1.0 / s, v);
525 }
526 else
527 {
528 gsl_blas_dscal (GSL_DBL_EPSILON / s, v);
529 gsl_blas_dscal (1.0 / GSL_DBL_EPSILON, v);
530 }
531
532 *alpha = beta;
533 }
534
535 return tau;
536 }
537
538 /*
539 cod_householder_hv
540 Apply Householder reflection H = (I - tau*v*v') to vector v from the left,
541
542 w' = H * w
543
544 Inputs: tau - Householder scalar
545 v - Householder vector, size M
546 w - on input, w vector, size M
547 on output, H * w
548
549 Notes:
550 1) Based on LAPACK routine DLARZ
551 */
552
553 static int
cod_householder_hv(const double tau,const gsl_vector * v,gsl_vector * w)554 cod_householder_hv(const double tau, const gsl_vector * v, gsl_vector * w)
555 {
556 if (tau == 0)
557 {
558 return GSL_SUCCESS; /* H = I */
559 }
560 else
561 {
562 const size_t M = w->size;
563 const size_t L = v->size;
564 double w0 = gsl_vector_get(w, 0);
565 gsl_vector_view w1 = gsl_vector_subvector(w, M - L, L);
566 double d1, d;
567
568 /* d1 := v . w(M-L:M) */
569 gsl_blas_ddot(v, &w1.vector, &d1);
570
571 /* d := w(1) + v . w(M-L:M) */
572 d = w0 + d1;
573
574 /* w(1) = w(1) - tau * d */
575 gsl_vector_set(w, 0, w0 - tau * d);
576
577 /* w(M-L:M) = w(M-L:M) - tau * d * v */
578 gsl_blas_daxpy(-tau * d, v, &w1.vector);
579
580 return GSL_SUCCESS;
581 }
582 }
583
584 /*
585 cod_householder_mh
586 Apply Householder reflection H = (I - tau*v*v') to matrix A from the right
587
588 Inputs: tau - Householder scalar
589 v - Householder vector, size N-M
590 A - matrix, size M-by-N
591 work - workspace, size M
592
593 Notes:
594 1) Based on LAPACK routine DLARZ
595 */
596
597 static int
cod_householder_mh(const double tau,const gsl_vector * v,gsl_matrix * A,gsl_vector * work)598 cod_householder_mh(const double tau, const gsl_vector * v, gsl_matrix * A,
599 gsl_vector * work)
600 {
601 if (tau == 0)
602 {
603 return GSL_SUCCESS; /* H = I */
604 }
605 else
606 {
607 const size_t M = A->size1;
608 const size_t N = A->size2;
609 const size_t L = v->size;
610 gsl_vector_view A1 = gsl_matrix_subcolumn(A, 0, 0, M);
611 gsl_matrix_view C = gsl_matrix_submatrix(A, 0, N - L, M, L);
612
613 /* work(1:M) = A(1:M,1) */
614 gsl_vector_memcpy(work, &A1.vector);
615
616 /* work(1:M) = work(1:M) + A(1:M,M+1:N) * v(1:N-M) */
617 gsl_blas_dgemv(CblasNoTrans, 1.0, &C.matrix, v, 1.0, work);
618
619 /* A(1:M,1) = A(1:M,1) - tau * work(1:M) */
620 gsl_blas_daxpy(-tau, work, &A1.vector);
621
622 /* A(1:M,M+1:N) = A(1:M,M+1:N) - tau * work(1:M) * v(1:N-M)' */
623 gsl_blas_dger(-tau, work, v, &C.matrix);
624
625 return GSL_SUCCESS;
626 }
627 }
628
629 /*
630 cod_householder_Zvec
631 Multiply a vector by Z
632
633 Inputs: QRZT - encoded COD matrix
634 tau_Z - Householder scalars for Z
635 rank - matrix rank
636 v - on input, vector of length N
637 on output, Z * v
638 */
639
640 static int
cod_householder_Zvec(const gsl_matrix * QRZT,const gsl_vector * tau_Z,const size_t rank,gsl_vector * v)641 cod_householder_Zvec(const gsl_matrix * QRZT, const gsl_vector * tau_Z, const size_t rank,
642 gsl_vector * v)
643 {
644 const size_t M = QRZT->size1;
645 const size_t N = QRZT->size2;
646
647 if (tau_Z->size != GSL_MIN (M, N))
648 {
649 GSL_ERROR("tau_Z must be GSL_MIN(M,N)", GSL_EBADLEN);
650 }
651 else if (v->size != N)
652 {
653 GSL_ERROR("v must be length N", GSL_EBADLEN);
654 }
655 else
656 {
657 if (rank < N)
658 {
659 size_t i;
660
661 for (i = 0; i < rank; ++i)
662 {
663 gsl_vector_const_view h = gsl_matrix_const_subrow (QRZT, i, rank, N - rank);
664 gsl_vector_view w = gsl_vector_subvector (v, i, N - i);
665 double ti = gsl_vector_get (tau_Z, i);
666 cod_householder_hv(ti, &h.vector, &w.vector);
667 }
668 }
669
670 return GSL_SUCCESS;
671 }
672 }
673
674 /*
675 cod_trireg_solve()
676
677 This function computes the solution to the least squares system
678
679 [ R ] x = [ b ]
680 [ lambda*I ] [ 0 ]
681
682 where R is an N-by-N upper triangular matrix, lambda is a scalar parameter,
683 and b is a vector of length N. This is done by computing the QR factorization
684
685 [ R ] = W S^T
686 [ lambda*I ]
687
688 where S^T is upper triangular, and solving
689
690 S^T x = W^T [ b ]
691 [ 0 ]
692
693 Inputs: R - full rank upper triangular matrix; the diagonal
694 elements are modified but restored on output
695 lambda - scalar parameter lambda
696 b - right hand side vector b
697 S - workspace, N-by-N
698 x - (output) least squares solution of the system
699 work - workspace of length N
700 */
701
702 static int
cod_trireg_solve(const gsl_matrix * R,const double lambda,const gsl_vector * b,gsl_matrix * S,gsl_vector * x,gsl_vector * work)703 cod_trireg_solve (const gsl_matrix * R, const double lambda, const gsl_vector * b,
704 gsl_matrix * S, gsl_vector * x, gsl_vector * work)
705 {
706 const size_t N = R->size2;
707 gsl_vector_const_view diag = gsl_matrix_const_diagonal(R);
708 size_t i, j, k;
709
710 if (lambda <= 0.0)
711 {
712 GSL_ERROR("lambda must be positive", GSL_EINVAL);
713 }
714
715 /* copy R and b to preserve input and initialise S; store diag(R) in work */
716 gsl_matrix_transpose_tricpy(CblasUpper, CblasUnit, S, R);
717 gsl_vector_memcpy(work, &diag.vector);
718 gsl_vector_memcpy(x, b);
719
720 /* eliminate the diagonal matrix lambda*I using Givens rotations */
721
722 for (j = 0; j < N; j++)
723 {
724 double bj = 0.0;
725
726 gsl_matrix_set (S, j, j, lambda);
727
728 for (k = j + 1; k < N; k++)
729 {
730 gsl_matrix_set (S, k, k, 0.0);
731 }
732
733 /* the transformations to eliminate the row of lambda*I modify only a
734 single element of b beyond the first n, which is initially
735 zero */
736
737 for (k = j; k < N; k++)
738 {
739 /* determine a Givens rotation which eliminates the
740 appropriate element in the current row of lambda*I */
741
742 double sine, cosine;
743
744 double xk = gsl_vector_get (x, k);
745 double rkk = gsl_vector_get (work, k);
746 double skk = gsl_matrix_get (S, k, k);
747
748 if (skk == 0)
749 {
750 continue;
751 }
752
753 if (fabs (rkk) < fabs (skk))
754 {
755 double cotangent = rkk / skk;
756 sine = 0.5 / sqrt (0.25 + 0.25 * cotangent * cotangent);
757 cosine = sine * cotangent;
758 }
759 else
760 {
761 double tangent = skk / rkk;
762 cosine = 0.5 / sqrt (0.25 + 0.25 * tangent * tangent);
763 sine = cosine * tangent;
764 }
765
766 /* Compute the modified diagonal element of r and the
767 modified element of [b,0] */
768
769 {
770 double new_rkk = cosine * rkk + sine * skk;
771 double new_xk = cosine * xk + sine * bj;
772
773 bj = -sine * xk + cosine * bj;
774
775 gsl_vector_set(work, k, new_rkk);
776 gsl_matrix_set(S, k, k, new_rkk);
777 gsl_vector_set(x, k, new_xk);
778 }
779
780 /* Accumulate the transformation in the row of s */
781
782 for (i = k + 1; i < N; i++)
783 {
784 double sik = gsl_matrix_get (S, i, k);
785 double sii = gsl_matrix_get (S, i, i);
786
787 double new_sik = cosine * sik + sine * sii;
788 double new_sii = -sine * sik + cosine * sii;
789
790 gsl_matrix_set(S, i, k, new_sik);
791 gsl_matrix_set(S, i, i, new_sii);
792 }
793 }
794 }
795
796 /* solve: S^T x = rhs in place */
797 gsl_blas_dtrsv(CblasLower, CblasTrans, CblasNonUnit, S, x);
798
799 return GSL_SUCCESS;
800 }
801