svdstep_base.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2008-2020, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 /* linalg/svdstep.c
24  *
25  * Copyright (C) 2007 Brian Gough
26  *
27  * This program is free software; you can redistribute it and/or modify
28  * it under the terms of the GNU General Public License as published by
29  * the Free Software Foundation; either version 3 of the License, or (at
30  * your option) any later version.
31  *
32  * This program is distributed in the hope that it will be useful, but
33  * WITHOUT ANY WARRANTY; without even the implied warranty of
34  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
35  * General Public License for more details.
36  *
37  * You should have received a copy of the GNU General Public License
38  * along with this program; if not, write to the Free Software
39  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
40  * 02110-1301, USA.
41  */
42 /** \file svdstep_base.h
43  \brief File for SVD decomposition
44 */
45 
46 #ifdef DOXYGEN
47 namespace o2scl_linalg {
48 #endif
49 
50  /** \brief Zero out small elements in \c f according to the
51  scales set in \c d
52 
53  The parameter \c N is the size of \c d.
54 
55  Used in \ref SV_decomp().
56  */
57  template<class vec_t, class vec2_t>
58  void chop_small_elements(size_t N, vec_t &d, vec2_t &f) {
59 
60  double d_i=O2SCL_IX(d,0);
61 
62  for (size_t i=0;i<N-1;i++) {
63 
64  double f_i=O2SCL_IX(f,i);
65  double d_ip1=O2SCL_IX(d,i+1);
66 
67  double dbl_eps=std::numeric_limits<double>::epsilon();
68 
69  if (fabs (f_i)<dbl_eps*(fabs(d_i)+fabs(d_ip1))) {
70  O2SCL_IX(f,i)=0.0;
71  }
72  d_i=d_ip1;
73  }
74 
75  return;
76  }
77 
78  /** \brief Desc
79 
80  The parameter \c n is the size of the vector \c d.
81 
82  Used in \ref qrstep() and \ref qrstep_sub().
83  */
84  template<class vec_t, class vec2_t>
85  double trailing_eigenvalue(size_t n, const vec_t &d, const vec2_t &f) {
86 
87  double da=O2SCL_IX(d,n-2);
88  double db=O2SCL_IX(d,n-1);
89  double fa=(n>2) ? O2SCL_IX(f,n-3) : 0.0;
90  double fb=O2SCL_IX(f,n-2);
91 
92  double ta=da*da+fa*fa;
93  double tb=db*db+fb*fb;
94  double tab=da*fb;
95 
96  double dt=(ta-tb)/2.0;
97 
98  double S=ta+tb;
99  double da2=da*da,db2=db*db;
100  double fa2=fa*fa,fb2=fb*fb;
101  double P=(da2*db2)+(fa2*db2)+(fa2*fb2);
102  double D=hypot(dt,tab);
103  double r1=S/2+D;
104 
105  double mu;
106  if (dt>=0) {
107  // [GSL] tb < ta, choose smaller root
108  mu=(r1>0) ? P / r1 : 0.0;
109  } else {
110  // [GSL] tb > ta, choose larger root
111  mu=r1;
112  }
113  return mu;
114  }
115 
116  /** \brief Desc
117 
118  Used in \ref svd2() and \ref svd2_sub().
119  */
120  void create_schur(double d0, double f0, double d1, double &c,
121  double &s) {
122 
123  double apq=2.0*d0*f0;
124 
125  if (d0 == 0 || f0 == 0) {
126  c=1.0;
127  s=0.0;
128  return;
129  }
130 
131  // [GSL] Check if we need to rescale to avoid underflow/overflow
132  if (fabs(d0) < GSL_SQRT_DBL_MIN || fabs(d0) > GSL_SQRT_DBL_MAX
133  || fabs(f0) < GSL_SQRT_DBL_MIN || fabs(f0) > GSL_SQRT_DBL_MAX
134  || fabs(d1) < GSL_SQRT_DBL_MIN || fabs(d1) > GSL_SQRT_DBL_MAX) {
135 
136  double scale;
137  int d0_exp,f0_exp;
138  frexp(d0,&d0_exp);
139  frexp(f0,&f0_exp);
140  // [GSL] Bring |d0*f0| into the range GSL_DBL_MIN to GSL_DBL_MAX
141  scale=ldexp(1.0,-(d0_exp+f0_exp)/4);
142  d0*=scale;
143  f0*=scale;
144  d1*=scale;
145  apq=2.0*d0*f0;
146  }
147 
148  if (apq != 0.0) {
149  double t;
150  double tau=(f0*f0+(d1+d0)*(d1-d0))/apq;
151 
152  if (tau >= 0.0) {
153  t=1.0/(tau+hypot(1.0,tau));
154  } else {
155  t=-1.0/(-tau+hypot(1.0,tau));
156  }
157 
158  c=1.0/hypot(1.0,t);
159  s=t*(c);
160  } else {
161  c=1.0;
162  s=0.0;
163  }
164  return;
165  }
166 
167  /** \brief 2-variable SVD
168 
169  The parameter \c M is the number of rows in \c U and \c N
170  is the number of rows in \c V. Both U and V assumed to have
171  two columns.
172 
173  Used in \ref qrstep().
174  */
175  template<class vec_t, class vec2_t, class mat_t, class mat2_t>
176  void svd2(size_t M, size_t N, vec_t &d, vec2_t &f, mat_t &U, mat2_t &V) {
177 
178  double c,s,a11,a12,a21,a22;
179 
180  double d0=O2SCL_IX(d,0);
181  double f0=O2SCL_IX(f,0);
182 
183  double d1=O2SCL_IX(d,1);
184 
185  if (d0 == 0.0) {
186 
187  // [GSL] Eliminate off-diagonal element in [0,f0;0,d1] to
188  // make [d,0;0,0]
189  o2scl_linalg::create_givens(f0,d1,c,s);
190 
191  // [GSL] compute B <= G^T B X, where X=[0,1;1,0]
192 
193  O2SCL_IX(d,0)=c*f0-s*d1;
194  O2SCL_IX(f,0)=s*f0+c*d1;
195  O2SCL_IX(d,1)=0.0;
196 
197  // [GSL] Compute U <= U G
198 
199  for (size_t i=0;i<M;i++) {
200 
201  double Uip=O2SCL_IX2(U,i,0);
202  double Uiq=O2SCL_IX2(U,i,1);
203  O2SCL_IX2(U,i,0)=c*Uip-s*Uiq;
204  O2SCL_IX2(U,i,1)=s*Uip+c*Uiq;
205  }
206 
207  // [GSL] Compute V <= V X
208 
209  double temp;
210  for(size_t ik=0;ik<N;ik++) {
211  temp=O2SCL_IX2(V,ik,0);
212  O2SCL_IX2(V,ik,0)=O2SCL_IX2(V,ik,1);
213  O2SCL_IX2(V,ik,1)=temp;
214  }
215 
216  return;
217 
218  } else if (d1 == 0.0) {
219 
220  // [GSL] Eliminate off-diagonal element in [d0,f0;0,0]
221 
222  o2scl_linalg::create_givens(d0,f0,c,s);
223 
224  // [GSL] compute B <= B G
225 
226  O2SCL_IX(d,0)=d0*c-f0*s;
227  O2SCL_IX(f,0)=0.0;
228 
229  // [GSL] Compute V <= V G
230 
231  for (size_t i=0;i<N;i++) {
232  double Vip=O2SCL_IX2(V,i,0);
233  double Viq=O2SCL_IX2(V,i,1);
234  O2SCL_IX2(V,i,0)=c*Vip-s*Viq;
235  O2SCL_IX2(V,i,1)=s*Vip+c*Viq;
236  }
237 
238  return;
239 
240  } else {
241 
242  // [GSL] Make columns orthogonal, A = [d0, f0; 0, d1] * G
243 
244  create_schur(d0,f0,d1,c,s);
245 
246  // [GSL] compute B <= B G
247 
248  a11=c*d0-s*f0;
249  a21=-s*d1;
250  a12=s*d0+c*f0;
251  a22=c*d1;
252 
253  // [GSL] Compute V <= V G
254 
255  for (size_t i=0;i<N;i++) {
256 
257  double Vip=O2SCL_IX2(V,i,0);
258  double Viq=O2SCL_IX2(V,i,1);
259  O2SCL_IX2(V,i,0)=c*Vip-s*Viq;
260  O2SCL_IX2(V,i,1)=s*Vip+c*Viq;
261  }
262 
263  // [GSL] Eliminate off-diagonal elements, bring column with largest
264  // norm to first column
265 
266  if (hypot(a11,a21) < hypot(a12,a22)) {
267 
268  double t1,t2;
269 
270  // [GSL] B <= B X
271 
272  t1=a11; a11=a12; a12=t1;
273  t2=a21; a21=a22; a22=t2;
274 
275  // [GSL] V <= V X
276 
277  double temp;
278  for(size_t ik=0;ik<N;ik++) {
279  temp=O2SCL_IX2(V,ik,0);
280  O2SCL_IX2(V,ik,0)=O2SCL_IX2(V,ik,1);
281  O2SCL_IX2(V,ik,1)=temp;
282  }
283  }
284 
285  o2scl_linalg::create_givens(a11,a21,c,s);
286 
287  // [GSL] compute B <= G^T B
288 
289  O2SCL_IX(d,0)=c*a11-s*a21;
290  O2SCL_IX(f,0)=c*a12-s*a22;
291  O2SCL_IX(d,1)=s*a12+c*a22;
292 
293  // [GSL] Compute U <= U G
294 
295  for (size_t i=0;i<M;i++) {
296  double Uip=O2SCL_IX2(U,i,0);
297  double Uiq=O2SCL_IX2(U,i,1);
298  O2SCL_IX2(U,i,0)=c*Uip-s*Uiq;
299  O2SCL_IX2(U,i,1)=s*Uip+c*Uiq;
300  }
301 
302  return;
303  }
304  }
305 
306  /** \brief Shifted 2-variable SVD
307 
308  The parameter \c M is the number of rows in \c U and \c N
309  is the number of rows in \c V. Both U and V assumed to have
310  two columns.
311 
312  Used in \ref qrstep_sub().
313  */
314  template<class vec_t, class vec2_t, class mat_t, class mat2_t>
315  void svd2_sub(size_t M, size_t N, vec_t &d, vec2_t &f, mat_t &U,
316  mat2_t &V, size_t a) {
317 
318  double c,s,a11,a12,a21,a22;
319 
320  double d0=O2SCL_IX(d,a);
321  double f0=O2SCL_IX(f,a);
322 
323  double d1=O2SCL_IX(d,a+1);
324 
325  if (d0 == 0.0) {
326 
327  // [GSL] Eliminate off-diagonal element in [0,f0;0,d1] to
328  // make [d,0;0,0]
329  o2scl_linalg::create_givens(f0,d1,c,s);
330 
331  // [GSL] compute B <= G^T B X, where X=[0,1;1,0]
332 
333  O2SCL_IX(d,a)=c*f0-s*d1;
334  O2SCL_IX(f,a)=s*f0+c*d1;
335  O2SCL_IX(d,a+1)=0.0;
336 
337  // [GSL] Compute U <= U G
338 
339  for (size_t i=0;i<M;i++) {
340 
341  double Uip=O2SCL_IX2(U,i,a);
342  double Uiq=O2SCL_IX2(U,i,a+1);
343  O2SCL_IX2(U,i,a)=c*Uip-s*Uiq;
344  O2SCL_IX2(U,i,a+1)=s*Uip+c*Uiq;
345  }
346 
347  // [GSL] Compute V <= V X
348 
349  double temp;
350  for(size_t ik=0;ik<N;ik++) {
351  temp=O2SCL_IX2(V,ik,a);
352  O2SCL_IX2(V,ik,a)=O2SCL_IX2(V,ik,a+1);
353  O2SCL_IX2(V,ik,a+1)=temp;
354  }
355 
356  return;
357 
358  } else if (d1 == 0.0) {
359 
360  // [GSL] Eliminate off-diagonal element in [d0,f0;0,0]
361 
362  o2scl_linalg::create_givens(d0,f0,c,s);
363 
364  // [GSL] compute B <= B G
365 
366  O2SCL_IX(d,a)=d0*c-f0*s;
367  O2SCL_IX(f,a)=0.0;
368 
369  // [GSL] Compute V <= V G
370 
371  for (size_t i=0;i<N;i++) {
372  double Vip=O2SCL_IX2(V,i,a);
373  double Viq=O2SCL_IX2(V,i,a+1);
374  O2SCL_IX2(V,i,a)=c*Vip-s*Viq;
375  O2SCL_IX2(V,i,a+1)=s*Vip+c*Viq;
376  }
377 
378  return;
379 
380  } else {
381 
382  // [GSL] Make columns orthogonal, A = [d0, f0; 0, d1] * G
383 
384  create_schur(d0,f0,d1,c,s);
385 
386  // [GSL] compute B <= B G
387 
388  a11=c*d0-s*f0;
389  a21=-s*d1;
390  a12=s*d0+c*f0;
391  a22=c*d1;
392 
393  // [GSL] Compute V <= V G
394 
395  for (size_t i=0;i<N;i++) {
396 
397  double Vip=O2SCL_IX2(V,i,a);
398  double Viq=O2SCL_IX2(V,i,a+1);
399  O2SCL_IX2(V,i,a)=c*Vip-s*Viq;
400  O2SCL_IX2(V,i,a+1)=s*Vip+c*Viq;
401  }
402 
403  // [GSL] Eliminate off-diagonal elements, bring column with largest
404  // norm to first column
405 
406  if (hypot(a11,a21)<hypot(a12,a22)) {
407 
408  double t1, t2;
409 
410  // [GSL] B <= B X
411 
412  t1=a11; a11=a12; a12=t1;
413  t2=a21; a21=a22; a22=t2;
414 
415  // [GSL] V <= V X
416 
417  double temp;
418  for(size_t ik=0;ik<N;ik++) {
419  temp=O2SCL_IX2(V,ik,a);
420  O2SCL_IX2(V,ik,a)=O2SCL_IX2(V,ik,a+1);
421  O2SCL_IX2(V,ik,a+1)=temp;
422  }
423  }
424 
425  o2scl_linalg::create_givens(a11,a21,c,s);
426 
427  // [GSL] compute B <= G^T B
428 
429  O2SCL_IX(d,a)=c*a11-s*a21;
430  O2SCL_IX(f,a)=c*a12-s*a22;
431  O2SCL_IX(d,a+1)=s*a12+c*a22;
432 
433  // [GSL] Compute U <= U G
434 
435  for (size_t i=0;i<M;i++) {
436  double Uip=O2SCL_IX2(U,i,a);
437  double Uiq=O2SCL_IX2(U,i,a+1);
438  O2SCL_IX2(U,i,a)=c*Uip-s*Uiq;
439  O2SCL_IX2(U,i,a+1)=s*Uip+c*Uiq;
440  }
441 
442  return;
443  }
444  }
445 
446  /** \brief Desc
447 
448  The vector \c d should be of size <tt>n</tt>, the vector \c f
449  should be of size <tt>n</tt>, and the matrix U should be of size
450  <tt>(M,n)</tt>
451 
452  Used in \ref qrstep() and \ref qrstep_sub().
453  */
454  template<class vec_t, class vec2_t, class mat_t>
455  void chase_out_intermediate_zero(size_t M, size_t n, vec_t &d,
456  vec2_t &f, mat_t &U, size_t k0) {
457 
458  double c, s;
459 
460  double x=O2SCL_IX(f,k0);
461  double y=O2SCL_IX(d,k0+1);
462 
463  for (size_t k=k0;k<n-1;k++) {
464 
465  o2scl_linalg::create_givens(y,-x,c,s);
466 
467  // [GSL] Compute U <= U G
468  for (size_t i=0; i < M; i++) {
469  double Uip=O2SCL_IX2(U,i,k0);
470  double Uiq=O2SCL_IX2(U,i,k+1);
471  //std::cout << "Uip,Uiq: " << Uip << " " << Uiq << std::endl;
472  O2SCL_IX2(U,i,k0)=c*Uip-s*Uiq;
473  O2SCL_IX2(U,i,k+1)=s*Uip+c*Uiq;
474  }
475 
476  // [GSL] compute B <= G^T B
477 
478  O2SCL_IX(d,k+1)=s*x+c*y;
479 
480  if (k == k0) {
481  O2SCL_IX(f,k)=c*x-s*y;
482  }
483 
484  if (k<n-2) {
485  double z=O2SCL_IX(f,k+1);
486  O2SCL_IX(f,k+1)=c*z;
487 
488  x=-s*z;
489  y=O2SCL_IX(d,k+2);
490  }
491  }
492 
493  return;
494  }
495 
496  /** \brief Desc
497 
498  The vector \c d should be of size <tt>n</tt>, the vector \c f
499  should be of size <tt>n</tt>, and the matrix V should be of size
500  <tt>(N,n)</tt>
501 
502  Used in \ref qrstep().
503  */
504  template<class vec_t, class vec2_t, class mat_t>
505  void chase_out_trailing_zero(size_t N, size_t n, vec_t &d,
506  vec2_t &f, mat_t &V) {
507 
508  double c, s;
509 
510  double x=O2SCL_IX(d,n-2);
511  double y=O2SCL_IX(f,n-2);
512 
513  for (size_t k=n-1;k-- > 0;) {
514 
516 
517  // [GSL] Compute V <= V G where G = [c, s ; -s, c]
518 
519  for (size_t i=0;i<N;i++) {
520  double Vip=O2SCL_IX2(V,i,k);
521  double Viq=O2SCL_IX2(V,i,n-1);
522  O2SCL_IX2(V,i,k)=c*Vip-s*Viq;
523  O2SCL_IX2(V,i,n-1)=s*Vip+c*Viq;
524  }
525 
526  // [GSL] compute B <= B G
527 
528  O2SCL_IX(d,k)=c*x-s*y;
529 
530  if (k==n-2) {
531  O2SCL_IX(f,k)=s*x+c*y;
532  }
533 
534  if (k>0) {
535  double z=O2SCL_IX(f,k-1);
536  O2SCL_IX(f,k-1)=c*z;
537  x=O2SCL_IX(d,k-1);
538  y=s*z;
539  }
540  }
541 
542  return;
543  }
544 
545  /** \brief Desc
546 
547  The vector \c d should be of size <tt>n</tt>, the vector \c f
548  should be of size <tt>n</tt>, and the matrix V should be of size
549  <tt>(N,n)</tt>
550 
551  Used in \ref qrstep_sub().
552  */
553  template<class vec_t, class vec2_t, class mat_t>
554  void chase_out_trailing_zero_sub(size_t N, size_t n, vec_t &d,
555  vec2_t &f, mat_t &V, size_t a) {
556 
557  double c, s;
558 
559  double x=O2SCL_IX(d,n-2);
560  double y=O2SCL_IX(f,n-2);
561 
562  for (int k=((int)n)-1;k>=((int)a);k--) {
563 
565 
566  // [GSL] Compute V <= V G where G = [c, s ; -s, c]
567 
568  for (size_t i=0;i<N;i++) {
569  double Vip=O2SCL_IX2(V,i,k);
570  double Viq=O2SCL_IX2(V,i,n-1);
571  O2SCL_IX2(V,i,k)=c*Vip-s*Viq;
572  O2SCL_IX2(V,i,n-1)=s*Vip+c*Viq;
573  }
574 
575  // [GSL] compute B <= B G
576 
577  O2SCL_IX(d,k)=c*x-s*y;
578 
579  if (k==((int)n)-2) {
580  O2SCL_IX(f,k)=s*x+c*y;
581  }
582 
583  if (k>((int)a)) {
584  double z=O2SCL_IX(f,k-1);
585  O2SCL_IX(f,k-1)=c*z;
586  x=O2SCL_IX(d,k-1);
587  y=s*z;
588  }
589  }
590 
591  return;
592  }
593 
594  /** \brief Desc
595 
596  The vector \c d should be of size \c n, the vector \c f should
597  be of size \c n, the matrix U should be of size <tt>(M,N)</tt>,
598  and the matrix \c V should be of size <tt>(N,N)</tt>.
599  */
600  template<class vec_t, class vec2_t, class mat_t, class mat2_t>
601  void qrstep(size_t M, size_t N, size_t n,
602  vec_t &d, vec2_t &f, mat_t &U, mat2_t &V) {
603 
604  double y, z;
605  double ak, bk, zk, ap, bp, aq, bq;
606  size_t i, k;
607 
608  if (n==1) {
609  // [GSL] shouldn't happen
610  return;
611  }
612 
613  // [GSL] Compute 2x2 svd directly
614 
615  if (n==2) {
616  svd2(M,N,d,f,U,V);
617  return;
618  }
619 
620  // [GSL] Chase out any zeroes on the diagonal
621 
622  for (i=0;i<n-1;i++) {
623  double d_i=O2SCL_IX(d,i);
624  if (d_i==0.0) {
625  chase_out_intermediate_zero(M,n,d,f,U,i);
626  return;
627  }
628  }
629 
630  // [GSL] Chase out any zero at the end of the diagonal
631  double d_nm1=O2SCL_IX(d,n-1);
632 
633  if (d_nm1==0.0) {
634  chase_out_trailing_zero(N,n,d,f,V);
635  return;
636  }
637 
638  // [GSL] Apply QR reduction steps to the diagonal and offdiagonal
639 
640  double d0=O2SCL_IX(d,0);
641  double f0=O2SCL_IX(f,0);
642 
643  double d1=O2SCL_IX(d,1);
644  double f1=O2SCL_IX(f,1);
645 
646  double mu=trailing_eigenvalue(n,d,f);
647  y=d0*d0-mu;
648  z=d0*f0;
649 
650  // [GSL] Set up the recurrence for Givens rotations on a bidiagonal
651  // matrix
652 
653  ak=0;
654  bk=0;
655 
656  ap=d0;
657  bp=f0;
658 
659  aq=d1;
660  bq=f1;
661 
662  for (k=0; k < n-1; k++) {
663 
664  double c, s;
666 
667  // [GSL] Compute V <= V G
668 
669  for (i=0; i < N; i++) {
670  double Vip=O2SCL_IX2(V,i,k);
671  double Viq=O2SCL_IX2(V,i,k+1);
672  O2SCL_IX2(V,i,k)=c*Vip-s*Viq;
673  O2SCL_IX2(V,i,k+1)=s*Vip+c*Viq;
674  }
675 
676  // [GSL] compute B <= B G
677 
678  {
679  double bk1=c*bk-s*z;
680 
681  double ap1=c*ap-s*bp;
682  double bp1=s*ap+c*bp;
683  double zp1=-s*aq;
684 
685  double aq1=c*aq;
686 
687  if (k > 0) O2SCL_IX(f,k-1)=bk1;
688 
689  ak=ap1;
690  bk=bp1;
691  zk=zp1;
692 
693  ap=aq1;
694 
695  if (k<n-2) bp=O2SCL_IX(f,k+1);
696  else bp=0.0;
697 
698  y=ak;
699  z=zk;
700  }
701 
703 
704  // [GSL] Compute U <= U G
705 
706  for (i=0;i<M;i++) {
707  double Uip=O2SCL_IX2(U,i,k);
708  double Uiq=O2SCL_IX2(U,i,k+1);
709  O2SCL_IX2(U,i,k)=c*Uip-s*Uiq;
710  O2SCL_IX2(U,i,k+1)=s*Uip+c*Uiq;
711  }
712 
713  // [GSL] compute B <= G^T B
714 
715  double ak1=c*ak-s*zk;
716  double bk1=c*bk-s*ap;
717  double zk1=-s*bp;
718 
719  double ap1=s*bk+c*ap;
720  double bp1=c*bp;
721 
722  O2SCL_IX(d,k)=ak1;
723 
724  ak=ak1;
725  bk=bk1;
726  zk=zk1;
727 
728  ap=ap1;
729  bp=bp1;
730 
731  if (k < n-2) {
732  aq=O2SCL_IX(d,k+2);
733  } else {
734  aq=0.0;
735  }
736 
737  y=bk;
738  z=zk;
739  }
740 
741  O2SCL_IX(f,n-2)=bk;
742  O2SCL_IX(d,n-1)=ap;
743 
744  return;
745  }
746 
747  /** \brief A special form of qrstep() for SV_decomp()
748 
749  The vector \c d should be of size <tt>n</tt>, the vector \c f
750  should be of size <tt>n</tt>, the matrix U should be of size
751  <tt>(M,n)</tt>, and the matrix \c V should be of size
752  <tt>(N,n)</tt>.
753 
754  A version of qrstep(), but starting at the a'th column of U, the
755  a'th column of V, and the a'th entries of \c d and \c f.
756 
757  This function is used by \ref SV_decomp().
758  */
759  template<class vec_t, class vec2_t, class mat_t, class mat2_t>
760  void qrstep_sub(size_t M, size_t N, size_t n,
761  vec_t &d, vec2_t &f, mat_t &U, mat2_t &V, size_t a) {
762 
763  double y, z;
764  double ak, bk, zk, ap, bp, aq, bq;
765  size_t i, k;
766 
767  //std::cout << "M,N,n: " << M << " " << N << " " << n << std::endl;
768 
769  if (n-a==1) {
770  // [GSL] shouldn't happen
771  return;
772  }
773 
774  // [GSL] Compute 2x2 svd directly
775 
776  if (n-a==2) {
777  svd2_sub(M,N,d,f,U,V,a);
778  return;
779  }
780 
781  // [GSL] Chase out any zeroes on the diagonal
782 
783  for (i=a;i<n-1;i++) {
784  double d_i=O2SCL_IX(d,i);
785  //std::cout << "d_i: " << i << " " << n << " "
786  //<< d_i << std::endl;
787  if (d_i==0.0) {
788  chase_out_intermediate_zero(M,n,d,f,U,i);
789  return;
790  }
791  }
792 
793  // [GSL] Chase out any zero at the end of the diagonal
794  double d_nm1=O2SCL_IX(d,n-1);
795  //std::cout << "d_nm1: " << d_nm1 << std::endl;
796  if (d_nm1==0.0) {
797  chase_out_trailing_zero_sub(N,n,d,f,V,a);
798  return;
799  }
800 
801  // [GSL] Apply QR reduction steps to the diagonal and offdiagonal
802 
803  double d0=O2SCL_IX(d,a);
804  double f0=O2SCL_IX(f,a);
805 
806  double d1=O2SCL_IX(d,a+1);
807  double f1=O2SCL_IX(f,a+1);
808 
809  //std::cout << "d0,f0,d1,f1: " << d0 << " " << f0 << " " << d1 << " "
810  //<< f1 << std::endl;
811 
812  double mu=trailing_eigenvalue(n,d,f);
813  y=d0*d0-mu;
814  z=d0*f0;
815 
816  // [GSL] Set up the recurrence for Givens rotations on a bidiagonal
817  // matrix
818 
819  ak=0;
820  bk=0;
821 
822  ap=d0;
823  bp=f0;
824 
825  aq=d1;
826  bq=f1;
827 
828  for (k=a;k<n-1;k++) {
829 
830  double c, s;
832 
833  // [GSL] Compute V <= V G
834 
835  for (i=0;i<N;i++) {
836  double Vip=O2SCL_IX2(V,i,k);
837  double Viq=O2SCL_IX2(V,i,k+1);
838  //std::cout << "Vip,Viq: " << Vip << " " << Viq << std::endl;
839  O2SCL_IX2(V,i,k)=c*Vip-s*Viq;
840  O2SCL_IX2(V,i,k+1)=s*Vip+c*Viq;
841  }
842 
843  // [GSL] compute B <= B G
844 
845  {
846  double bk1=c*bk-s*z;
847 
848  double ap1=c*ap-s*bp;
849  double bp1=s*ap+c*bp;
850  double zp1=-s*aq;
851 
852  double aq1=c*aq;
853 
854  if (k>a) O2SCL_IX(f,k-1)=bk1;
855 
856  ak=ap1;
857  bk=bp1;
858  zk=zp1;
859 
860  ap=aq1;
861 
862  if (k<n-2) bp=O2SCL_IX(f,k+1);
863  else bp=0.0;
864 
865  y=ak;
866  z=zk;
867  }
868 
870 
871  // [GSL] Compute U <= U G
872 
873  for (i=0;i<M;i++) {
874  double Uip=O2SCL_IX2(U,i,k);
875  double Uiq=O2SCL_IX2(U,i,k+1);
876  //std::cout << "Uip2,Uiq2: " << Uip << " " << Uiq << std::endl;
877  O2SCL_IX2(U,i,k)=c*Uip-s*Uiq;
878  O2SCL_IX2(U,i,k+1)=s*Uip+c*Uiq;
879  }
880 
881  // [GSL] compute B <= G^T B
882 
883  //std::cout << "k,bk,ap2: " << k << " " << bk << " " << ap << std::endl;
884  //std::cout << "ak,zk,bp: " << ak << " " << zk << " "
885  //<< bp << std::endl;
886 
887  {
888  //std::cout << "prod1: " << c*ak << " " << s*zk << std::endl;
889  //std::cout << "prod2: " << c*bk << " " << s*ap << std::endl;
890  //std::cout << "prod3: " << s*bk << " " << c*ap << std::endl;
891  double ak1=c*ak-s*zk;
892  double bk1=c*bk-s*ap;
893  double zk1=-s*bp;
894 
895  double ap1=s*bk+c*ap;
896  double bp1=c*bp;
897 
898  O2SCL_IX(d,k)=ak1;
899 
900  ak=ak1;
901  bk=bk1;
902  zk=zk1;
903 
904  ap=ap1;
905  bp=bp1;
906  //std::cout << "c,s: " << c << " " << s << std::endl;
907  //std::cout << "k,bk,ap: " << k << " " << bk << " " << ap << std::endl;
908 
909  if (k < n-2) {
910  aq=O2SCL_IX(d,k+2);
911  } else {
912  aq=0.0;
913  }
914 
915  y=bk;
916  z=zk;
917  }
918  }
919 
920  O2SCL_IX(f,n-2)=bk;
921  O2SCL_IX(d,n-1)=ap;
922  //std::cout << "bk,ap: " << bk << " " << ap << std::endl;
923 
924  return;
925  }
926 
927 #ifdef DOXYGEN
928 }
929 #endif
o2scl_linalg
The namespace for linear algebra classes and functions.
Definition: bidiag.h:36
o2scl_linalg::svd2
void svd2(size_t M, size_t N, vec_t &d, vec2_t &f, mat_t &U, mat2_t &V)
2-variable SVD
Definition: svdstep_base.h:176
o2scl_linalg::create_schur
void create_schur(double d0, double f0, double d1, double &c, double &s)
Desc.
Definition: svdstep_base.h:120
o2scl_linalg::svd2_sub
void svd2_sub(size_t M, size_t N, vec_t &d, vec2_t &f, mat_t &U, mat2_t &V, size_t a)
Shifted 2-variable SVD.
Definition: svdstep_base.h:315
o2scl_linalg::chase_out_intermediate_zero
void chase_out_intermediate_zero(size_t M, size_t n, vec_t &d, vec2_t &f, mat_t &U, size_t k0)
Desc.
Definition: svdstep_base.h:455
o2scl_linalg::trailing_eigenvalue
double trailing_eigenvalue(size_t n, const vec_t &d, const vec2_t &f)
Desc.
Definition: svdstep_base.h:85
o2scl_linalg::create_givens
void create_givens(const double a, const double b, double &c, double &s)
Create a Givens rotation matrix.
o2scl_linalg::chop_small_elements
void chop_small_elements(size_t N, vec_t &d, vec2_t &f)
Zero out small elements in f according to the scales set in d.
Definition: svdstep_base.h:58
o2scl_linalg::qrstep
void qrstep(size_t M, size_t N, size_t n, vec_t &d, vec2_t &f, mat_t &U, mat2_t &V)
Desc.
Definition: svdstep_base.h:601
o2scl_linalg::qrstep_sub
void qrstep_sub(size_t M, size_t N, size_t n, vec_t &d, vec2_t &f, mat_t &U, mat2_t &V, size_t a)
A special form of qrstep() for SV_decomp()
Definition: svdstep_base.h:760
o2scl_linalg::chase_out_trailing_zero
void chase_out_trailing_zero(size_t N, size_t n, vec_t &d, vec2_t &f, mat_t &V)
Desc.
Definition: svdstep_base.h:505
o2scl_linalg::chase_out_trailing_zero_sub
void chase_out_trailing_zero_sub(size_t N, size_t n, vec_t &d, vec2_t &f, mat_t &V, size_t a)
Desc.
Definition: svdstep_base.h:554

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).