267 {
270 int idir;
271 int isub;
274 int ll;
275 int lll;
277 int maxit;
278 int oldll;
279 int oldm;
313 int maxitr;
314 bool matrixsplitflag;
315 bool iterflag;
320 bool rightside;
321 bool fwddir;
323 int mm1;
324 int mm0;
325 bool bchangedir;
326 int uend;
327 int cend;
328 int vend;
329
330
332 if( n==0 )
333 {
335 }
336 if( n==1 )
337 {
338 if( d(1)<0 )
339 {
340 d(1) = -d(1);
341 if( ncvt>0 )
342 {
344 }
345 }
347 }
348
349
350
351
362 maxitr = 12;
363 rightside = true;
364 fwddir = true;
365
366
367
368
370 for(
i=1;
i<=n-1;
i++)
371 {
373 }
375 for(
i=1;
i<=n-1;
i++)
376 {
378 }
379 e(n) = 0;
380 idir = 0;
381
382
383
384
387
388
389
390
391
392 if( !isupper )
393 {
394 for(
i=1;
i<=n-1;
i++)
395 {
402 }
403
404
405
406
407 if( nru>0 )
408 {
410 }
411 if( ncc>0 )
412 {
414 }
415 }
416
417
418
419
420
421
423 tol = tolmul*eps;
424 if( !isfractionalaccuracyrequired )
425 {
426 tol = -tol;
427 }
428
429
430
431
432 smax = 0;
434 {
436 }
437 for(
i=1;
i<=n-1;
i++)
438 {
440 }
441 sminl = 0;
442 if( tol>=0 )
443 {
444
445
446
447
449 if( sminoa!=0 )
450 {
453 {
456 if( sminoa==0 )
457 {
458 break;
459 }
460 }
461 }
464 }
465 else
466 {
467
468
469
470
472 }
473
474
475
476
477
478
479 maxit = maxitr*n*n;
481 oldll = -1;
482 oldm = -1;
483
484
485
486
488
489
490
491
492 while( true )
493 {
494
495
496
497
499 {
500 break;
501 }
503 {
506 }
507
508
509
510
512 {
514 }
516 smin = smax;
517 matrixsplitflag = false;
518 for(lll=1; lll<=
m-1; lll++)
519 {
523 if( tol<0 && abss<=thresh )
524 {
525 d(ll) = 0;
526 }
527 if( abse<=thresh )
528 {
529 matrixsplitflag = true;
530 break;
531 }
534 }
535 if( !matrixsplitflag )
536 {
537 ll = 0;
538 }
539 else
540 {
541
542
543
544
545 e(ll) = 0;
547 {
548
549
550
551
553 continue;
554 }
555 }
556 ll = ll+1;
557
558
559
560
562 {
563
564
565
566
571
572
573
574
575 if( ncvt>0 )
576 {
578 mm1 =
m-1+(vstart-1);
584 }
585 if( nru>0 )
586 {
594 }
595 if( ncc>0 )
596 {
604 }
606 continue;
607 }
608
609
610
611
612
613
614
615
616
617
618 bchangedir = false;
620 {
621 bchangedir = true;
622 }
624 {
625 bchangedir = true;
626 }
627 if( ll!=oldll ||
m!=oldm || bchangedir )
628 {
630 {
631
632
633
634
635 idir = 1;
636 }
637 else
638 {
639
640
641
642
643 idir = 2;
644 }
645 }
646
647
648
649
650 if( idir==1 )
651 {
652
653
654
655
656
658 {
660 continue;
661 }
662 if( tol>=0 )
663 {
664
665
666
667
668
671 iterflag = false;
672 for(lll=ll; lll<=
m-1; lll++)
673 {
675 {
676 e(lll) = 0;
677 iterflag = true;
678 break;
679 }
680 sminlo = sminl;
683 }
684 if( iterflag )
685 {
686 continue;
687 }
688 }
689 }
690 else
691 {
692
693
694
695
696
698 {
699 e(ll) = 0;
700 continue;
701 }
702 if( tol>=0 )
703 {
704
705
706
707
708
711 iterflag = false;
712 for(lll=
m-1; lll>=ll; lll--)
713 {
715 {
716 e(lll) = 0;
717 iterflag = true;
718 break;
719 }
720 sminlo = sminl;
723 }
724 if( iterflag )
725 {
726 continue;
727 }
728 }
729 }
730 oldll = ll;
732
733
734
735
736
738 {
739
740
741
742
743 shift = 0;
744 }
745 else
746 {
747
748
749
750
751 if( idir==1 )
752 {
755 }
756 else
757 {
760 }
761
762
763
764
765 if( sll>0 )
766 {
768 {
769 shift = 0;
770 }
771 }
772 }
773
774
775
776
778
779
780
781
782 if( shift==0 )
783 {
784 if( idir==1 )
785 {
786
787
788
789
790
791 cs = 1;
792 oldcs = 1;
793 for(
i=ll;
i<=
m-1;
i++)
794 {
797 {
799 }
804 work2(
i-ll+1) = oldcs;
805 work3(
i-ll+1) = oldsn;
806 }
810
811
812
813
814 if( ncvt>0 )
815 {
817 }
818 if( nru>0 )
819 {
821 }
822 if( ncc>0 )
823 {
825 }
826
827
828
829
831 {
833 }
834 }
835 else
836 {
837
838
839
840
841
842 cs = 1;
843 oldcs = 1;
844 for(
i=
m;
i>=ll+1;
i--)
845 {
848 {
850 }
856 work3(
i-ll) = -oldsn;
857 }
861
862
863
864
865 if( ncvt>0 )
866 {
868 }
869 if( nru>0 )
870 {
872 }
873 if( ncc>0 )
874 {
876 }
877
878
879
880
882 {
883 e(ll) = 0;
884 }
885 }
886 }
887 else
888 {
889
890
891
892
893 if( idir==1 )
894 {
895
896
897
898
899
902 for(
i=ll;
i<=
m-1;
i++)
903 {
906 {
908 }
909 f = cosr*d(
i)+sinr*e(
i);
910 e(
i) = cosr*e(
i)-sinr*d(
i);
912 d(
i+1) = cosr*d(
i+1);
915 f = cosl*e(
i)+sinl*d(
i+1);
916 d(
i+1) = cosl*d(
i+1)-sinl*e(
i);
918 {
920 e(
i+1) = cosl*e(
i+1);
921 }
922 work0(
i-ll+1) = cosr;
923 work1(
i-ll+1) = sinr;
924 work2(
i-ll+1) = cosl;
925 work3(
i-ll+1) = sinl;
926 }
928
929
930
931
932 if( ncvt>0 )
933 {
935 }
936 if( nru>0 )
937 {
939 }
940 if( ncc>0 )
941 {
943 }
944
945
946
947
949 {
951 }
952 }
953 else
954 {
955
956
957
958
959
962 for(
i=
m;
i>=ll+1;
i--)
963 {
966 {
968 }
969 f = cosr*d(
i)+sinr*e(
i-1);
970 e(
i-1) = cosr*e(
i-1)-sinr*d(
i);
972 d(
i-1) = cosr*d(
i-1);
975 f = cosl*e(
i-1)+sinl*d(
i-1);
976 d(
i-1) = cosl*d(
i-1)-sinl*e(
i-1);
978 {
980 e(
i-2) = cosl*e(
i-2);
981 }
986 }
988
989
990
991
993 {
994 e(ll) = 0;
995 }
996
997
998
999
1000 if( ncvt>0 )
1001 {
1003 }
1004 if( nru>0 )
1005 {
1007 }
1008 if( ncc>0 )
1009 {
1011 }
1012 }
1013 }
1014
1015
1016
1017
1018 continue;
1019 }
1020
1021
1022
1023
1025 {
1027 {
1029
1030
1031
1032
1033 if( ncvt>0 )
1034 {
1036 }
1037 }
1038 }
1039
1040
1041
1042
1043
1044 for(
i=1;
i<=n-1;
i++)
1045 {
1046
1047
1048
1049
1050 isub = 1;
1051 smin = d(1);
1052 for(
j=2;
j<=n+1-
i;
j++)
1053 {
1055 {
1058 }
1059 }
1061 {
1062
1063
1064
1065
1068 if( ncvt>0 )
1069 {
1074 }
1075 if( nru>0 )
1076 {
1081 }
1082 if( ncc>0 )
1083 {
1088 }
1089 }
1090 }
1092 }
static const ampf getAlgoPascalEpsilon()
static const ampf getAlgoPascalMinNumber()
raw_vector< T > getvector(int iStart, int iEnd)
void setbounds(int iLow, int iHigh)
raw_vector< T > getrow(int iRow, int iColumnStart, int iColumnEnd)
raw_vector< T > getcolumn(int iColumn, int iRowStart, int iRowEnd)
static matrix mu(matrix A, const ring R)
const ampf< Precision > abs(const ampf< Precision > &x)
const ampf< Precision > maximum(const ampf< Precision > &x, const ampf< Precision > &y)
const ampf< Precision > sqrt(const ampf< Precision > &x)
const ampf< Precision > pow(const ampf< Precision > &x, const ampf< Precision > &y)
const ampf< Precision > sqr(const ampf< Precision > &x)
const ampf< Precision > minimum(const ampf< Precision > &x, const ampf< Precision > &y)
int maxint(int m1, int m2)
void vmove(raw_vector< T > vdst, const_raw_vector< T > vsrc)
void vadd(raw_vector< T > vdst, const_raw_vector< T > vsrc)
void vmul(raw_vector< T > vdst, T2 alpha)
void vsub(raw_vector< T > vdst, const_raw_vector< T > vsrc)
amp::ampf< Precision > extsignbdsqr(amp::ampf< Precision > a, amp::ampf< Precision > b)
void svdv2x2(amp::ampf< Precision > f, amp::ampf< Precision > g, amp::ampf< Precision > h, amp::ampf< Precision > &ssmin, amp::ampf< Precision > &ssmax, amp::ampf< Precision > &snr, amp::ampf< Precision > &csr, amp::ampf< Precision > &snl, amp::ampf< Precision > &csl)
void svd2x2(amp::ampf< Precision > f, amp::ampf< Precision > g, amp::ampf< Precision > h, amp::ampf< Precision > &ssmin, amp::ampf< Precision > &ssmax)
void applyrotationsfromtheright(bool isforward, int m1, int m2, int n1, int n2, const ap::template_1d_array< amp::ampf< Precision > > &c, const ap::template_1d_array< amp::ampf< Precision > > &s, ap::template_2d_array< amp::ampf< Precision > > &a, ap::template_1d_array< amp::ampf< Precision > > &work)
void generaterotation(amp::ampf< Precision > f, amp::ampf< Precision > g, amp::ampf< Precision > &cs, amp::ampf< Precision > &sn, amp::ampf< Precision > &r)
void applyrotationsfromtheleft(bool isforward, int m1, int m2, int n1, int n2, const ap::template_1d_array< amp::ampf< Precision > > &c, const ap::template_1d_array< amp::ampf< Precision > > &s, ap::template_2d_array< amp::ampf< Precision > > &a, ap::template_1d_array< amp::ampf< Precision > > &work)