57 template<
class vec_t,
class vec2_t>
60 double d_i=O2SCL_IX(d,0);
62 for (
size_t i=0;i<N-1;i++) {
64 double f_i=O2SCL_IX(f,i);
65 double d_ip1=O2SCL_IX(d,i+1);
67 double dbl_eps=std::numeric_limits<double>::epsilon();
69 if (fabs (f_i)<dbl_eps*(fabs(d_i)+fabs(d_ip1))) {
84 template<
class vec_t,
class vec2_t>
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);
92 double ta=da*da+fa*fa;
93 double tb=db*db+fb*fb;
96 double dt=(ta-tb)/2.0;
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);
108 mu=(r1>0) ? P / r1 : 0.0;
123 double apq=2.0*d0*f0;
125 if (d0 == 0 || f0 == 0) {
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) {
141 scale=ldexp(1.0,-(d0_exp+f0_exp)/4);
150 double tau=(f0*f0+(d1+d0)*(d1-d0))/apq;
153 t=1.0/(tau+hypot(1.0,tau));
155 t=-1.0/(-tau+hypot(1.0,tau));
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) {
179 double c,s,a11,a12,a21,a22;
181 double d0=O2SCL_IX(d,0);
182 double f0=O2SCL_IX(f,0);
184 double d1=O2SCL_IX(d,1);
194 O2SCL_IX(d,0)=c*f0-s*d1;
195 O2SCL_IX(f,0)=s*f0+c*d1;
200 for (
size_t i=0;i<M;i++) {
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;
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;
219 }
else if (d1 == 0.0) {
227 O2SCL_IX(d,0)=d0*c-f0*s;
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;
256 for (
size_t i=0;i<N;i++) {
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;
267 if (hypot(a11,a21) < hypot(a12,a22)) {
273 t1=a11; a11=a12; a12=t1;
274 t2=a21; a21=a22; a22=t2;
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;
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;
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;
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) {
320 double c,s,a11,a12,a21,a22;
322 double d0=O2SCL_IX(d,a);
323 double f0=O2SCL_IX(f,a);
325 double d1=O2SCL_IX(d,a+1);
335 O2SCL_IX(d,a)=c*f0-s*d1;
336 O2SCL_IX(f,a)=s*f0+c*d1;
341 for (
size_t i=0;i<M;i++) {
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;
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;
360 }
else if (d1 == 0.0) {
368 O2SCL_IX(d,a)=d0*c-f0*s;
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;
397 for (
size_t i=0;i<N;i++) {
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;
408 if (hypot(a11,a21)<hypot(a12,a22)) {
414 t1=a11; a11=a12; a12=t1;
415 t2=a21; a21=a22; a22=t2;
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;
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;
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;
456 template<
class vec_t,
class vec2_t,
class mat_t>
458 vec2_t &f, mat_t &U,
size_t k0) {
462 double x=O2SCL_IX(f,k0);
463 double y=O2SCL_IX(d,k0+1);
465 for (
size_t k=k0;k<n-1;k++) {
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);
474 O2SCL_IX2(U,i,k0)=c*Uip-s*Uiq;
475 O2SCL_IX2(U,i,k+1)=s*Uip+c*Uiq;
480 O2SCL_IX(d,k+1)=s*x+c*y;
483 O2SCL_IX(f,k)=c*x-s*y;
487 double z=O2SCL_IX(f,k+1);
506 template<
class vec_t,
class vec2_t,
class mat_t>
508 vec2_t &f, mat_t &V) {
512 double x=O2SCL_IX(d,n-2);
513 double y=O2SCL_IX(f,n-2);
515 for (
size_t k=n-1;k-- > 0;) {
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;
530 O2SCL_IX(d,k)=c*x-s*y;
533 O2SCL_IX(f,k)=s*x+c*y;
537 double z=O2SCL_IX(f,k-1);
555 template<
class vec_t,
class vec2_t,
class mat_t>
557 vec2_t &f, mat_t &V,
size_t a) {
561 double x=O2SCL_IX(d,n-2);
562 double y=O2SCL_IX(f,n-2);
564 for (
int k=((
int)n)-1;k>=((int)a);k--) {
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;
579 O2SCL_IX(d,k)=c*x-s*y;
582 O2SCL_IX(f,k)=s*x+c*y;
586 double z=O2SCL_IX(f,k-1);
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) {
607 double ak, bk, zk, ap, bp, aq, bq;
624 for (i=0;i<n-1;i++) {
625 double d_i=O2SCL_IX(d,i);
633 double d_nm1=O2SCL_IX(d,n-1);
642 double d0=O2SCL_IX(d,0);
643 double f0=O2SCL_IX(f,0);
645 double d1=O2SCL_IX(d,1);
646 double f1=O2SCL_IX(f,1);
664 for (k=0; k < n-1; k++) {
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;
683 double ap1=c*ap-s*bp;
684 double bp1=s*ap+c*bp;
689 if (k > 0) O2SCL_IX(f,k-1)=bk1;
697 if (k<n-2) bp=O2SCL_IX(f,k+1);
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;
717 double ak1=c*ak-s*zk;
718 double bk1=c*bk-s*ap;
721 double ap1=s*bk+c*ap;
761 template<
class vec_t,
class vec2_t,
class mat_t,
class mat2_t>
763 vec_t &d, vec2_t &f, mat_t &U, mat2_t &V,
size_t a) {
766 double ak, bk, zk, ap, bp, aq, bq;
785 for (i=a;i<n-1;i++) {
786 double d_i=O2SCL_IX(d,i);
796 double d_nm1=O2SCL_IX(d,n-1);
805 double d0=O2SCL_IX(d,a);
806 double f0=O2SCL_IX(f,a);
808 double d1=O2SCL_IX(d,a+1);
809 double f1=O2SCL_IX(f,a+1);
830 for (k=a;k<n-1;k++) {
838 double Vip=O2SCL_IX2(V,i,k);
839 double Viq=O2SCL_IX2(V,i,k+1);
841 O2SCL_IX2(V,i,k)=c*Vip-s*Viq;
842 O2SCL_IX2(V,i,k+1)=s*Vip+c*Viq;
850 double ap1=c*ap-s*bp;
851 double bp1=s*ap+c*bp;
856 if (k>a) O2SCL_IX(f,k-1)=bk1;
864 if (k<n-2) bp=O2SCL_IX(f,k+1);
876 double Uip=O2SCL_IX2(U,i,k);
877 double Uiq=O2SCL_IX2(U,i,k+1);
879 O2SCL_IX2(U,i,k)=c*Uip-s*Uiq;
880 O2SCL_IX2(U,i,k+1)=s*Uip+c*Uiq;
893 double ak1=c*ak-s*zk;
894 double bk1=c*bk-s*ap;
897 double ap1=s*bk+c*ap;
void create_schur(double d0, double f0, double d1, double &c, double &s)
Desc.
void chase_out_trailing_zero(size_t N, size_t n, vec_t &d, vec2_t &f, mat_t &V)
Desc.
void svd2(size_t M, size_t N, vec_t &d, vec2_t &f, mat_t &U, mat2_t &V)
2-variable SVD
void chase_out_intermediate_zero(size_t M, size_t n, vec_t &d, vec2_t &f, mat_t &U, size_t k0)
Desc.
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.
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.
The namespace for linear algebra classes and functions.
void qrstep(size_t M, size_t N, size_t n, vec_t &d, vec2_t &f, mat_t &U, mat2_t &V)
Desc.
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()
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.
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.