svdstep_base.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2008-2017, 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  size_t i;
179  double c,s,a11,a12,a21,a22;
180 
181  double d0=O2SCL_IX(d,0);
182  double f0=O2SCL_IX(f,0);
183 
184  double d1=O2SCL_IX(d,1);
185 
186  if (d0 == 0.0) {
187 
188  // [GSL] Eliminate off-diagonal element in [0,f0;0,d1] to
189  // make [d,0;0,0]
190  o2scl_linalg::create_givens(f0,d1,c,s);
191 
192  // [GSL] compute B <= G^T B X, where X=[0,1;1,0]
193 
194  O2SCL_IX(d,0)=c*f0-s*d1;
195  O2SCL_IX(f,0)=s*f0+c*d1;
196  O2SCL_IX(d,1)=0.0;
197 
198  // [GSL] Compute U <= U G
199 
200  for (size_t i=0;i<M;i++) {
201 
202  double Uip=O2SCL_IX2(U,i,0);
203  double Uiq=O2SCL_IX2(U,i,1);
204  O2SCL_IX2(U,i,0)=c*Uip-s*Uiq;
205  O2SCL_IX2(U,i,1)=s*Uip+c*Uiq;
206  }
207 
208  // [GSL] Compute V <= V X
209 
210  double temp;
211  for(size_t ik=0;ik<N;ik++) {
212  temp=O2SCL_IX2(V,ik,0);
213  O2SCL_IX2(V,ik,0)=O2SCL_IX2(V,ik,1);
214  O2SCL_IX2(V,ik,1)=temp;
215  }
216 
217  return;
218 
219  } else if (d1 == 0.0) {
220 
221  // [GSL] Eliminate off-diagonal element in [d0,f0;0,0]
222 
223  o2scl_linalg::create_givens(d0,f0,c,s);
224 
225  // [GSL] compute B <= B G
226 
227  O2SCL_IX(d,0)=d0*c-f0*s;
228  O2SCL_IX(f,0)=0.0;
229 
230  // [GSL] Compute V <= V G
231 
232  for (size_t i=0;i<N;i++) {
233  double Vip=O2SCL_IX2(V,i,0);
234  double Viq=O2SCL_IX2(V,i,1);
235  O2SCL_IX2(V,i,0)=c*Vip-s*Viq;
236  O2SCL_IX2(V,i,1)=s*Vip+c*Viq;
237  }
238 
239  return;
240 
241  } else {
242 
243  // [GSL] Make columns orthogonal, A = [d0, f0; 0, d1] * G
244 
245  create_schur(d0,f0,d1,c,s);
246 
247  // [GSL] compute B <= B G
248 
249  a11=c*d0-s*f0;
250  a21=-s*d1;
251  a12=s*d0+c*f0;
252  a22=c*d1;
253 
254  // [GSL] Compute V <= V G
255 
256  for (size_t i=0;i<N;i++) {
257 
258  double Vip=O2SCL_IX2(V,i,0);
259  double Viq=O2SCL_IX2(V,i,1);
260  O2SCL_IX2(V,i,0)=c*Vip-s*Viq;
261  O2SCL_IX2(V,i,1)=s*Vip+c*Viq;
262  }
263 
264  // [GSL] Eliminate off-diagonal elements, bring column with largest
265  // norm to first column
266 
267  if (hypot(a11,a21) < hypot(a12,a22)) {
268 
269  double t1,t2;
270 
271  // [GSL] B <= B X
272 
273  t1=a11; a11=a12; a12=t1;
274  t2=a21; a21=a22; a22=t2;
275 
276  // [GSL] V <= V X
277 
278  double temp;
279  for(size_t ik=0;ik<N;ik++) {
280  temp=O2SCL_IX2(V,ik,0);
281  O2SCL_IX2(V,ik,0)=O2SCL_IX2(V,ik,1);
282  O2SCL_IX2(V,ik,1)=temp;
283  }
284  }
285 
286  o2scl_linalg::create_givens(a11,a21,c,s);
287 
288  // [GSL] compute B <= G^T B
289 
290  O2SCL_IX(d,0)=c*a11-s*a21;
291  O2SCL_IX(f,0)=c*a12-s*a22;
292  O2SCL_IX(d,1)=s*a12+c*a22;
293 
294  // [GSL] Compute U <= U G
295 
296  for (size_t i=0;i<M;i++) {
297  double Uip=O2SCL_IX2(U,i,0);
298  double Uiq=O2SCL_IX2(U,i,1);
299  O2SCL_IX2(U,i,0)=c*Uip-s*Uiq;
300  O2SCL_IX2(U,i,1)=s*Uip+c*Uiq;
301  }
302 
303  return;
304  }
305  }
306 
307  /** \brief Shifted 2-variable SVD
308 
309  The parameter \c M is the number of rows in \c U and \c N
310  is the number of rows in \c V. Both U and V assumed to have
311  two columns.
312 
313  Used in \ref qrstep_sub().
314  */
315  template<class vec_t, class vec2_t, class mat_t, class mat2_t>
316  void svd2_sub(size_t M, size_t N, vec_t &d, vec2_t &f, mat_t &U,
317  mat2_t &V, size_t a) {
318 
319  size_t i;
320  double c,s,a11,a12,a21,a22;
321 
322  double d0=O2SCL_IX(d,a);
323  double f0=O2SCL_IX(f,a);
324 
325  double d1=O2SCL_IX(d,a+1);
326 
327  if (d0 == 0.0) {
328 
329  // [GSL] Eliminate off-diagonal element in [0,f0;0,d1] to
330  // make [d,0;0,0]
331  o2scl_linalg::create_givens(f0,d1,c,s);
332 
333  // [GSL] compute B <= G^T B X, where X=[0,1;1,0]
334 
335  O2SCL_IX(d,a)=c*f0-s*d1;
336  O2SCL_IX(f,a)=s*f0+c*d1;
337  O2SCL_IX(d,a+1)=0.0;
338 
339  // [GSL] Compute U <= U G
340 
341  for (size_t i=0;i<M;i++) {
342 
343  double Uip=O2SCL_IX2(U,i,a);
344  double Uiq=O2SCL_IX2(U,i,a+1);
345  O2SCL_IX2(U,i,a)=c*Uip-s*Uiq;
346  O2SCL_IX2(U,i,a+1)=s*Uip+c*Uiq;
347  }
348 
349  // [GSL] Compute V <= V X
350 
351  double temp;
352  for(size_t ik=0;ik<N;ik++) {
353  temp=O2SCL_IX2(V,ik,a);
354  O2SCL_IX2(V,ik,a)=O2SCL_IX2(V,ik,a+1);
355  O2SCL_IX2(V,ik,a+1)=temp;
356  }
357 
358  return;
359 
360  } else if (d1 == 0.0) {
361 
362  // [GSL] Eliminate off-diagonal element in [d0,f0;0,0]
363 
364  o2scl_linalg::create_givens(d0,f0,c,s);
365 
366  // [GSL] compute B <= B G
367 
368  O2SCL_IX(d,a)=d0*c-f0*s;
369  O2SCL_IX(f,a)=0.0;
370 
371  // [GSL] Compute V <= V G
372 
373  for (size_t i=0;i<N;i++) {
374  double Vip=O2SCL_IX2(V,i,a);
375  double Viq=O2SCL_IX2(V,i,a+1);
376  O2SCL_IX2(V,i,a)=c*Vip-s*Viq;
377  O2SCL_IX2(V,i,a+1)=s*Vip+c*Viq;
378  }
379 
380  return;
381 
382  } else {
383 
384  // [GSL] Make columns orthogonal, A = [d0, f0; 0, d1] * G
385 
386  create_schur(d0,f0,d1,c,s);
387 
388  // [GSL] compute B <= B G
389 
390  a11=c*d0-s*f0;
391  a21=-s*d1;
392  a12=s*d0+c*f0;
393  a22=c*d1;
394 
395  // [GSL] Compute V <= V G
396 
397  for (size_t i=0;i<N;i++) {
398 
399  double Vip=O2SCL_IX2(V,i,a);
400  double Viq=O2SCL_IX2(V,i,a+1);
401  O2SCL_IX2(V,i,a)=c*Vip-s*Viq;
402  O2SCL_IX2(V,i,a+1)=s*Vip+c*Viq;
403  }
404 
405  // [GSL] Eliminate off-diagonal elements, bring column with largest
406  // norm to first column
407 
408  if (hypot(a11,a21)<hypot(a12,a22)) {
409 
410  double t1, t2;
411 
412  // [GSL] B <= B X
413 
414  t1=a11; a11=a12; a12=t1;
415  t2=a21; a21=a22; a22=t2;
416 
417  // [GSL] V <= V X
418 
419  double temp;
420  for(size_t ik=0;ik<N;ik++) {
421  temp=O2SCL_IX2(V,ik,a);
422  O2SCL_IX2(V,ik,a)=O2SCL_IX2(V,ik,a+1);
423  O2SCL_IX2(V,ik,a+1)=temp;
424  }
425  }
426 
427  o2scl_linalg::create_givens(a11,a21,c,s);
428 
429  // [GSL] compute B <= G^T B
430 
431  O2SCL_IX(d,a)=c*a11-s*a21;
432  O2SCL_IX(f,a)=c*a12-s*a22;
433  O2SCL_IX(d,a+1)=s*a12+c*a22;
434 
435  // [GSL] Compute U <= U G
436 
437  for (size_t i=0;i<M;i++) {
438  double Uip=O2SCL_IX2(U,i,a);
439  double Uiq=O2SCL_IX2(U,i,a+1);
440  O2SCL_IX2(U,i,a)=c*Uip-s*Uiq;
441  O2SCL_IX2(U,i,a+1)=s*Uip+c*Uiq;
442  }
443 
444  return;
445  }
446  }
447 
448  /** \brief Desc
449 
450  The vector \c d should be of size <tt>n</tt>, the vector \c f
451  should be of size <tt>n</tt>, and the matrix U should be of size
452  <tt>(M,n)</tt>
453 
454  Used in \ref qrstep() and \ref qrstep_sub().
455  */
456  template<class vec_t, class vec2_t, class mat_t>
457  void chase_out_intermediate_zero(size_t M, size_t n, vec_t &d,
458  vec2_t &f, mat_t &U, size_t k0) {
459 
460  double c, s;
461 
462  double x=O2SCL_IX(f,k0);
463  double y=O2SCL_IX(d,k0+1);
464 
465  for (size_t k=k0;k<n-1;k++) {
466 
467  o2scl_linalg::create_givens(y,-x,c,s);
468 
469  // [GSL] Compute U <= U G
470  for (size_t i=0; i < M; i++) {
471  double Uip=O2SCL_IX2(U,i,k0);
472  double Uiq=O2SCL_IX2(U,i,k+1);
473  //std::cout << "Uip,Uiq: " << Uip << " " << Uiq << std::endl;
474  O2SCL_IX2(U,i,k0)=c*Uip-s*Uiq;
475  O2SCL_IX2(U,i,k+1)=s*Uip+c*Uiq;
476  }
477 
478  // [GSL] compute B <= G^T B
479 
480  O2SCL_IX(d,k+1)=s*x+c*y;
481 
482  if (k == k0) {
483  O2SCL_IX(f,k)=c*x-s*y;
484  }
485 
486  if (k<n-2) {
487  double z=O2SCL_IX(f,k+1);
488  O2SCL_IX(f,k+1)=c*z;
489 
490  x=-s*z;
491  y=O2SCL_IX(d,k+2);
492  }
493  }
494 
495  return;
496  }
497 
498  /** \brief Desc
499 
500  The vector \c d should be of size <tt>n</tt>, the vector \c f
501  should be of size <tt>n</tt>, and the matrix V should be of size
502  <tt>(N,n)</tt>
503 
504  Used in \ref qrstep().
505  */
506  template<class vec_t, class vec2_t, class mat_t>
507  void chase_out_trailing_zero(size_t N, size_t n, vec_t &d,
508  vec2_t &f, mat_t &V) {
509 
510  double c, s;
511 
512  double x=O2SCL_IX(d,n-2);
513  double y=O2SCL_IX(f,n-2);
514 
515  for (size_t k=n-1;k-- > 0;) {
516 
518 
519  // [GSL] Compute V <= V G where G = [c, s ; -s, c]
520 
521  for (size_t i=0;i<N;i++) {
522  double Vip=O2SCL_IX2(V,i,k);
523  double Viq=O2SCL_IX2(V,i,n-1);
524  O2SCL_IX2(V,i,k)=c*Vip-s*Viq;
525  O2SCL_IX2(V,i,n-1)=s*Vip+c*Viq;
526  }
527 
528  // [GSL] compute B <= B G
529 
530  O2SCL_IX(d,k)=c*x-s*y;
531 
532  if (k==n-2) {
533  O2SCL_IX(f,k)=s*x+c*y;
534  }
535 
536  if (k>0) {
537  double z=O2SCL_IX(f,k-1);
538  O2SCL_IX(f,k-1)=c*z;
539  x=O2SCL_IX(d,k-1);
540  y=s*z;
541  }
542  }
543 
544  return;
545  }
546 
547  /** \brief Desc
548 
549  The vector \c d should be of size <tt>n</tt>, the vector \c f
550  should be of size <tt>n</tt>, and the matrix V should be of size
551  <tt>(N,n)</tt>
552 
553  Used in \ref qrstep_sub().
554  */
555  template<class vec_t, class vec2_t, class mat_t>
556  void chase_out_trailing_zero_sub(size_t N, size_t n, vec_t &d,
557  vec2_t &f, mat_t &V, size_t a) {
558 
559  double c, s;
560 
561  double x=O2SCL_IX(d,n-2);
562  double y=O2SCL_IX(f,n-2);
563 
564  for (int k=((int)n)-1;k>=((int)a);k--) {
565 
567 
568  // [GSL] Compute V <= V G where G = [c, s ; -s, c]
569 
570  for (size_t i=0;i<N;i++) {
571  double Vip=O2SCL_IX2(V,i,k);
572  double Viq=O2SCL_IX2(V,i,n-1);
573  O2SCL_IX2(V,i,k)=c*Vip-s*Viq;
574  O2SCL_IX2(V,i,n-1)=s*Vip+c*Viq;
575  }
576 
577  // [GSL] compute B <= B G
578 
579  O2SCL_IX(d,k)=c*x-s*y;
580 
581  if (k==((int)n)-2) {
582  O2SCL_IX(f,k)=s*x+c*y;
583  }
584 
585  if (k>((int)a)) {
586  double z=O2SCL_IX(f,k-1);
587  O2SCL_IX(f,k-1)=c*z;
588  x=O2SCL_IX(d,k-1);
589  y=s*z;
590  }
591  }
592 
593  return;
594  }
595 
596  /** \brief Desc
597 
598  The vector \c d should be of size \c n, the vector \c f should
599  be of size \c n, the matrix U should be of size <tt>(M,N)</tt>,
600  and the matrix \c V should be of size <tt>(N,N)</tt>.
601  */
602  template<class vec_t, class vec2_t, class mat_t, class mat2_t>
603  void qrstep(size_t M, size_t N, size_t n,
604  vec_t &d, vec2_t &f, mat_t &U, mat2_t &V) {
605 
606  double y, z;
607  double ak, bk, zk, ap, bp, aq, bq;
608  size_t i, k;
609 
610  if (n==1) {
611  // [GSL] shouldn't happen
612  return;
613  }
614 
615  // [GSL] Compute 2x2 svd directly
616 
617  if (n==2) {
618  svd2(M,N,d,f,U,V);
619  return;
620  }
621 
622  // [GSL] Chase out any zeroes on the diagonal
623 
624  for (i=0;i<n-1;i++) {
625  double d_i=O2SCL_IX(d,i);
626  if (d_i==0.0) {
627  chase_out_intermediate_zero(M,n,d,f,U,i);
628  return;
629  }
630  }
631 
632  // [GSL] Chase out any zero at the end of the diagonal
633  double d_nm1=O2SCL_IX(d,n-1);
634 
635  if (d_nm1==0.0) {
636  chase_out_trailing_zero(N,n,d,f,V);
637  return;
638  }
639 
640  // [GSL] Apply QR reduction steps to the diagonal and offdiagonal
641 
642  double d0=O2SCL_IX(d,0);
643  double f0=O2SCL_IX(f,0);
644 
645  double d1=O2SCL_IX(d,1);
646  double f1=O2SCL_IX(f,1);
647 
648  double mu=trailing_eigenvalue(n,d,f);
649  y=d0*d0-mu;
650  z=d0*f0;
651 
652  // [GSL] Set up the recurrence for Givens rotations on a bidiagonal
653  // matrix
654 
655  ak=0;
656  bk=0;
657 
658  ap=d0;
659  bp=f0;
660 
661  aq=d1;
662  bq=f1;
663 
664  for (k=0; k < n-1; k++) {
665 
666  double c, s;
668 
669  // [GSL] Compute V <= V G
670 
671  for (i=0; i < N; i++) {
672  double Vip=O2SCL_IX2(V,i,k);
673  double Viq=O2SCL_IX2(V,i,k+1);
674  O2SCL_IX2(V,i,k)=c*Vip-s*Viq;
675  O2SCL_IX2(V,i,k+1)=s*Vip+c*Viq;
676  }
677 
678  // [GSL] compute B <= B G
679 
680  {
681  double bk1=c*bk-s*z;
682 
683  double ap1=c*ap-s*bp;
684  double bp1=s*ap+c*bp;
685  double zp1=-s*aq;
686 
687  double aq1=c*aq;
688 
689  if (k > 0) O2SCL_IX(f,k-1)=bk1;
690 
691  ak=ap1;
692  bk=bp1;
693  zk=zp1;
694 
695  ap=aq1;
696 
697  if (k<n-2) bp=O2SCL_IX(f,k+1);
698  else bp=0.0;
699 
700  y=ak;
701  z=zk;
702  }
703 
705 
706  // [GSL] Compute U <= U G
707 
708  for (i=0;i<M;i++) {
709  double Uip=O2SCL_IX2(U,i,k);
710  double Uiq=O2SCL_IX2(U,i,k+1);
711  O2SCL_IX2(U,i,k)=c*Uip-s*Uiq;
712  O2SCL_IX2(U,i,k+1)=s*Uip+c*Uiq;
713  }
714 
715  // [GSL] compute B <= G^T B
716 
717  double ak1=c*ak-s*zk;
718  double bk1=c*bk-s*ap;
719  double zk1=-s*bp;
720 
721  double ap1=s*bk+c*ap;
722  double bp1=c*bp;
723 
724  O2SCL_IX(d,k)=ak1;
725 
726  ak=ak1;
727  bk=bk1;
728  zk=zk1;
729 
730  ap=ap1;
731  bp=bp1;
732 
733  if (k < n-2) {
734  aq=O2SCL_IX(d,k+2);
735  } else {
736  aq=0.0;
737  }
738 
739  y=bk;
740  z=zk;
741  }
742 
743  O2SCL_IX(f,n-2)=bk;
744  O2SCL_IX(d,n-1)=ap;
745 
746  return;
747  }
748 
749  /** \brief A special form of qrstep() for SV_decomp()
750 
751  The vector \c d should be of size <tt>n</tt>, the vector \c f
752  should be of size <tt>n</tt>, the matrix U should be of size
753  <tt>(M,n)</tt>, and the matrix \c V should be of size
754  <tt>(N,n)</tt>.
755 
756  A version of qrstep(), but starting at the a'th column of U, the
757  a'th column of V, and the a'th entries of \c d and \c f.
758 
759  This function is used by \ref SV_decomp().
760  */
761  template<class vec_t, class vec2_t, class mat_t, class mat2_t>
762  void qrstep_sub(size_t M, size_t N, size_t n,
763  vec_t &d, vec2_t &f, mat_t &U, mat2_t &V, size_t a) {
764 
765  double y, z;
766  double ak, bk, zk, ap, bp, aq, bq;
767  size_t i, k;
768 
769  //std::cout << "M,N,n: " << M << " " << N << " " << n << std::endl;
770 
771  if (n-a==1) {
772  // [GSL] shouldn't happen
773  return;
774  }
775 
776  // [GSL] Compute 2x2 svd directly
777 
778  if (n-a==2) {
779  svd2_sub(M,N,d,f,U,V,a);
780  return;
781  }
782 
783  // [GSL] Chase out any zeroes on the diagonal
784 
785  for (i=a;i<n-1;i++) {
786  double d_i=O2SCL_IX(d,i);
787  //std::cout << "d_i: " << i << " " << n << " "
788  //<< d_i << std::endl;
789  if (d_i==0.0) {
790  chase_out_intermediate_zero(M,n,d,f,U,i);
791  return;
792  }
793  }
794 
795  // [GSL] Chase out any zero at the end of the diagonal
796  double d_nm1=O2SCL_IX(d,n-1);
797  //std::cout << "d_nm1: " << d_nm1 << std::endl;
798  if (d_nm1==0.0) {
799  chase_out_trailing_zero_sub(N,n,d,f,V,a);
800  return;
801  }
802 
803  // [GSL] Apply QR reduction steps to the diagonal and offdiagonal
804 
805  double d0=O2SCL_IX(d,a);
806  double f0=O2SCL_IX(f,a);
807 
808  double d1=O2SCL_IX(d,a+1);
809  double f1=O2SCL_IX(f,a+1);
810 
811  //std::cout << "d0,f0,d1,f1: " << d0 << " " << f0 << " " << d1 << " "
812  //<< f1 << std::endl;
813 
814  double mu=trailing_eigenvalue(n,d,f);
815  y=d0*d0-mu;
816  z=d0*f0;
817 
818  // [GSL] Set up the recurrence for Givens rotations on a bidiagonal
819  // matrix
820 
821  ak=0;
822  bk=0;
823 
824  ap=d0;
825  bp=f0;
826 
827  aq=d1;
828  bq=f1;
829 
830  for (k=a;k<n-1;k++) {
831 
832  double c, s;
834 
835  // [GSL] Compute V <= V G
836 
837  for (i=0;i<N;i++) {
838  double Vip=O2SCL_IX2(V,i,k);
839  double Viq=O2SCL_IX2(V,i,k+1);
840  //std::cout << "Vip,Viq: " << Vip << " " << Viq << std::endl;
841  O2SCL_IX2(V,i,k)=c*Vip-s*Viq;
842  O2SCL_IX2(V,i,k+1)=s*Vip+c*Viq;
843  }
844 
845  // [GSL] compute B <= B G
846 
847  {
848  double bk1=c*bk-s*z;
849 
850  double ap1=c*ap-s*bp;
851  double bp1=s*ap+c*bp;
852  double zp1=-s*aq;
853 
854  double aq1=c*aq;
855 
856  if (k>a) O2SCL_IX(f,k-1)=bk1;
857 
858  ak=ap1;
859  bk=bp1;
860  zk=zp1;
861 
862  ap=aq1;
863 
864  if (k<n-2) bp=O2SCL_IX(f,k+1);
865  else bp=0.0;
866 
867  y=ak;
868  z=zk;
869  }
870 
872 
873  // [GSL] Compute U <= U G
874 
875  for (i=0;i<M;i++) {
876  double Uip=O2SCL_IX2(U,i,k);
877  double Uiq=O2SCL_IX2(U,i,k+1);
878  //std::cout << "Uip2,Uiq2: " << Uip << " " << Uiq << std::endl;
879  O2SCL_IX2(U,i,k)=c*Uip-s*Uiq;
880  O2SCL_IX2(U,i,k+1)=s*Uip+c*Uiq;
881  }
882 
883  // [GSL] compute B <= G^T B
884 
885  //std::cout << "k,bk,ap2: " << k << " " << bk << " " << ap << std::endl;
886  //std::cout << "ak,zk,bp: " << ak << " " << zk << " "
887  //<< bp << std::endl;
888 
889  {
890  //std::cout << "prod1: " << c*ak << " " << s*zk << std::endl;
891  //std::cout << "prod2: " << c*bk << " " << s*ap << std::endl;
892  //std::cout << "prod3: " << s*bk << " " << c*ap << std::endl;
893  double ak1=c*ak-s*zk;
894  double bk1=c*bk-s*ap;
895  double zk1=-s*bp;
896 
897  double ap1=s*bk+c*ap;
898  double bp1=c*bp;
899 
900  O2SCL_IX(d,k)=ak1;
901 
902  ak=ak1;
903  bk=bk1;
904  zk=zk1;
905 
906  ap=ap1;
907  bp=bp1;
908  //std::cout << "c,s: " << c << " " << s << std::endl;
909  //std::cout << "k,bk,ap: " << k << " " << bk << " " << ap << std::endl;
910 
911  if (k < n-2) {
912  aq=O2SCL_IX(d,k+2);
913  } else {
914  aq=0.0;
915  }
916 
917  y=bk;
918  z=zk;
919  }
920  }
921 
922  O2SCL_IX(f,n-2)=bk;
923  O2SCL_IX(d,n-1)=ap;
924  //std::cout << "bk,ap: " << bk << " " << ap << std::endl;
925 
926  return;
927  }
928 
929 #ifdef DOXYGEN
930 }
931 #endif
void create_schur(double d0, double f0, double d1, double &c, double &s)
Desc.
Definition: svdstep_base.h:120
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:507
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
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:457
void create_givens(const double a, const double b, double &c, double &s)
Create a Givens rotation matrix.
double trailing_eigenvalue(size_t n, const vec_t &d, const vec2_t &f)
Desc.
Definition: svdstep_base.h:85
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:316
The namespace for linear algebra classes and functions.
Definition: bidiag.h:36
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:603
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:762
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
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:556

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