Mercurial > jwe > octave

comparison liboctave/array/CMatrix.cc @ 21136:7cac4e7458f2

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: clean up code around calls to current_liboctave_error_handler. Remove statements after call to handler that are no longer reachable. Place input validation first and immediately call handler if necessary. Change if/error_handler/else to if/error_handler and re-indent code. * Array-util.cc, Array.cc, CColVector.cc, CDiagMatrix.cc, CMatrix.cc, CNDArray.cc, CRowVector.cc, CSparse.cc, DiagArray2.cc, MArray.cc, PermMatrix.cc, Sparse.cc, Sparse.h, chMatrix.cc, chNDArray.cc, dColVector.cc, dDiagMatrix.cc, dMatrix.cc, dNDArray.cc, dRowVector.cc, dSparse.cc, fCColVector.cc, fCDiagMatrix.cc, fCMatrix.cc, fCNDArray.cc, fCRowVector.cc, fColVector.cc, fDiagMatrix.cc, fMatrix.cc, fNDArray.cc, fRowVector.cc, idx-vector.cc, CmplxAEPBAL.cc, CmplxCHOL.cc, CmplxGEPBAL.cc, CmplxHESS.cc, CmplxLU.cc, CmplxQR.cc, CmplxSCHUR.cc, CmplxSVD.cc, DASPK.cc, EIG.cc, LSODE.cc, Quad.cc, SparseCmplxCHOL.cc, SparseCmplxLU.cc, SparseCmplxQR.cc, SparseQR.cc, SparsedbleCHOL.cc, SparsedbleLU.cc, base-lu.cc, bsxfun-defs.cc, dbleAEPBAL.cc, dbleCHOL.cc, dbleGEPBAL.cc, dbleHESS.cc, dbleLU.cc, dbleQR.cc, dbleSCHUR.cc, dbleSVD.cc, eigs-base.cc, fCmplxAEPBAL.cc, fCmplxCHOL.cc, fCmplxLU.cc, fCmplxQR.cc, fCmplxSCHUR.cc, fEIG.cc, floatAEPBAL.cc, floatCHOL.cc, floatGEPBAL.cc, floatHESS.cc, floatLU.cc, floatQR.cc, floatSCHUR.cc, floatSVD.cc, lo-specfun.cc, oct-fftw.cc, oct-rand.cc, oct-spparms.cc, sparse-base-chol.cc, sparse-dmsolve.cc, file-ops.cc, lo-sysdep.cc, mach-info.cc, oct-env.cc, oct-syscalls.cc, cmd-edit.cc, cmd-hist.cc, data-conv.cc, lo-ieee.cc, lo-regexp.cc, oct-base64.cc, oct-shlib.cc, pathsearch.cc, singleton-cleanup.cc, sparse-util.cc, unwind-prot.cc: Remove statements after call to handler that are no longer reachable. Place input validation first and immediately call handler if necessary. Change if/error_handler/else to if/error_handler and re-indent code.
author Rik <rik@octave.org>
date Sat, 23 Jan 2016 13:52:03 -0800
parents 228b65504557
children 26f85aa072de
comparison
equal deleted inserted replaced
21135:95da3bc8a281 21136:7cac4e7458f2
401 { 401 {
402 octave_idx_type a_nr = a.rows (); 402 octave_idx_type a_nr = a.rows ();
403 octave_idx_type a_nc = a.cols (); 403 octave_idx_type a_nc = a.cols ();
404 404
405 if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ()) 405 if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ())
406 { 406 (*current_liboctave_error_handler) ("range error for insert");
407 (*current_liboctave_error_handler) ("range error for insert");
408 return *this;
409 }
410 407
411 if (a_nr >0 && a_nc > 0) 408 if (a_nr >0 && a_nc > 0)
412 { 409 {
413 make_unique (); 410 make_unique ();
414 411
424 ComplexMatrix::insert (const RowVector& a, octave_idx_type r, octave_idx_type c) 421 ComplexMatrix::insert (const RowVector& a, octave_idx_type r, octave_idx_type c)
425 { 422 {
426 octave_idx_type a_len = a.numel (); 423 octave_idx_type a_len = a.numel ();
427 424
428 if (r < 0 || r >= rows () || c < 0 || c + a_len > cols ()) 425 if (r < 0 || r >= rows () || c < 0 || c + a_len > cols ())
429 { 426 (*current_liboctave_error_handler) ("range error for insert");
430 (*current_liboctave_error_handler) ("range error for insert");
431 return *this;
432 }
433 427
434 if (a_len > 0) 428 if (a_len > 0)
435 { 429 {
436 make_unique (); 430 make_unique ();
437 431
447 octave_idx_type r, octave_idx_type c) 441 octave_idx_type r, octave_idx_type c)
448 { 442 {
449 octave_idx_type a_len = a.numel (); 443 octave_idx_type a_len = a.numel ();
450 444
451 if (r < 0 || r + a_len > rows () || c < 0 || c >= cols ()) 445 if (r < 0 || r + a_len > rows () || c < 0 || c >= cols ())
452 { 446 (*current_liboctave_error_handler) ("range error for insert");
453 (*current_liboctave_error_handler) ("range error for insert");
454 return *this;
455 }
456 447
457 if (a_len > 0) 448 if (a_len > 0)
458 { 449 {
459 make_unique (); 450 make_unique ();
460 451
471 { 462 {
472 octave_idx_type a_nr = a.rows (); 463 octave_idx_type a_nr = a.rows ();
473 octave_idx_type a_nc = a.cols (); 464 octave_idx_type a_nc = a.cols ();
474 465
475 if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ()) 466 if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ())
476 { 467 (*current_liboctave_error_handler) ("range error for insert");
477 (*current_liboctave_error_handler) ("range error for insert");
478 return *this;
479 }
480 468
481 fill (0.0, r, c, r + a_nr - 1, c + a_nc - 1); 469 fill (0.0, r, c, r + a_nr - 1, c + a_nc - 1);
482 470
483 octave_idx_type a_len = a.length (); 471 octave_idx_type a_len = a.length ();
484 472
505 ComplexMatrix::insert (const ComplexRowVector& a, 493 ComplexMatrix::insert (const ComplexRowVector& a,
506 octave_idx_type r, octave_idx_type c) 494 octave_idx_type r, octave_idx_type c)
507 { 495 {
508 octave_idx_type a_len = a.numel (); 496 octave_idx_type a_len = a.numel ();
509 if (r < 0 || r >= rows () || c < 0 || c + a_len > cols ()) 497 if (r < 0 || r >= rows () || c < 0 || c + a_len > cols ())
510 { 498 (*current_liboctave_error_handler) ("range error for insert");
511 (*current_liboctave_error_handler) ("range error for insert");
512 return *this;
513 }
514 499
515 for (octave_idx_type i = 0; i < a_len; i++) 500 for (octave_idx_type i = 0; i < a_len; i++)
516 elem (r, c+i) = a.elem (i); 501 elem (r, c+i) = a.elem (i);
517 502
518 return *this; 503 return *this;
523 octave_idx_type r, octave_idx_type c) 508 octave_idx_type r, octave_idx_type c)
524 { 509 {
525 octave_idx_type a_len = a.numel (); 510 octave_idx_type a_len = a.numel ();
526 511
527 if (r < 0 || r + a_len > rows () || c < 0 || c >= cols ()) 512 if (r < 0 || r + a_len > rows () || c < 0 || c >= cols ())
528 { 513 (*current_liboctave_error_handler) ("range error for insert");
529 (*current_liboctave_error_handler) ("range error for insert");
530 return *this;
531 }
532 514
533 if (a_len > 0) 515 if (a_len > 0)
534 { 516 {
535 make_unique (); 517 make_unique ();
536 518
547 { 529 {
548 octave_idx_type a_nr = a.rows (); 530 octave_idx_type a_nr = a.rows ();
549 octave_idx_type a_nc = a.cols (); 531 octave_idx_type a_nc = a.cols ();
550 532
551 if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ()) 533 if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ())
552 { 534 (*current_liboctave_error_handler) ("range error for insert");
553 (*current_liboctave_error_handler) ("range error for insert");
554 return *this;
555 }
556 535
557 fill (0.0, r, c, r + a_nr - 1, c + a_nc - 1); 536 fill (0.0, r, c, r + a_nr - 1, c + a_nc - 1);
558 537
559 octave_idx_type a_len = a.length (); 538 octave_idx_type a_len = a.length ();
560 539
612 octave_idx_type nr = rows (); 591 octave_idx_type nr = rows ();
613 octave_idx_type nc = cols (); 592 octave_idx_type nc = cols ();
614 593
615 if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0 594 if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0
616 || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc) 595 || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc)
617 { 596 (*current_liboctave_error_handler) ("range error for fill");
618 (*current_liboctave_error_handler) ("range error for fill");
619 return *this;
620 }
621 597
622 if (r1 > r2) { std::swap (r1, r2); } 598 if (r1 > r2) { std::swap (r1, r2); }
623 if (c1 > c2) { std::swap (c1, c2); } 599 if (c1 > c2) { std::swap (c1, c2); }
624 600
625 if (r2 >= r1 && c2 >= c1) 601 if (r2 >= r1 && c2 >= c1)
641 octave_idx_type nr = rows (); 617 octave_idx_type nr = rows ();
642 octave_idx_type nc = cols (); 618 octave_idx_type nc = cols ();
643 619
644 if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0 620 if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0
645 || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc) 621 || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc)
646 { 622 (*current_liboctave_error_handler) ("range error for fill");
647 (*current_liboctave_error_handler) ("range error for fill");
648 return *this;
649 }
650 623
651 if (r1 > r2) { std::swap (r1, r2); } 624 if (r1 > r2) { std::swap (r1, r2); }
652 if (c1 > c2) { std::swap (c1, c2); } 625 if (c1 > c2) { std::swap (c1, c2); }
653 626
654 if (r2 >= r1 && c2 >=c1) 627 if (r2 >= r1 && c2 >=c1)
667 ComplexMatrix::append (const Matrix& a) const 640 ComplexMatrix::append (const Matrix& a) const
668 { 641 {
669 octave_idx_type nr = rows (); 642 octave_idx_type nr = rows ();
670 octave_idx_type nc = cols (); 643 octave_idx_type nc = cols ();
671 if (nr != a.rows ()) 644 if (nr != a.rows ())
672 { 645 (*current_liboctave_error_handler) ("row dimension mismatch for append");
673 (*current_liboctave_error_handler) ("row dimension mismatch for append");
674 return *this;
675 }
676 646
677 octave_idx_type nc_insert = nc; 647 octave_idx_type nc_insert = nc;
678 ComplexMatrix retval (nr, nc + a.cols ()); 648 ComplexMatrix retval (nr, nc + a.cols ());
679 retval.insert (*this, 0, 0); 649 retval.insert (*this, 0, 0);
680 retval.insert (a, 0, nc_insert); 650 retval.insert (a, 0, nc_insert);
685 ComplexMatrix::append (const RowVector& a) const 655 ComplexMatrix::append (const RowVector& a) const
686 { 656 {
687 octave_idx_type nr = rows (); 657 octave_idx_type nr = rows ();
688 octave_idx_type nc = cols (); 658 octave_idx_type nc = cols ();
689 if (nr != 1) 659 if (nr != 1)
690 { 660 (*current_liboctave_error_handler) ("row dimension mismatch for append");
691 (*current_liboctave_error_handler) ("row dimension mismatch for append");
692 return *this;
693 }
694 661
695 octave_idx_type nc_insert = nc; 662 octave_idx_type nc_insert = nc;
696 ComplexMatrix retval (nr, nc + a.numel ()); 663 ComplexMatrix retval (nr, nc + a.numel ());
697 retval.insert (*this, 0, 0); 664 retval.insert (*this, 0, 0);
698 retval.insert (a, 0, nc_insert); 665 retval.insert (a, 0, nc_insert);
703 ComplexMatrix::append (const ColumnVector& a) const 670 ComplexMatrix::append (const ColumnVector& a) const
704 { 671 {
705 octave_idx_type nr = rows (); 672 octave_idx_type nr = rows ();
706 octave_idx_type nc = cols (); 673 octave_idx_type nc = cols ();
707 if (nr != a.numel ()) 674 if (nr != a.numel ())
708 { 675 (*current_liboctave_error_handler) ("row dimension mismatch for append");
709 (*current_liboctave_error_handler) ("row dimension mismatch for append");
710 return *this;
711 }
712 676
713 octave_idx_type nc_insert = nc; 677 octave_idx_type nc_insert = nc;
714 ComplexMatrix retval (nr, nc + 1); 678 ComplexMatrix retval (nr, nc + 1);
715 retval.insert (*this, 0, 0); 679 retval.insert (*this, 0, 0);
716 retval.insert (a, 0, nc_insert); 680 retval.insert (a, 0, nc_insert);
721 ComplexMatrix::append (const DiagMatrix& a) const 685 ComplexMatrix::append (const DiagMatrix& a) const
722 { 686 {
723 octave_idx_type nr = rows (); 687 octave_idx_type nr = rows ();
724 octave_idx_type nc = cols (); 688 octave_idx_type nc = cols ();
725 if (nr != a.rows ()) 689 if (nr != a.rows ())
726 { 690 (*current_liboctave_error_handler) ("row dimension mismatch for append");
727 (*current_liboctave_error_handler) ("row dimension mismatch for append");
728 return *this;
729 }
730 691
731 octave_idx_type nc_insert = nc; 692 octave_idx_type nc_insert = nc;
732 ComplexMatrix retval (nr, nc + a.cols ()); 693 ComplexMatrix retval (nr, nc + a.cols ());
733 retval.insert (*this, 0, 0); 694 retval.insert (*this, 0, 0);
734 retval.insert (a, 0, nc_insert); 695 retval.insert (a, 0, nc_insert);
739 ComplexMatrix::append (const ComplexMatrix& a) const 700 ComplexMatrix::append (const ComplexMatrix& a) const
740 { 701 {
741 octave_idx_type nr = rows (); 702 octave_idx_type nr = rows ();
742 octave_idx_type nc = cols (); 703 octave_idx_type nc = cols ();
743 if (nr != a.rows ()) 704 if (nr != a.rows ())
744 { 705 (*current_liboctave_error_handler) ("row dimension mismatch for append");
745 (*current_liboctave_error_handler) ("row dimension mismatch for append");
746 return *this;
747 }
748 706
749 octave_idx_type nc_insert = nc; 707 octave_idx_type nc_insert = nc;
750 ComplexMatrix retval (nr, nc + a.cols ()); 708 ComplexMatrix retval (nr, nc + a.cols ());
751 retval.insert (*this, 0, 0); 709 retval.insert (*this, 0, 0);
752 retval.insert (a, 0, nc_insert); 710 retval.insert (a, 0, nc_insert);
757 ComplexMatrix::append (const ComplexRowVector& a) const 715 ComplexMatrix::append (const ComplexRowVector& a) const
758 { 716 {
759 octave_idx_type nr = rows (); 717 octave_idx_type nr = rows ();
760 octave_idx_type nc = cols (); 718 octave_idx_type nc = cols ();
761 if (nr != 1) 719 if (nr != 1)
762 { 720 (*current_liboctave_error_handler) ("row dimension mismatch for append");
763 (*current_liboctave_error_handler) ("row dimension mismatch for append");
764 return *this;
765 }
766 721
767 octave_idx_type nc_insert = nc; 722 octave_idx_type nc_insert = nc;
768 ComplexMatrix retval (nr, nc + a.numel ()); 723 ComplexMatrix retval (nr, nc + a.numel ());
769 retval.insert (*this, 0, 0); 724 retval.insert (*this, 0, 0);
770 retval.insert (a, 0, nc_insert); 725 retval.insert (a, 0, nc_insert);
775 ComplexMatrix::append (const ComplexColumnVector& a) const 730 ComplexMatrix::append (const ComplexColumnVector& a) const
776 { 731 {
777 octave_idx_type nr = rows (); 732 octave_idx_type nr = rows ();
778 octave_idx_type nc = cols (); 733 octave_idx_type nc = cols ();
779 if (nr != a.numel ()) 734 if (nr != a.numel ())
780 { 735 (*current_liboctave_error_handler) ("row dimension mismatch for append");
781 (*current_liboctave_error_handler) ("row dimension mismatch for append");
782 return *this;
783 }
784 736
785 octave_idx_type nc_insert = nc; 737 octave_idx_type nc_insert = nc;
786 ComplexMatrix retval (nr, nc + 1); 738 ComplexMatrix retval (nr, nc + 1);
787 retval.insert (*this, 0, 0); 739 retval.insert (*this, 0, 0);
788 retval.insert (a, 0, nc_insert); 740 retval.insert (a, 0, nc_insert);
793 ComplexMatrix::append (const ComplexDiagMatrix& a) const 745 ComplexMatrix::append (const ComplexDiagMatrix& a) const
794 { 746 {
795 octave_idx_type nr = rows (); 747 octave_idx_type nr = rows ();
796 octave_idx_type nc = cols (); 748 octave_idx_type nc = cols ();
797 if (nr != a.rows ()) 749 if (nr != a.rows ())
798 { 750 (*current_liboctave_error_handler) ("row dimension mismatch for append");
799 (*current_liboctave_error_handler) ("row dimension mismatch for append");
800 return *this;
801 }
802 751
803 octave_idx_type nc_insert = nc; 752 octave_idx_type nc_insert = nc;
804 ComplexMatrix retval (nr, nc + a.cols ()); 753 ComplexMatrix retval (nr, nc + a.cols ());
805 retval.insert (*this, 0, 0); 754 retval.insert (*this, 0, 0);
806 retval.insert (a, 0, nc_insert); 755 retval.insert (a, 0, nc_insert);
811 ComplexMatrix::stack (const Matrix& a) const 760 ComplexMatrix::stack (const Matrix& a) const
812 { 761 {
813 octave_idx_type nr = rows (); 762 octave_idx_type nr = rows ();
814 octave_idx_type nc = cols (); 763 octave_idx_type nc = cols ();
815 if (nc != a.cols ()) 764 if (nc != a.cols ())
816 { 765 (*current_liboctave_error_handler) ("column dimension mismatch for stack");
817 (*current_liboctave_error_handler)
818 ("column dimension mismatch for stack");
819 return *this;
820 }
821 766
822 octave_idx_type nr_insert = nr; 767 octave_idx_type nr_insert = nr;
823 ComplexMatrix retval (nr + a.rows (), nc); 768 ComplexMatrix retval (nr + a.rows (), nc);
824 retval.insert (*this, 0, 0); 769 retval.insert (*this, 0, 0);
825 retval.insert (a, nr_insert, 0); 770 retval.insert (a, nr_insert, 0);
830 ComplexMatrix::stack (const RowVector& a) const 775 ComplexMatrix::stack (const RowVector& a) const
831 { 776 {
832 octave_idx_type nr = rows (); 777 octave_idx_type nr = rows ();
833 octave_idx_type nc = cols (); 778 octave_idx_type nc = cols ();
834 if (nc != a.numel ()) 779 if (nc != a.numel ())
835 { 780 (*current_liboctave_error_handler) ("column dimension mismatch for stack");
836 (*current_liboctave_error_handler)
837 ("column dimension mismatch for stack");
838 return *this;
839 }
840 781
841 octave_idx_type nr_insert = nr; 782 octave_idx_type nr_insert = nr;
842 ComplexMatrix retval (nr + 1, nc); 783 ComplexMatrix retval (nr + 1, nc);
843 retval.insert (*this, 0, 0); 784 retval.insert (*this, 0, 0);
844 retval.insert (a, nr_insert, 0); 785 retval.insert (a, nr_insert, 0);
849 ComplexMatrix::stack (const ColumnVector& a) const 790 ComplexMatrix::stack (const ColumnVector& a) const
850 { 791 {
851 octave_idx_type nr = rows (); 792 octave_idx_type nr = rows ();
852 octave_idx_type nc = cols (); 793 octave_idx_type nc = cols ();
853 if (nc != 1) 794 if (nc != 1)
854 { 795 (*current_liboctave_error_handler) ("column dimension mismatch for stack");
855 (*current_liboctave_error_handler)
856 ("column dimension mismatch for stack");
857 return *this;
858 }
859 796
860 octave_idx_type nr_insert = nr; 797 octave_idx_type nr_insert = nr;
861 ComplexMatrix retval (nr + a.numel (), nc); 798 ComplexMatrix retval (nr + a.numel (), nc);
862 retval.insert (*this, 0, 0); 799 retval.insert (*this, 0, 0);
863 retval.insert (a, nr_insert, 0); 800 retval.insert (a, nr_insert, 0);
868 ComplexMatrix::stack (const DiagMatrix& a) const 805 ComplexMatrix::stack (const DiagMatrix& a) const
869 { 806 {
870 octave_idx_type nr = rows (); 807 octave_idx_type nr = rows ();
871 octave_idx_type nc = cols (); 808 octave_idx_type nc = cols ();
872 if (nc != a.cols ()) 809 if (nc != a.cols ())
873 { 810 (*current_liboctave_error_handler) ("column dimension mismatch for stack");
874 (*current_liboctave_error_handler)
875 ("column dimension mismatch for stack");
876 return *this;
877 }
878 811
879 octave_idx_type nr_insert = nr; 812 octave_idx_type nr_insert = nr;
880 ComplexMatrix retval (nr + a.rows (), nc); 813 ComplexMatrix retval (nr + a.rows (), nc);
881 retval.insert (*this, 0, 0); 814 retval.insert (*this, 0, 0);
882 retval.insert (a, nr_insert, 0); 815 retval.insert (a, nr_insert, 0);
887 ComplexMatrix::stack (const ComplexMatrix& a) const 820 ComplexMatrix::stack (const ComplexMatrix& a) const
888 { 821 {
889 octave_idx_type nr = rows (); 822 octave_idx_type nr = rows ();
890 octave_idx_type nc = cols (); 823 octave_idx_type nc = cols ();
891 if (nc != a.cols ()) 824 if (nc != a.cols ())
892 { 825 (*current_liboctave_error_handler) ("column dimension mismatch for stack");
893 (*current_liboctave_error_handler)
894 ("column dimension mismatch for stack");
895 return *this;
896 }
897 826
898 octave_idx_type nr_insert = nr; 827 octave_idx_type nr_insert = nr;
899 ComplexMatrix retval (nr + a.rows (), nc); 828 ComplexMatrix retval (nr + a.rows (), nc);
900 retval.insert (*this, 0, 0); 829 retval.insert (*this, 0, 0);
901 retval.insert (a, nr_insert, 0); 830 retval.insert (a, nr_insert, 0);
906 ComplexMatrix::stack (const ComplexRowVector& a) const 835 ComplexMatrix::stack (const ComplexRowVector& a) const
907 { 836 {
908 octave_idx_type nr = rows (); 837 octave_idx_type nr = rows ();
909 octave_idx_type nc = cols (); 838 octave_idx_type nc = cols ();
910 if (nc != a.numel ()) 839 if (nc != a.numel ())
911 { 840 (*current_liboctave_error_handler) ("column dimension mismatch for stack");
912 (*current_liboctave_error_handler)
913 ("column dimension mismatch for stack");
914 return *this;
915 }
916 841
917 octave_idx_type nr_insert = nr; 842 octave_idx_type nr_insert = nr;
918 ComplexMatrix retval (nr + 1, nc); 843 ComplexMatrix retval (nr + 1, nc);
919 retval.insert (*this, 0, 0); 844 retval.insert (*this, 0, 0);
920 retval.insert (a, nr_insert, 0); 845 retval.insert (a, nr_insert, 0);
925 ComplexMatrix::stack (const ComplexColumnVector& a) const 850 ComplexMatrix::stack (const ComplexColumnVector& a) const
926 { 851 {
927 octave_idx_type nr = rows (); 852 octave_idx_type nr = rows ();
928 octave_idx_type nc = cols (); 853 octave_idx_type nc = cols ();
929 if (nc != 1) 854 if (nc != 1)
930 { 855 (*current_liboctave_error_handler) ("column dimension mismatch for stack");
931 (*current_liboctave_error_handler)
932 ("column dimension mismatch for stack");
933 return *this;
934 }
935 856
936 octave_idx_type nr_insert = nr; 857 octave_idx_type nr_insert = nr;
937 ComplexMatrix retval (nr + a.numel (), nc); 858 ComplexMatrix retval (nr + a.numel (), nc);
938 retval.insert (*this, 0, 0); 859 retval.insert (*this, 0, 0);
939 retval.insert (a, nr_insert, 0); 860 retval.insert (a, nr_insert, 0);
944 ComplexMatrix::stack (const ComplexDiagMatrix& a) const 865 ComplexMatrix::stack (const ComplexDiagMatrix& a) const
945 { 866 {
946 octave_idx_type nr = rows (); 867 octave_idx_type nr = rows ();
947 octave_idx_type nc = cols (); 868 octave_idx_type nc = cols ();
948 if (nc != a.cols ()) 869 if (nc != a.cols ())
949 { 870 (*current_liboctave_error_handler) ("column dimension mismatch for stack");
950 (*current_liboctave_error_handler)
951 ("column dimension mismatch for stack");
952 return *this;
953 }
954 871
955 octave_idx_type nr_insert = nr; 872 octave_idx_type nr_insert = nr;
956 ComplexMatrix retval (nr + a.rows (), nc); 873 ComplexMatrix retval (nr + a.rows (), nc);
957 retval.insert (*this, 0, 0); 874 retval.insert (*this, 0, 0);
958 retval.insert (a, nr_insert, 0); 875 retval.insert (a, nr_insert, 0);
1047 octave_idx_type nr = rows (); 964 octave_idx_type nr = rows ();
1048 octave_idx_type nc = cols (); 965 octave_idx_type nc = cols ();
1049 966
1050 if (nr != nc || nr == 0 || nc == 0) 967 if (nr != nc || nr == 0 || nc == 0)
1051 (*current_liboctave_error_handler) ("inverse requires square matrix"); 968 (*current_liboctave_error_handler) ("inverse requires square matrix");
1052 else 969
1053 { 970 int typ = mattype.type ();
1054 int typ = mattype.type (); 971 char uplo = (typ == MatrixType::Lower ? 'L' : 'U');
1055 char uplo = (typ == MatrixType::Lower ? 'L' : 'U'); 972 char udiag = 'N';
1056 char udiag = 'N'; 973 retval = *this;
1057 retval = *this; 974 Complex *tmp_data = retval.fortran_vec ();
1058 Complex *tmp_data = retval.fortran_vec (); 975
1059 976 F77_XFCN (ztrtri, ZTRTRI, (F77_CONST_CHAR_ARG2 (&uplo, 1),
1060 F77_XFCN (ztrtri, ZTRTRI, (F77_CONST_CHAR_ARG2 (&uplo, 1), 977 F77_CONST_CHAR_ARG2 (&udiag, 1),
978 nr, tmp_data, nr, info
979 F77_CHAR_ARG_LEN (1)
980 F77_CHAR_ARG_LEN (1)));
981
982 // Throw-away extra info LAPACK gives so as to not change output.
983 rcon = 0.0;
984 if (info != 0)
985 info = -1;
986 else if (calc_cond)
987 {
988 octave_idx_type ztrcon_info = 0;
989 char job = '1';
990
991 OCTAVE_LOCAL_BUFFER (Complex, cwork, 2*nr);
992 OCTAVE_LOCAL_BUFFER (double, rwork, nr);
993
994 F77_XFCN (ztrcon, ZTRCON, (F77_CONST_CHAR_ARG2 (&job, 1),
995 F77_CONST_CHAR_ARG2 (&uplo, 1),
1061 F77_CONST_CHAR_ARG2 (&udiag, 1), 996 F77_CONST_CHAR_ARG2 (&udiag, 1),
1062 nr, tmp_data, nr, info 997 nr, tmp_data, nr, rcon,
998 cwork, rwork, ztrcon_info
999 F77_CHAR_ARG_LEN (1)
1063 F77_CHAR_ARG_LEN (1) 1000 F77_CHAR_ARG_LEN (1)
1064 F77_CHAR_ARG_LEN (1))); 1001 F77_CHAR_ARG_LEN (1)));
1065 1002
1066 // Throw-away extra info LAPACK gives so as to not change output. 1003 if (ztrcon_info != 0)
1067 rcon = 0.0;
1068 if (info != 0)
1069 info = -1; 1004 info = -1;
1070 else if (calc_cond) 1005 }
1071 { 1006
1072 octave_idx_type ztrcon_info = 0; 1007 if (info == -1 && ! force)
1073 char job = '1'; 1008 retval = *this; // Restore matrix contents.
1074
1075 OCTAVE_LOCAL_BUFFER (Complex, cwork, 2*nr);
1076 OCTAVE_LOCAL_BUFFER (double, rwork, nr);
1077
1078 F77_XFCN (ztrcon, ZTRCON, (F77_CONST_CHAR_ARG2 (&job, 1),
1079 F77_CONST_CHAR_ARG2 (&uplo, 1),
1080 F77_CONST_CHAR_ARG2 (&udiag, 1),
1081 nr, tmp_data, nr, rcon,
1082 cwork, rwork, ztrcon_info
1083 F77_CHAR_ARG_LEN (1)
1084 F77_CHAR_ARG_LEN (1)
1085 F77_CHAR_ARG_LEN (1)));
1086
1087 if (ztrcon_info != 0)
1088 info = -1;
1089 }
1090
1091 if (info == -1 && ! force)
1092 retval = *this; // Restore matrix contents.
1093 }
1094 1009
1095 return retval; 1010 return retval;
1096 } 1011 }
1097 1012
1098 ComplexMatrix 1013 ComplexMatrix
1104 octave_idx_type nr = rows (); 1019 octave_idx_type nr = rows ();
1105 octave_idx_type nc = cols (); 1020 octave_idx_type nc = cols ();
1106 1021
1107 if (nr != nc) 1022 if (nr != nc)
1108 (*current_liboctave_error_handler) ("inverse requires square matrix"); 1023 (*current_liboctave_error_handler) ("inverse requires square matrix");
1024
1025 Array<octave_idx_type> ipvt (dim_vector (nr, 1));
1026 octave_idx_type *pipvt = ipvt.fortran_vec ();
1027
1028 retval = *this;
1029 Complex *tmp_data = retval.fortran_vec ();
1030
1031 Array<Complex> z (dim_vector (1, 1));
1032 octave_idx_type lwork = -1;
1033
1034 // Query the optimum work array size.
1035
1036 F77_XFCN (zgetri, ZGETRI, (nc, tmp_data, nr, pipvt,
1037 z.fortran_vec (), lwork, info));
1038
1039 lwork = static_cast<octave_idx_type> (std::real (z(0)));
1040 lwork = (lwork < 2 *nc ? 2*nc : lwork);
1041 z.resize (dim_vector (lwork, 1));
1042 Complex *pz = z.fortran_vec ();
1043
1044 info = 0;
1045
1046 // Calculate the norm of the matrix, for later use.
1047 double anorm;
1048 //if (calc_cond) // Must always calculate anorm for bug #45577
1049 anorm = retval.abs ().sum ().row (static_cast<octave_idx_type>(0)).max ();
1050
1051 // Work around bug #45577, LAPACK crashes Octave if norm is NaN
1052 if (xisnan (anorm))
1053 info = -1;
1109 else 1054 else
1110 { 1055 F77_XFCN (zgetrf, ZGETRF, (nc, nc, tmp_data, nr, pipvt, info));
1111 Array<octave_idx_type> ipvt (dim_vector (nr, 1)); 1056
1112 octave_idx_type *pipvt = ipvt.fortran_vec (); 1057 // Throw-away extra info LAPACK gives so as to not change output.
1113 1058 rcon = 0.0;
1114 retval = *this; 1059 if (info != 0)
1115 Complex *tmp_data = retval.fortran_vec (); 1060 info = -1;
1116 1061 else if (calc_cond)
1117 Array<Complex> z (dim_vector (1, 1)); 1062 {
1118 octave_idx_type lwork = -1; 1063 // Now calculate the condition number for non-singular matrix.
1119 1064 octave_idx_type zgecon_info = 0;
1120 // Query the optimum work array size. 1065 char job = '1';
1066 Array<double> rz (dim_vector (2 * nc, 1));
1067 double *prz = rz.fortran_vec ();
1068 F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
1069 nc, tmp_data, nr, anorm,
1070 rcon, pz, prz, zgecon_info
1071 F77_CHAR_ARG_LEN (1)));
1072
1073 if (zgecon_info != 0)
1074 info = -1;
1075 }
1076
1077 if (info == -1 && ! force)
1078 retval = *this; // Restore contents.
1079 else
1080 {
1081 octave_idx_type zgetri_info = 0;
1121 1082
1122 F77_XFCN (zgetri, ZGETRI, (nc, tmp_data, nr, pipvt, 1083 F77_XFCN (zgetri, ZGETRI, (nc, tmp_data, nr, pipvt,
1123 z.fortran_vec (), lwork, info)); 1084 pz, lwork, zgetri_info));
1124 1085
1125 lwork = static_cast<octave_idx_type> (std::real (z(0))); 1086 if (zgetri_info != 0)
1126 lwork = (lwork < 2 *nc ? 2*nc : lwork);
1127 z.resize (dim_vector (lwork, 1));
1128 Complex *pz = z.fortran_vec ();
1129
1130 info = 0;
1131
1132 // Calculate the norm of the matrix, for later use.
1133 double anorm;
1134 //if (calc_cond) // Must always calculate anorm for bug #45577
1135 anorm = retval.abs ().sum ().row (static_cast<octave_idx_type>(0)).max ();
1136
1137 // Work around bug #45577, LAPACK crashes Octave if norm is NaN
1138 if (xisnan (anorm))
1139 info = -1; 1087 info = -1;
1140 else 1088 }
1141 F77_XFCN (zgetrf, ZGETRF, (nc, nc, tmp_data, nr, pipvt, info)); 1089
1142 1090 if (info != 0)
1143 // Throw-away extra info LAPACK gives so as to not change output. 1091 mattype.mark_as_rectangular ();
1144 rcon = 0.0;
1145 if (info != 0)
1146 info = -1;
1147 else if (calc_cond)
1148 {
1149 // Now calculate the condition number for non-singular matrix.
1150 octave_idx_type zgecon_info = 0;
1151 char job = '1';
1152 Array<double> rz (dim_vector (2 * nc, 1));
1153 double *prz = rz.fortran_vec ();
1154 F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
1155 nc, tmp_data, nr, anorm,
1156 rcon, pz, prz, zgecon_info
1157 F77_CHAR_ARG_LEN (1)));
1158
1159 if (zgecon_info != 0)
1160 info = -1;
1161 }
1162
1163 if (info == -1 && ! force)
1164 retval = *this; // Restore contents.
1165 else
1166 {
1167 octave_idx_type zgetri_info = 0;
1168
1169 F77_XFCN (zgetri, ZGETRI, (nc, tmp_data, nr, pipvt,
1170 pz, lwork, zgetri_info));
1171
1172 if (zgetri_info != 0)
1173 info = -1;
1174 }
1175
1176 if (info != 0)
1177 mattype.mark_as_rectangular ();
1178 }
1179 1092
1180 return retval; 1093 return retval;
1181 } 1094 }
1182 1095
1183 ComplexMatrix 1096 ComplexMatrix
1625 octave_idx_type nr = rows (); 1538 octave_idx_type nr = rows ();
1626 octave_idx_type nc = cols (); 1539 octave_idx_type nc = cols ();
1627 1540
1628 if (nr != nc) 1541 if (nr != nc)
1629 (*current_liboctave_error_handler) ("matrix must be square"); 1542 (*current_liboctave_error_handler) ("matrix must be square");
1630 else 1543
1631 { 1544 volatile int typ = mattype.type ();
1632 volatile int typ = mattype.type (); 1545
1633 1546 // Even though the matrix is marked as singular (Rectangular), we may
1634 // Even though the matrix is marked as singular (Rectangular), we may 1547 // still get a useful number from the LU factorization, because it always
1635 // still get a useful number from the LU factorization, because it always 1548 // completes.
1636 // completes. 1549
1637 1550 if (typ == MatrixType::Unknown)
1638 if (typ == MatrixType::Unknown) 1551 typ = mattype.type (*this);
1639 typ = mattype.type (*this); 1552 else if (typ == MatrixType::Rectangular)
1640 else if (typ == MatrixType::Rectangular) 1553 typ = MatrixType::Full;
1641 typ = MatrixType::Full; 1554
1642 1555 if (typ == MatrixType::Lower || typ == MatrixType::Upper)
1643 if (typ == MatrixType::Lower || typ == MatrixType::Upper) 1556 {
1557 for (octave_idx_type i = 0; i < nc; i++)
1558 retval *= elem (i,i);
1559 }
1560 else if (typ == MatrixType::Hermitian)
1561 {
1562 ComplexMatrix atmp = *this;
1563 Complex *tmp_data = atmp.fortran_vec ();
1564
1565 double anorm = 0;
1566 if (calc_cond) anorm = xnorm (*this, 1);
1567
1568
1569 char job = 'L';
1570 F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr,
1571 tmp_data, nr, info
1572 F77_CHAR_ARG_LEN (1)));
1573
1574 if (info != 0)
1644 { 1575 {
1576 rcon = 0.0;
1577 mattype.mark_as_unsymmetric ();
1578 typ = MatrixType::Full;
1579 }
1580 else
1581 {
1582 Array<Complex> z (dim_vector (2 * nc, 1));
1583 Complex *pz = z.fortran_vec ();
1584 Array<double> rz (dim_vector (nc, 1));
1585 double *prz = rz.fortran_vec ();
1586
1587 F77_XFCN (zpocon, ZPOCON, (F77_CONST_CHAR_ARG2 (&job, 1),
1588 nr, tmp_data, nr, anorm,
1589 rcon, pz, prz, info
1590 F77_CHAR_ARG_LEN (1)));
1591
1592 if (info != 0)
1593 rcon = 0.0;
1594
1645 for (octave_idx_type i = 0; i < nc; i++) 1595 for (octave_idx_type i = 0; i < nc; i++)
1646 retval *= elem (i,i); 1596 retval *= atmp (i,i);
1597
1598 retval = retval.square ();
1647 } 1599 }
1648 else if (typ == MatrixType::Hermitian) 1600 }
1601 else if (typ != MatrixType::Full)
1602 (*current_liboctave_error_handler) ("det: invalid dense matrix type");
1603
1604 if (typ == MatrixType::Full)
1605 {
1606 Array<octave_idx_type> ipvt (dim_vector (nr, 1));
1607 octave_idx_type *pipvt = ipvt.fortran_vec ();
1608
1609 ComplexMatrix atmp = *this;
1610 Complex *tmp_data = atmp.fortran_vec ();
1611
1612 info = 0;
1613
1614 // Calculate the norm of the matrix, for later use.
1615 double anorm = 0;
1616 //if (calc_cond) // Must always calculate anorm for bug #45577
1617 anorm = xnorm (*this, 1);
1618
1619 // Work around bug #45577, LAPACK crashes Octave if norm is NaN
1620 if (xisnan (anorm))
1621 info = -1;
1622 else
1623 F77_XFCN (zgetrf, ZGETRF, (nr, nr, tmp_data, nr, pipvt, info));
1624
1625 // Throw-away extra info LAPACK gives so as to not change output.
1626 rcon = 0.0;
1627 if (info != 0)
1649 { 1628 {
1650 ComplexMatrix atmp = *this; 1629 info = -1;
1651 Complex *tmp_data = atmp.fortran_vec (); 1630 retval = ComplexDET ();
1652 1631 }
1653 double anorm = 0; 1632 else
1654 if (calc_cond) anorm = xnorm (*this, 1); 1633 {
1655 1634 if (calc_cond)
1656
1657 char job = 'L';
1658 F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr,
1659 tmp_data, nr, info
1660 F77_CHAR_ARG_LEN (1)));
1661
1662 if (info != 0)
1663 { 1635 {
1664 rcon = 0.0; 1636 // Now calc the condition number for non-singular matrix.
1665 mattype.mark_as_unsymmetric (); 1637 char job = '1';
1666 typ = MatrixType::Full;
1667 }
1668 else
1669 {
1670 Array<Complex> z (dim_vector (2 * nc, 1)); 1638 Array<Complex> z (dim_vector (2 * nc, 1));
1671 Complex *pz = z.fortran_vec (); 1639 Complex *pz = z.fortran_vec ();
1672 Array<double> rz (dim_vector (nc, 1)); 1640 Array<double> rz (dim_vector (2 * nc, 1));
1673 double *prz = rz.fortran_vec (); 1641 double *prz = rz.fortran_vec ();
1674 1642
1675 F77_XFCN (zpocon, ZPOCON, (F77_CONST_CHAR_ARG2 (&job, 1), 1643 F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
1676 nr, tmp_data, nr, anorm, 1644 nc, tmp_data, nr, anorm,
1677 rcon, pz, prz, info 1645 rcon, pz, prz, info
1678 F77_CHAR_ARG_LEN (1))); 1646 F77_CHAR_ARG_LEN (1)));
1679
1680 if (info != 0)
1681 rcon = 0.0;
1682
1683 for (octave_idx_type i = 0; i < nc; i++)
1684 retval *= atmp (i,i);
1685
1686 retval = retval.square ();
1687 } 1647 }
1688 } 1648
1689 else if (typ != MatrixType::Full)
1690 (*current_liboctave_error_handler) ("det: invalid dense matrix type");
1691
1692 if (typ == MatrixType::Full)
1693 {
1694 Array<octave_idx_type> ipvt (dim_vector (nr, 1));
1695 octave_idx_type *pipvt = ipvt.fortran_vec ();
1696
1697 ComplexMatrix atmp = *this;
1698 Complex *tmp_data = atmp.fortran_vec ();
1699
1700 info = 0;
1701
1702 // Calculate the norm of the matrix, for later use.
1703 double anorm = 0;
1704 //if (calc_cond) // Must always calculate anorm for bug #45577
1705 anorm = xnorm (*this, 1);
1706
1707 // Work around bug #45577, LAPACK crashes Octave if norm is NaN
1708 if (xisnan (anorm))
1709 info = -1;
1710 else
1711 F77_XFCN (zgetrf, ZGETRF, (nr, nr, tmp_data, nr, pipvt, info));
1712
1713 // Throw-away extra info LAPACK gives so as to not change output.
1714 rcon = 0.0;
1715 if (info != 0) 1649 if (info != 0)
1716 { 1650 {
1717 info = -1; 1651 info = -1;
1718 retval = ComplexDET (); 1652 retval = ComplexDET ();
1719 } 1653 }
1720 else 1654 else
1721 { 1655 {
1722 if (calc_cond) 1656 for (octave_idx_type i = 0; i < nc; i++)
1723 { 1657 {
1724 // Now calc the condition number for non-singular matrix. 1658 Complex c = atmp(i,i);
1725 char job = '1'; 1659 retval *= (ipvt(i) != (i+1)) ? -c : c;
1726 Array<Complex> z (dim_vector (2 * nc, 1));
1727 Complex *pz = z.fortran_vec ();
1728 Array<double> rz (dim_vector (2 * nc, 1));
1729 double *prz = rz.fortran_vec ();
1730
1731 F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
1732 nc, tmp_data, nr, anorm,
1733 rcon, pz, prz, info
1734 F77_CHAR_ARG_LEN (1)));
1735 }
1736
1737 if (info != 0)
1738 {
1739 info = -1;
1740 retval = ComplexDET ();
1741 }
1742 else
1743 {
1744 for (octave_idx_type i = 0; i < nc; i++)
1745 {
1746 Complex c = atmp(i,i);
1747 retval *= (ipvt(i) != (i+1)) ? -c : c;
1748 }
1749 } 1660 }
1750 } 1661 }
1751 } 1662 }
1752 } 1663 }
1753 1664
1768 octave_idx_type nr = rows (); 1679 octave_idx_type nr = rows ();
1769 octave_idx_type nc = cols (); 1680 octave_idx_type nc = cols ();
1770 1681
1771 if (nr != nc) 1682 if (nr != nc)
1772 (*current_liboctave_error_handler) ("matrix must be square"); 1683 (*current_liboctave_error_handler) ("matrix must be square");
1773 else if (nr == 0 || nc == 0) 1684
1685 if (nr == 0 || nc == 0)
1774 rcon = octave_Inf; 1686 rcon = octave_Inf;
1775 else 1687 else
1776 { 1688 {
1777 volatile int typ = mattype.type (); 1689 volatile int typ = mattype.type ();
1778 1690
1942 octave_idx_type nc = cols (); 1854 octave_idx_type nc = cols ();
1943 1855
1944 if (nr != b.rows ()) 1856 if (nr != b.rows ())
1945 (*current_liboctave_error_handler) 1857 (*current_liboctave_error_handler)
1946 ("matrix dimension mismatch solution of linear equations"); 1858 ("matrix dimension mismatch solution of linear equations");
1947 else if (nr == 0 || nc == 0 || b.cols () == 0) 1859
1860 if (nr == 0 || nc == 0 || b.cols () == 0)
1948 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0)); 1861 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
1949 else 1862 else
1950 { 1863 {
1951 volatile int typ = mattype.type (); 1864 volatile int typ = mattype.type ();
1952 1865
1953 if (typ == MatrixType::Permuted_Upper || typ == MatrixType::Upper) 1866 if (typ != MatrixType::Permuted_Upper && typ != MatrixType::Upper)
1867 (*current_liboctave_error_handler) ("incorrect matrix type");
1868
1869 octave_idx_type b_nc = b.cols ();
1870 rcon = 1.;
1871 info = 0;
1872
1873 if (typ == MatrixType::Permuted_Upper)
1874 (*current_liboctave_error_handler)
1875 ("permuted triangular matrix not implemented");
1876
1877 const Complex *tmp_data = fortran_vec ();
1878
1879 retval = b;
1880 Complex *result = retval.fortran_vec ();
1881
1882 char uplo = 'U';
1883 char trans = get_blas_char (transt);
1884 char dia = 'N';
1885
1886 F77_XFCN (ztrtrs, ZTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1),
1887 F77_CONST_CHAR_ARG2 (&trans, 1),
1888 F77_CONST_CHAR_ARG2 (&dia, 1),
1889 nr, b_nc, tmp_data, nr,
1890 result, nr, info
1891 F77_CHAR_ARG_LEN (1)
1892 F77_CHAR_ARG_LEN (1)
1893 F77_CHAR_ARG_LEN (1)));
1894
1895 if (calc_cond)
1954 { 1896 {
1955 octave_idx_type b_nc = b.cols (); 1897 char norm = '1';
1956 rcon = 1.; 1898 uplo = 'U';
1957 info = 0; 1899 dia = 'N';
1958 1900
1959 if (typ == MatrixType::Permuted_Upper) 1901 Array<Complex> z (dim_vector (2 * nc, 1));
1902 Complex *pz = z.fortran_vec ();
1903 Array<double> rz (dim_vector (nc, 1));
1904 double *prz = rz.fortran_vec ();
1905
1906 F77_XFCN (ztrcon, ZTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),
1907 F77_CONST_CHAR_ARG2 (&uplo, 1),
1908 F77_CONST_CHAR_ARG2 (&dia, 1),
1909 nr, tmp_data, nr, rcon,
1910 pz, prz, info
1911 F77_CHAR_ARG_LEN (1)
1912 F77_CHAR_ARG_LEN (1)
1913 F77_CHAR_ARG_LEN (1)));
1914
1915 if (info != 0)
1916 info = -2;
1917
1918 volatile double rcond_plus_one = rcon + 1.0;
1919
1920 if (rcond_plus_one == 1.0 || xisnan (rcon))
1960 { 1921 {
1961 (*current_liboctave_error_handler) 1922 info = -2;
1962 ("permuted triangular matrix not implemented"); 1923
1963 } 1924 if (sing_handler)
1964 else 1925 sing_handler (rcon);
1965 { 1926 else
1966 const Complex *tmp_data = fortran_vec (); 1927 warn_singular_matrix (rcon);
1967
1968 retval = b;
1969 Complex *result = retval.fortran_vec ();
1970
1971 char uplo = 'U';
1972 char trans = get_blas_char (transt);
1973 char dia = 'N';
1974
1975 F77_XFCN (ztrtrs, ZTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1),
1976 F77_CONST_CHAR_ARG2 (&trans, 1),
1977 F77_CONST_CHAR_ARG2 (&dia, 1),
1978 nr, b_nc, tmp_data, nr,
1979 result, nr, info
1980 F77_CHAR_ARG_LEN (1)
1981 F77_CHAR_ARG_LEN (1)
1982 F77_CHAR_ARG_LEN (1)));
1983
1984 if (calc_cond)
1985 {
1986 char norm = '1';
1987 uplo = 'U';
1988 dia = 'N';
1989
1990 Array<Complex> z (dim_vector (2 * nc, 1));
1991 Complex *pz = z.fortran_vec ();
1992 Array<double> rz (dim_vector (nc, 1));
1993 double *prz = rz.fortran_vec ();
1994
1995 F77_XFCN (ztrcon, ZTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),
1996 F77_CONST_CHAR_ARG2 (&uplo, 1),
1997 F77_CONST_CHAR_ARG2 (&dia, 1),
1998 nr, tmp_data, nr, rcon,
1999 pz, prz, info
2000 F77_CHAR_ARG_LEN (1)
2001 F77_CHAR_ARG_LEN (1)
2002 F77_CHAR_ARG_LEN (1)));
2003
2004 if (info != 0)
2005 info = -2;
2006
2007 volatile double rcond_plus_one = rcon + 1.0;
2008
2009 if (rcond_plus_one == 1.0 || xisnan (rcon))
2010 {
2011 info = -2;
2012
2013 if (sing_handler)
2014 sing_handler (rcon);
2015 else
2016 warn_singular_matrix (rcon);
2017 }
2018 }
2019 } 1928 }
2020 } 1929 }
2021 else
2022 (*current_liboctave_error_handler) ("incorrect matrix type");
2023 } 1930 }
2024 1931
2025 return retval; 1932 return retval;
2026 } 1933 }
2027 1934
2037 octave_idx_type nc = cols (); 1944 octave_idx_type nc = cols ();
2038 1945
2039 if (nr != b.rows ()) 1946 if (nr != b.rows ())
2040 (*current_liboctave_error_handler) 1947 (*current_liboctave_error_handler)
2041 ("matrix dimension mismatch solution of linear equations"); 1948 ("matrix dimension mismatch solution of linear equations");
2042 else if (nr == 0 || nc == 0 || b.cols () == 0) 1949
1950 if (nr == 0 || nc == 0 || b.cols () == 0)
2043 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0)); 1951 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
2044 else 1952 else
2045 { 1953 {
2046 volatile int typ = mattype.type (); 1954 volatile int typ = mattype.type ();
2047 1955
2048 if (typ == MatrixType::Permuted_Lower || typ == MatrixType::Lower) 1956 if (typ != MatrixType::Permuted_Lower && typ != MatrixType::Lower)
1957 (*current_liboctave_error_handler) ("incorrect matrix type");
1958
1959 octave_idx_type b_nc = b.cols ();
1960 rcon = 1.;
1961 info = 0;
1962
1963 if (typ == MatrixType::Permuted_Lower)
1964 (*current_liboctave_error_handler)
1965 ("permuted triangular matrix not implemented");
1966
1967 const Complex *tmp_data = fortran_vec ();
1968
1969 retval = b;
1970 Complex *result = retval.fortran_vec ();
1971
1972 char uplo = 'L';
1973 char trans = get_blas_char (transt);
1974 char dia = 'N';
1975
1976 F77_XFCN (ztrtrs, ZTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1),
1977 F77_CONST_CHAR_ARG2 (&trans, 1),
1978 F77_CONST_CHAR_ARG2 (&dia, 1),
1979 nr, b_nc, tmp_data, nr,
1980 result, nr, info
1981 F77_CHAR_ARG_LEN (1)
1982 F77_CHAR_ARG_LEN (1)
1983 F77_CHAR_ARG_LEN (1)));
1984
1985 if (calc_cond)
2049 { 1986 {
2050 octave_idx_type b_nc = b.cols (); 1987 char norm = '1';
2051 rcon = 1.; 1988 uplo = 'L';
2052 info = 0; 1989 dia = 'N';
2053 1990
2054 if (typ == MatrixType::Permuted_Lower) 1991 Array<Complex> z (dim_vector (2 * nc, 1));
1992 Complex *pz = z.fortran_vec ();
1993 Array<double> rz (dim_vector (nc, 1));
1994 double *prz = rz.fortran_vec ();
1995
1996 F77_XFCN (ztrcon, ZTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),
1997 F77_CONST_CHAR_ARG2 (&uplo, 1),
1998 F77_CONST_CHAR_ARG2 (&dia, 1),
1999 nr, tmp_data, nr, rcon,
2000 pz, prz, info
2001 F77_CHAR_ARG_LEN (1)
2002 F77_CHAR_ARG_LEN (1)
2003 F77_CHAR_ARG_LEN (1)));
2004
2005 if (info != 0)
2006 info = -2;
2007
2008 volatile double rcond_plus_one = rcon + 1.0;
2009
2010 if (rcond_plus_one == 1.0 || xisnan (rcon))
2055 { 2011 {
2056 (*current_liboctave_error_handler) 2012 info = -2;
2057 ("permuted triangular matrix not implemented"); 2013
2058 } 2014 if (sing_handler)
2059 else 2015 sing_handler (rcon);
2060 { 2016 else
2061 const Complex *tmp_data = fortran_vec (); 2017 warn_singular_matrix (rcon);
2062
2063 retval = b;
2064 Complex *result = retval.fortran_vec ();
2065
2066 char uplo = 'L';
2067 char trans = get_blas_char (transt);
2068 char dia = 'N';
2069
2070 F77_XFCN (ztrtrs, ZTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1),
2071 F77_CONST_CHAR_ARG2 (&trans, 1),
2072 F77_CONST_CHAR_ARG2 (&dia, 1),
2073 nr, b_nc, tmp_data, nr,
2074 result, nr, info
2075 F77_CHAR_ARG_LEN (1)
2076 F77_CHAR_ARG_LEN (1)
2077 F77_CHAR_ARG_LEN (1)));
2078
2079 if (calc_cond)
2080 {
2081 char norm = '1';
2082 uplo = 'L';
2083 dia = 'N';
2084
2085 Array<Complex> z (dim_vector (2 * nc, 1));
2086 Complex *pz = z.fortran_vec ();
2087 Array<double> rz (dim_vector (nc, 1));
2088 double *prz = rz.fortran_vec ();
2089
2090 F77_XFCN (ztrcon, ZTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),
2091 F77_CONST_CHAR_ARG2 (&uplo, 1),
2092 F77_CONST_CHAR_ARG2 (&dia, 1),
2093 nr, tmp_data, nr, rcon,
2094 pz, prz, info
2095 F77_CHAR_ARG_LEN (1)
2096 F77_CHAR_ARG_LEN (1)
2097 F77_CHAR_ARG_LEN (1)));
2098
2099 if (info != 0)
2100 info = -2;
2101
2102 volatile double rcond_plus_one = rcon + 1.0;
2103
2104 if (rcond_plus_one == 1.0 || xisnan (rcon))
2105 {
2106 info = -2;
2107
2108 if (sing_handler)
2109 sing_handler (rcon);
2110 else
2111 warn_singular_matrix (rcon);
2112 }
2113 }
2114 } 2018 }
2115 } 2019 }
2116 else
2117 (*current_liboctave_error_handler) ("incorrect matrix type");
2118 } 2020 }
2119 2021
2120 return retval; 2022 return retval;
2121 } 2023 }
2122 2024
2133 2035
2134 2036
2135 if (nr != nc || nr != b.rows ()) 2037 if (nr != nc || nr != b.rows ())
2136 (*current_liboctave_error_handler) 2038 (*current_liboctave_error_handler)
2137 ("matrix dimension mismatch solution of linear equations"); 2039 ("matrix dimension mismatch solution of linear equations");
2138 else if (nr == 0 || b.cols () == 0) 2040
2041 if (nr == 0 || b.cols () == 0)
2139 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0)); 2042 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
2140 else 2043 else
2141 { 2044 {
2142 volatile int typ = mattype.type (); 2045 volatile int typ = mattype.type ();
2143 2046
2385 retval = hermitian ().solve (mattype, b, info, rcon, sing_handler, 2288 retval = hermitian ().solve (mattype, b, info, rcon, sing_handler,
2386 singular_fallback); 2289 singular_fallback);
2387 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian) 2290 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
2388 retval = fsolve (mattype, b, info, rcon, sing_handler, true); 2291 retval = fsolve (mattype, b, info, rcon, sing_handler, true);
2389 else if (typ != MatrixType::Rectangular) 2292 else if (typ != MatrixType::Rectangular)
2390 { 2293 (*current_liboctave_error_handler) ("unknown matrix type");
2391 (*current_liboctave_error_handler) ("unknown matrix type");
2392 return ComplexMatrix ();
2393 }
2394 2294
2395 // Rectangular or one of the above solvers flags a singular matrix 2295 // Rectangular or one of the above solvers flags a singular matrix
2396 if (singular_fallback && mattype.type () == MatrixType::Rectangular) 2296 if (singular_fallback && mattype.type () == MatrixType::Rectangular)
2397 { 2297 {
2398 octave_idx_type rank; 2298 octave_idx_type rank;
2664 octave_idx_type n = cols (); 2564 octave_idx_type n = cols ();
2665 2565
2666 if (m != b.rows ()) 2566 if (m != b.rows ())
2667 (*current_liboctave_error_handler) 2567 (*current_liboctave_error_handler)
2668 ("matrix dimension mismatch solution of linear equations"); 2568 ("matrix dimension mismatch solution of linear equations");
2669 else if (m== 0 || n == 0 || b.cols () == 0) 2569
2570 if (m== 0 || n == 0 || b.cols () == 0)
2670 retval = ComplexMatrix (n, b.cols (), Complex (0.0, 0.0)); 2571 retval = ComplexMatrix (n, b.cols (), Complex (0.0, 0.0));
2671 else 2572 else
2672 { 2573 {
2673 volatile octave_idx_type minmn = (m < n ? m : n); 2574 volatile octave_idx_type minmn = (m < n ? m : n);
2674 octave_idx_type maxmn = m > n ? m : n; 2575 octave_idx_type maxmn = m > n ? m : n;
2861 octave_idx_type n = cols (); 2762 octave_idx_type n = cols ();
2862 2763
2863 if (m != b.numel ()) 2764 if (m != b.numel ())
2864 (*current_liboctave_error_handler) 2765 (*current_liboctave_error_handler)
2865 ("matrix dimension mismatch solution of linear equations"); 2766 ("matrix dimension mismatch solution of linear equations");
2866 else if (m == 0 || n == 0 || b.cols () == 0) 2767
2768 if (m == 0 || n == 0 || b.cols () == 0)
2867 retval = ComplexColumnVector (n, Complex (0.0, 0.0)); 2769 retval = ComplexColumnVector (n, Complex (0.0, 0.0));
2868 else 2770 else
2869 { 2771 {
2870 volatile octave_idx_type minmn = (m < n ? m : n); 2772 volatile octave_idx_type minmn = (m < n ? m : n);
2871 octave_idx_type maxmn = m > n ? m : n; 2773 octave_idx_type maxmn = m > n ? m : n;
3168 } 3070 }
3169 3071
3170 ComplexDiagMatrix 3072 ComplexDiagMatrix
3171 ComplexMatrix::diag (octave_idx_type m, octave_idx_type n) const 3073 ComplexMatrix::diag (octave_idx_type m, octave_idx_type n) const
3172 { 3074 {
3173 ComplexDiagMatrix retval; 3075 octave_idx_type nr = rows ();
3174 3076 octave_idx_type nc = cols ();
3175 octave_idx_type nr = rows (); 3077
3176 octave_idx_type nc = cols (); 3078 if (nr != 1 && nc != 1)
3177 3079 (*current_liboctave_error_handler) ("diag: expecting vector argument");
3178 if (nr == 1 || nc == 1) 3080
3179 retval = ComplexDiagMatrix (*this, m, n); 3081 return ComplexDiagMatrix (*this, m, n);
3180 else
3181 (*current_liboctave_error_handler)
3182 ("diag: expecting vector argument");
3183
3184 return retval;
3185 } 3082 }
3186 3083
3187 bool 3084 bool
3188 ComplexMatrix::row_is_real_only (octave_idx_type i) const 3085 ComplexMatrix::row_is_real_only (octave_idx_type i) const
3189 { 3086 {
3856 { 3753 {
3857 octave_idx_type nr = a.rows (); 3754 octave_idx_type nr = a.rows ();
3858 octave_idx_type nc = a.columns (); 3755 octave_idx_type nc = a.columns ();
3859 3756
3860 if (nr != b.rows () || nc != b.columns ()) 3757 if (nr != b.rows () || nc != b.columns ())
3861 { 3758 (*current_liboctave_error_handler)
3862 (*current_liboctave_error_handler) 3759 ("two-arg min requires same size arguments");
3863 ("two-arg min requires same size arguments");
3864 return ComplexMatrix ();
3865 }
3866 3760
3867 EMPTY_RETURN_CHECK (ComplexMatrix); 3761 EMPTY_RETURN_CHECK (ComplexMatrix);
3868 3762
3869 ComplexMatrix result (nr, nc); 3763 ComplexMatrix result (nr, nc);
3870 3764
3944 { 3838 {
3945 octave_idx_type nr = a.rows (); 3839 octave_idx_type nr = a.rows ();
3946 octave_idx_type nc = a.columns (); 3840 octave_idx_type nc = a.columns ();
3947 3841
3948 if (nr != b.rows () || nc != b.columns ()) 3842 if (nr != b.rows () || nc != b.columns ())
3949 { 3843 (*current_liboctave_error_handler)
3950 (*current_liboctave_error_handler) 3844 ("two-arg max requires same size arguments");
3951 ("two-arg max requires same size arguments");
3952 return ComplexMatrix ();
3953 }
3954 3845
3955 EMPTY_RETURN_CHECK (ComplexMatrix); 3846 EMPTY_RETURN_CHECK (ComplexMatrix);
3956 3847
3957 ComplexMatrix result (nr, nc); 3848 ComplexMatrix result (nr, nc);
3958 3849