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) {
178 double c,s,a11,a12,a21,a22;
180 double d0=O2SCL_IX(d,0);
181 double f0=O2SCL_IX(f,0);
183 double d1=O2SCL_IX(d,1);
193 O2SCL_IX(d,0)=c*f0-s*d1;
194 O2SCL_IX(f,0)=s*f0+c*d1;
199 for (
size_t i=0;i<M;i++) {
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;
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;
218 }
else if (d1 == 0.0) {
226 O2SCL_IX(d,0)=d0*c-f0*s;
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;
255 for (
size_t i=0;i<N;i++) {
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;
266 if (hypot(a11,a21) < hypot(a12,a22)) {
272 t1=a11; a11=a12; a12=t1;
273 t2=a21; a21=a22; a22=t2;
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;
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;
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;
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) {
318 double c,s,a11,a12,a21,a22;
320 double d0=O2SCL_IX(d,a);
321 double f0=O2SCL_IX(f,a);
323 double d1=O2SCL_IX(d,a+1);
333 O2SCL_IX(d,a)=c*f0-s*d1;
334 O2SCL_IX(f,a)=s*f0+c*d1;
339 for (
size_t i=0;i<M;i++) {
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;
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;
358 }
else if (d1 == 0.0) {
366 O2SCL_IX(d,a)=d0*c-f0*s;
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;
395 for (
size_t i=0;i<N;i++) {
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;
406 if (hypot(a11,a21)<hypot(a12,a22)) {
412 t1=a11; a11=a12; a12=t1;
413 t2=a21; a21=a22; a22=t2;
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;
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;
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;
454 template<
class vec_t,
class vec2_t,
class mat_t>
456 vec2_t &f, mat_t &U,
size_t k0) {
460 double x=O2SCL_IX(f,k0);
461 double y=O2SCL_IX(d,k0+1);
463 for (
size_t k=k0;k<n-1;k++) {
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);
472 O2SCL_IX2(U,i,k0)=c*Uip-s*Uiq;
473 O2SCL_IX2(U,i,k+1)=s*Uip+c*Uiq;
478 O2SCL_IX(d,k+1)=s*x+c*y;
481 O2SCL_IX(f,k)=c*x-s*y;
485 double z=O2SCL_IX(f,k+1);
504 template<
class vec_t,
class vec2_t,
class mat_t>
506 vec2_t &f, mat_t &V) {
510 double x=O2SCL_IX(d,n-2);
511 double y=O2SCL_IX(f,n-2);
513 for (
size_t k=n-1;k-- > 0;) {
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;
528 O2SCL_IX(d,k)=c*x-s*y;
531 O2SCL_IX(f,k)=s*x+c*y;
535 double z=O2SCL_IX(f,k-1);
553 template<
class vec_t,
class vec2_t,
class mat_t>
555 vec2_t &f, mat_t &V,
size_t a) {
559 double x=O2SCL_IX(d,n-2);
560 double y=O2SCL_IX(f,n-2);
562 for (
int k=((
int)n)-1;k>=((int)a);k--) {
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;
577 O2SCL_IX(d,k)=c*x-s*y;
580 O2SCL_IX(f,k)=s*x+c*y;
584 double z=O2SCL_IX(f,k-1);
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) {
605 double ak, bk, zk, ap, bp, aq, bq;
622 for (i=0;i<n-1;i++) {
623 double d_i=O2SCL_IX(d,i);
631 double d_nm1=O2SCL_IX(d,n-1);
640 double d0=O2SCL_IX(d,0);
641 double f0=O2SCL_IX(f,0);
643 double d1=O2SCL_IX(d,1);
644 double f1=O2SCL_IX(f,1);
662 for (k=0; k < n-1; k++) {
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;
681 double ap1=c*ap-s*bp;
682 double bp1=s*ap+c*bp;
687 if (k > 0) O2SCL_IX(f,k-1)=bk1;
695 if (k<n-2) bp=O2SCL_IX(f,k+1);
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;
715 double ak1=c*ak-s*zk;
716 double bk1=c*bk-s*ap;
719 double ap1=s*bk+c*ap;
759 template<
class vec_t,
class vec2_t,
class mat_t,
class mat2_t>
761 vec_t &d, vec2_t &f, mat_t &U, mat2_t &V,
size_t a) {
764 double ak, bk, zk, ap, bp, aq, bq;
783 for (i=a;i<n-1;i++) {
784 double d_i=O2SCL_IX(d,i);
794 double d_nm1=O2SCL_IX(d,n-1);
803 double d0=O2SCL_IX(d,a);
804 double f0=O2SCL_IX(f,a);
806 double d1=O2SCL_IX(d,a+1);
807 double f1=O2SCL_IX(f,a+1);
828 for (k=a;k<n-1;k++) {
836 double Vip=O2SCL_IX2(V,i,k);
837 double Viq=O2SCL_IX2(V,i,k+1);
839 O2SCL_IX2(V,i,k)=c*Vip-s*Viq;
840 O2SCL_IX2(V,i,k+1)=s*Vip+c*Viq;
848 double ap1=c*ap-s*bp;
849 double bp1=s*ap+c*bp;
854 if (k>a) O2SCL_IX(f,k-1)=bk1;
862 if (k<n-2) bp=O2SCL_IX(f,k+1);
874 double Uip=O2SCL_IX2(U,i,k);
875 double Uiq=O2SCL_IX2(U,i,k+1);
877 O2SCL_IX2(U,i,k)=c*Uip-s*Uiq;
878 O2SCL_IX2(U,i,k+1)=s*Uip+c*Uiq;
891 double ak1=c*ak-s*zk;
892 double bk1=c*bk-s*ap;
895 double ap1=s*bk+c*ap;