Mercurial > jwe > octave

comparison liboctave/array/fMatrix.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 40051830f89b
comparison
equal deleted inserted replaced
21135:95da3bc8a281 21136:7cac4e7458f2
347 octave_idx_type r, octave_idx_type c) 347 octave_idx_type r, octave_idx_type c)
348 { 348 {
349 octave_idx_type a_len = a.numel (); 349 octave_idx_type a_len = a.numel ();
350 350
351 if (r < 0 || r >= rows () || c < 0 || c + a_len > cols ()) 351 if (r < 0 || r >= rows () || c < 0 || c + a_len > cols ())
352 { 352 (*current_liboctave_error_handler) ("range error for insert");
353 (*current_liboctave_error_handler) ("range error for insert");
354 return *this;
355 }
356 353
357 if (a_len > 0) 354 if (a_len > 0)
358 { 355 {
359 make_unique (); 356 make_unique ();
360 357
370 octave_idx_type r, octave_idx_type c) 367 octave_idx_type r, octave_idx_type c)
371 { 368 {
372 octave_idx_type a_len = a.numel (); 369 octave_idx_type a_len = a.numel ();
373 370
374 if (r < 0 || r + a_len > rows () || c < 0 || c >= cols ()) 371 if (r < 0 || r + a_len > rows () || c < 0 || c >= cols ())
375 { 372 (*current_liboctave_error_handler) ("range error for insert");
376 (*current_liboctave_error_handler) ("range error for insert");
377 return *this;
378 }
379 373
380 if (a_len > 0) 374 if (a_len > 0)
381 { 375 {
382 make_unique (); 376 make_unique ();
383 377
394 { 388 {
395 octave_idx_type a_nr = a.rows (); 389 octave_idx_type a_nr = a.rows ();
396 octave_idx_type a_nc = a.cols (); 390 octave_idx_type a_nc = a.cols ();
397 391
398 if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ()) 392 if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ())
399 { 393 (*current_liboctave_error_handler) ("range error for insert");
400 (*current_liboctave_error_handler) ("range error for insert");
401 return *this;
402 }
403 394
404 fill (0.0, r, c, r + a_nr - 1, c + a_nc - 1); 395 fill (0.0, r, c, r + a_nr - 1, c + a_nc - 1);
405 396
406 octave_idx_type a_len = a.length (); 397 octave_idx_type a_len = a.length ();
407 398
441 octave_idx_type nr = rows (); 432 octave_idx_type nr = rows ();
442 octave_idx_type nc = cols (); 433 octave_idx_type nc = cols ();
443 434
444 if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0 435 if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0
445 || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc) 436 || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc)
446 { 437 (*current_liboctave_error_handler) ("range error for fill");
447 (*current_liboctave_error_handler) ("range error for fill");
448 return *this;
449 }
450 438
451 if (r1 > r2) { std::swap (r1, r2); } 439 if (r1 > r2) { std::swap (r1, r2); }
452 if (c1 > c2) { std::swap (c1, c2); } 440 if (c1 > c2) { std::swap (c1, c2); }
453 441
454 if (r2 >= r1 && c2 >= c1) 442 if (r2 >= r1 && c2 >= c1)
467 FloatMatrix::append (const FloatMatrix& a) const 455 FloatMatrix::append (const FloatMatrix& a) const
468 { 456 {
469 octave_idx_type nr = rows (); 457 octave_idx_type nr = rows ();
470 octave_idx_type nc = cols (); 458 octave_idx_type nc = cols ();
471 if (nr != a.rows ()) 459 if (nr != a.rows ())
472 { 460 (*current_liboctave_error_handler) ("row dimension mismatch for append");
473 (*current_liboctave_error_handler) ("row dimension mismatch for append");
474 return FloatMatrix ();
475 }
476 461
477 octave_idx_type nc_insert = nc; 462 octave_idx_type nc_insert = nc;
478 FloatMatrix retval (nr, nc + a.cols ()); 463 FloatMatrix retval (nr, nc + a.cols ());
479 retval.insert (*this, 0, 0); 464 retval.insert (*this, 0, 0);
480 retval.insert (a, 0, nc_insert); 465 retval.insert (a, 0, nc_insert);
485 FloatMatrix::append (const FloatRowVector& a) const 470 FloatMatrix::append (const FloatRowVector& a) const
486 { 471 {
487 octave_idx_type nr = rows (); 472 octave_idx_type nr = rows ();
488 octave_idx_type nc = cols (); 473 octave_idx_type nc = cols ();
489 if (nr != 1) 474 if (nr != 1)
490 { 475 (*current_liboctave_error_handler) ("row dimension mismatch for append");
491 (*current_liboctave_error_handler) ("row dimension mismatch for append");
492 return FloatMatrix ();
493 }
494 476
495 octave_idx_type nc_insert = nc; 477 octave_idx_type nc_insert = nc;
496 FloatMatrix retval (nr, nc + a.numel ()); 478 FloatMatrix retval (nr, nc + a.numel ());
497 retval.insert (*this, 0, 0); 479 retval.insert (*this, 0, 0);
498 retval.insert (a, 0, nc_insert); 480 retval.insert (a, 0, nc_insert);
503 FloatMatrix::append (const FloatColumnVector& a) const 485 FloatMatrix::append (const FloatColumnVector& a) const
504 { 486 {
505 octave_idx_type nr = rows (); 487 octave_idx_type nr = rows ();
506 octave_idx_type nc = cols (); 488 octave_idx_type nc = cols ();
507 if (nr != a.numel ()) 489 if (nr != a.numel ())
508 { 490 (*current_liboctave_error_handler) ("row dimension mismatch for append");
509 (*current_liboctave_error_handler) ("row dimension mismatch for append");
510 return FloatMatrix ();
511 }
512 491
513 octave_idx_type nc_insert = nc; 492 octave_idx_type nc_insert = nc;
514 FloatMatrix retval (nr, nc + 1); 493 FloatMatrix retval (nr, nc + 1);
515 retval.insert (*this, 0, 0); 494 retval.insert (*this, 0, 0);
516 retval.insert (a, 0, nc_insert); 495 retval.insert (a, 0, nc_insert);
521 FloatMatrix::append (const FloatDiagMatrix& a) const 500 FloatMatrix::append (const FloatDiagMatrix& a) const
522 { 501 {
523 octave_idx_type nr = rows (); 502 octave_idx_type nr = rows ();
524 octave_idx_type nc = cols (); 503 octave_idx_type nc = cols ();
525 if (nr != a.rows ()) 504 if (nr != a.rows ())
526 { 505 (*current_liboctave_error_handler) ("row dimension mismatch for append");
527 (*current_liboctave_error_handler) ("row dimension mismatch for append");
528 return *this;
529 }
530 506
531 octave_idx_type nc_insert = nc; 507 octave_idx_type nc_insert = nc;
532 FloatMatrix retval (nr, nc + a.cols ()); 508 FloatMatrix retval (nr, nc + a.cols ());
533 retval.insert (*this, 0, 0); 509 retval.insert (*this, 0, 0);
534 retval.insert (a, 0, nc_insert); 510 retval.insert (a, 0, nc_insert);
539 FloatMatrix::stack (const FloatMatrix& a) const 515 FloatMatrix::stack (const FloatMatrix& a) const
540 { 516 {
541 octave_idx_type nr = rows (); 517 octave_idx_type nr = rows ();
542 octave_idx_type nc = cols (); 518 octave_idx_type nc = cols ();
543 if (nc != a.cols ()) 519 if (nc != a.cols ())
544 { 520 (*current_liboctave_error_handler) ("column dimension mismatch for stack");
545 (*current_liboctave_error_handler)
546 ("column dimension mismatch for stack");
547 return FloatMatrix ();
548 }
549 521
550 octave_idx_type nr_insert = nr; 522 octave_idx_type nr_insert = nr;
551 FloatMatrix retval (nr + a.rows (), nc); 523 FloatMatrix retval (nr + a.rows (), nc);
552 retval.insert (*this, 0, 0); 524 retval.insert (*this, 0, 0);
553 retval.insert (a, nr_insert, 0); 525 retval.insert (a, nr_insert, 0);
558 FloatMatrix::stack (const FloatRowVector& a) const 530 FloatMatrix::stack (const FloatRowVector& a) const
559 { 531 {
560 octave_idx_type nr = rows (); 532 octave_idx_type nr = rows ();
561 octave_idx_type nc = cols (); 533 octave_idx_type nc = cols ();
562 if (nc != a.numel ()) 534 if (nc != a.numel ())
563 { 535 (*current_liboctave_error_handler) ("column dimension mismatch for stack");
564 (*current_liboctave_error_handler)
565 ("column dimension mismatch for stack");
566 return FloatMatrix ();
567 }
568 536
569 octave_idx_type nr_insert = nr; 537 octave_idx_type nr_insert = nr;
570 FloatMatrix retval (nr + 1, nc); 538 FloatMatrix retval (nr + 1, nc);
571 retval.insert (*this, 0, 0); 539 retval.insert (*this, 0, 0);
572 retval.insert (a, nr_insert, 0); 540 retval.insert (a, nr_insert, 0);
577 FloatMatrix::stack (const FloatColumnVector& a) const 545 FloatMatrix::stack (const FloatColumnVector& a) const
578 { 546 {
579 octave_idx_type nr = rows (); 547 octave_idx_type nr = rows ();
580 octave_idx_type nc = cols (); 548 octave_idx_type nc = cols ();
581 if (nc != 1) 549 if (nc != 1)
582 { 550 (*current_liboctave_error_handler) ("column dimension mismatch for stack");
583 (*current_liboctave_error_handler)
584 ("column dimension mismatch for stack");
585 return FloatMatrix ();
586 }
587 551
588 octave_idx_type nr_insert = nr; 552 octave_idx_type nr_insert = nr;
589 FloatMatrix retval (nr + a.numel (), nc); 553 FloatMatrix retval (nr + a.numel (), nc);
590 retval.insert (*this, 0, 0); 554 retval.insert (*this, 0, 0);
591 retval.insert (a, nr_insert, 0); 555 retval.insert (a, nr_insert, 0);
596 FloatMatrix::stack (const FloatDiagMatrix& a) const 560 FloatMatrix::stack (const FloatDiagMatrix& a) const
597 { 561 {
598 octave_idx_type nr = rows (); 562 octave_idx_type nr = rows ();
599 octave_idx_type nc = cols (); 563 octave_idx_type nc = cols ();
600 if (nc != a.cols ()) 564 if (nc != a.cols ())
601 { 565 (*current_liboctave_error_handler) ("column dimension mismatch for stack");
602 (*current_liboctave_error_handler)
603 ("column dimension mismatch for stack");
604 return FloatMatrix ();
605 }
606 566
607 octave_idx_type nr_insert = nr; 567 octave_idx_type nr_insert = nr;
608 FloatMatrix retval (nr + a.rows (), nc); 568 FloatMatrix retval (nr + a.rows (), nc);
609 retval.insert (*this, 0, 0); 569 retval.insert (*this, 0, 0);
610 retval.insert (a, nr_insert, 0); 570 retval.insert (a, nr_insert, 0);
703 octave_idx_type nr = rows (); 663 octave_idx_type nr = rows ();
704 octave_idx_type nc = cols (); 664 octave_idx_type nc = cols ();
705 665
706 if (nr != nc || nr == 0 || nc == 0) 666 if (nr != nc || nr == 0 || nc == 0)
707 (*current_liboctave_error_handler) ("inverse requires square matrix"); 667 (*current_liboctave_error_handler) ("inverse requires square matrix");
708 else 668
709 { 669 int typ = mattype.type ();
710 int typ = mattype.type (); 670 char uplo = (typ == MatrixType::Lower ? 'L' : 'U');
711 char uplo = (typ == MatrixType::Lower ? 'L' : 'U'); 671 char udiag = 'N';
712 char udiag = 'N'; 672 retval = *this;
713 retval = *this; 673 float *tmp_data = retval.fortran_vec ();
714 float *tmp_data = retval.fortran_vec (); 674
715 675 F77_XFCN (strtri, STRTRI, (F77_CONST_CHAR_ARG2 (&uplo, 1),
716 F77_XFCN (strtri, STRTRI, (F77_CONST_CHAR_ARG2 (&uplo, 1), 676 F77_CONST_CHAR_ARG2 (&udiag, 1),
677 nr, tmp_data, nr, info
678 F77_CHAR_ARG_LEN (1)
679 F77_CHAR_ARG_LEN (1)));
680
681 // Throw-away extra info LAPACK gives so as to not change output.
682 rcon = 0.0;
683 if (info != 0)
684 info = -1;
685 else if (calc_cond)
686 {
687 octave_idx_type dtrcon_info = 0;
688 char job = '1';
689
690 OCTAVE_LOCAL_BUFFER (float, work, 3 * nr);
691 OCTAVE_LOCAL_BUFFER (octave_idx_type, iwork, nr);
692
693 F77_XFCN (strcon, STRCON, (F77_CONST_CHAR_ARG2 (&job, 1),
694 F77_CONST_CHAR_ARG2 (&uplo, 1),
717 F77_CONST_CHAR_ARG2 (&udiag, 1), 695 F77_CONST_CHAR_ARG2 (&udiag, 1),
718 nr, tmp_data, nr, info 696 nr, tmp_data, nr, rcon,
697 work, iwork, dtrcon_info
698 F77_CHAR_ARG_LEN (1)
719 F77_CHAR_ARG_LEN (1) 699 F77_CHAR_ARG_LEN (1)
720 F77_CHAR_ARG_LEN (1))); 700 F77_CHAR_ARG_LEN (1)));
721 701
722 // Throw-away extra info LAPACK gives so as to not change output. 702 if (dtrcon_info != 0)
723 rcon = 0.0;
724 if (info != 0)
725 info = -1; 703 info = -1;
726 else if (calc_cond) 704 }
727 { 705
728 octave_idx_type dtrcon_info = 0; 706 if (info == -1 && ! force)
729 char job = '1'; 707 retval = *this; // Restore matrix contents.
730
731 OCTAVE_LOCAL_BUFFER (float, work, 3 * nr);
732 OCTAVE_LOCAL_BUFFER (octave_idx_type, iwork, nr);
733
734 F77_XFCN (strcon, STRCON, (F77_CONST_CHAR_ARG2 (&job, 1),
735 F77_CONST_CHAR_ARG2 (&uplo, 1),
736 F77_CONST_CHAR_ARG2 (&udiag, 1),
737 nr, tmp_data, nr, rcon,
738 work, iwork, dtrcon_info
739 F77_CHAR_ARG_LEN (1)
740 F77_CHAR_ARG_LEN (1)
741 F77_CHAR_ARG_LEN (1)));
742
743 if (dtrcon_info != 0)
744 info = -1;
745 }
746
747 if (info == -1 && ! force)
748 retval = *this; // Restore matrix contents.
749 }
750 708
751 return retval; 709 return retval;
752 } 710 }
753 711
754 712
761 octave_idx_type nr = rows (); 719 octave_idx_type nr = rows ();
762 octave_idx_type nc = cols (); 720 octave_idx_type nc = cols ();
763 721
764 if (nr != nc || nr == 0 || nc == 0) 722 if (nr != nc || nr == 0 || nc == 0)
765 (*current_liboctave_error_handler) ("inverse requires square matrix"); 723 (*current_liboctave_error_handler) ("inverse requires square matrix");
724
725 Array<octave_idx_type> ipvt (dim_vector (nr, 1));
726 octave_idx_type *pipvt = ipvt.fortran_vec ();
727
728 retval = *this;
729 float *tmp_data = retval.fortran_vec ();
730
731 Array<float> z(dim_vector (1, 1));
732 octave_idx_type lwork = -1;
733
734 // Query the optimum work array size.
735 F77_XFCN (sgetri, SGETRI, (nc, tmp_data, nr, pipvt,
736 z.fortran_vec (), lwork, info));
737
738 lwork = static_cast<octave_idx_type> (z(0));
739 lwork = (lwork < 2 *nc ? 2*nc : lwork);
740 z.resize (dim_vector (lwork, 1));
741 float *pz = z.fortran_vec ();
742
743 info = 0;
744
745 // Calculate the norm of the matrix, for later use.
746 float anorm = 0;
747 if (calc_cond)
748 anorm = retval.abs ().sum ().row (static_cast<octave_idx_type>(0))
749 .max ();
750
751 F77_XFCN (sgetrf, SGETRF, (nc, nc, tmp_data, nr, pipvt, info));
752
753 // Throw-away extra info LAPACK gives so as to not change output.
754 rcon = 0.0;
755 if (info != 0)
756 info = -1;
757 else if (calc_cond)
758 {
759 octave_idx_type dgecon_info = 0;
760
761 // Now calculate the condition number for non-singular matrix.
762 char job = '1';
763 Array<octave_idx_type> iz (dim_vector (nc, 1));
764 octave_idx_type *piz = iz.fortran_vec ();
765 F77_XFCN (sgecon, SGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
766 nc, tmp_data, nr, anorm,
767 rcon, pz, piz, dgecon_info
768 F77_CHAR_ARG_LEN (1)));
769
770 if (dgecon_info != 0)
771 info = -1;
772 }
773
774 if (info == -1 && ! force)
775 retval = *this; // Restore matrix contents.
766 else 776 else
767 { 777 {
768 Array<octave_idx_type> ipvt (dim_vector (nr, 1)); 778 octave_idx_type dgetri_info = 0;
769 octave_idx_type *pipvt = ipvt.fortran_vec (); 779
770
771 retval = *this;
772 float *tmp_data = retval.fortran_vec ();
773
774 Array<float> z(dim_vector (1, 1));
775 octave_idx_type lwork = -1;
776
777 // Query the optimum work array size.
778 F77_XFCN (sgetri, SGETRI, (nc, tmp_data, nr, pipvt, 780 F77_XFCN (sgetri, SGETRI, (nc, tmp_data, nr, pipvt,
779 z.fortran_vec (), lwork, info)); 781 pz, lwork, dgetri_info));
780 782
781 lwork = static_cast<octave_idx_type> (z(0)); 783 if (dgetri_info != 0)
782 lwork = (lwork < 2 *nc ? 2*nc : lwork);
783 z.resize (dim_vector (lwork, 1));
784 float *pz = z.fortran_vec ();
785
786 info = 0;
787
788 // Calculate the norm of the matrix, for later use.
789 float anorm = 0;
790 if (calc_cond)
791 anorm = retval.abs ().sum ().row (static_cast<octave_idx_type>(0))
792 .max ();
793
794 F77_XFCN (sgetrf, SGETRF, (nc, nc, tmp_data, nr, pipvt, info));
795
796 // Throw-away extra info LAPACK gives so as to not change output.
797 rcon = 0.0;
798 if (info != 0)
799 info = -1; 784 info = -1;
800 else if (calc_cond) 785 }
801 { 786
802 octave_idx_type dgecon_info = 0; 787 if (info != 0)
803 788 mattype.mark_as_rectangular ();
804 // Now calculate the condition number for non-singular matrix.
805 char job = '1';
806 Array<octave_idx_type> iz (dim_vector (nc, 1));
807 octave_idx_type *piz = iz.fortran_vec ();
808 F77_XFCN (sgecon, SGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
809 nc, tmp_data, nr, anorm,
810 rcon, pz, piz, dgecon_info
811 F77_CHAR_ARG_LEN (1)));
812
813 if (dgecon_info != 0)
814 info = -1;
815 }
816
817 if (info == -1 && ! force)
818 retval = *this; // Restore matrix contents.
819 else
820 {
821 octave_idx_type dgetri_info = 0;
822
823 F77_XFCN (sgetri, SGETRI, (nc, tmp_data, nr, pipvt,
824 pz, lwork, dgetri_info));
825
826 if (dgetri_info != 0)
827 info = -1;
828 }
829
830 if (info != 0)
831 mattype.mark_as_rectangular ();
832 }
833 789
834 return retval; 790 return retval;
835 } 791 }
836 792
837 FloatMatrix 793 FloatMatrix
1275 octave_idx_type nr = rows (); 1231 octave_idx_type nr = rows ();
1276 octave_idx_type nc = cols (); 1232 octave_idx_type nc = cols ();
1277 1233
1278 if (nr != nc) 1234 if (nr != nc)
1279 (*current_liboctave_error_handler) ("matrix must be square"); 1235 (*current_liboctave_error_handler) ("matrix must be square");
1280 else 1236
1281 { 1237 volatile int typ = mattype.type ();
1282 volatile int typ = mattype.type (); 1238
1283 1239 // Even though the matrix is marked as singular (Rectangular), we may
1284 // Even though the matrix is marked as singular (Rectangular), we may 1240 // still get a useful number from the LU factorization, because it always
1285 // still get a useful number from the LU factorization, because it always 1241 // completes.
1286 // completes. 1242
1287 1243 if (typ == MatrixType::Unknown)
1288 if (typ == MatrixType::Unknown) 1244 typ = mattype.type (*this);
1289 typ = mattype.type (*this); 1245 else if (typ == MatrixType::Rectangular)
1290 else if (typ == MatrixType::Rectangular) 1246 typ = MatrixType::Full;
1291 typ = MatrixType::Full; 1247
1292 1248 if (typ == MatrixType::Lower || typ == MatrixType::Upper)
1293 if (typ == MatrixType::Lower || typ == MatrixType::Upper) 1249 {
1250 for (octave_idx_type i = 0; i < nc; i++)
1251 retval *= elem (i,i);
1252 }
1253 else if (typ == MatrixType::Hermitian)
1254 {
1255 FloatMatrix atmp = *this;
1256 float *tmp_data = atmp.fortran_vec ();
1257
1258 float anorm = 0;
1259 if (calc_cond) anorm = xnorm (*this, 1);
1260
1261
1262 char job = 'L';
1263 F77_XFCN (spotrf, SPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr,
1264 tmp_data, nr, info
1265 F77_CHAR_ARG_LEN (1)));
1266
1267 if (info != 0)
1294 { 1268 {
1269 rcon = 0.0;
1270 mattype.mark_as_unsymmetric ();
1271 typ = MatrixType::Full;
1272 }
1273 else
1274 {
1275 Array<float> z (dim_vector (3 * nc, 1));
1276 float *pz = z.fortran_vec ();
1277 Array<octave_idx_type> iz (dim_vector (nc, 1));
1278 octave_idx_type *piz = iz.fortran_vec ();
1279
1280 F77_XFCN (spocon, SPOCON, (F77_CONST_CHAR_ARG2 (&job, 1),
1281 nr, tmp_data, nr, anorm,
1282 rcon, pz, piz, info
1283 F77_CHAR_ARG_LEN (1)));
1284
1285 if (info != 0)
1286 rcon = 0.0;
1287
1295 for (octave_idx_type i = 0; i < nc; i++) 1288 for (octave_idx_type i = 0; i < nc; i++)
1296 retval *= elem (i,i); 1289 retval *= atmp (i,i);
1290
1291 retval = retval.square ();
1297 } 1292 }
1298 else if (typ == MatrixType::Hermitian) 1293 }
1294 else if (typ != MatrixType::Full)
1295 (*current_liboctave_error_handler) ("det: invalid dense matrix type");
1296
1297 if (typ == MatrixType::Full)
1298 {
1299 Array<octave_idx_type> ipvt (dim_vector (nr, 1));
1300 octave_idx_type *pipvt = ipvt.fortran_vec ();
1301
1302 FloatMatrix atmp = *this;
1303 float *tmp_data = atmp.fortran_vec ();
1304
1305 info = 0;
1306
1307 // Calculate the norm of the matrix, for later use.
1308 float anorm = 0;
1309 if (calc_cond) anorm = xnorm (*this, 1);
1310
1311 F77_XFCN (sgetrf, SGETRF, (nr, nr, tmp_data, nr, pipvt, info));
1312
1313 // Throw-away extra info LAPACK gives so as to not change output.
1314 rcon = 0.0;
1315 if (info != 0)
1299 { 1316 {
1300 FloatMatrix atmp = *this; 1317 info = -1;
1301 float *tmp_data = atmp.fortran_vec (); 1318 retval = FloatDET ();
1302 1319 }
1303 float anorm = 0; 1320 else
1304 if (calc_cond) anorm = xnorm (*this, 1); 1321 {
1305 1322 if (calc_cond)
1306
1307 char job = 'L';
1308 F77_XFCN (spotrf, SPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr,
1309 tmp_data, nr, info
1310 F77_CHAR_ARG_LEN (1)));
1311
1312 if (info != 0)
1313 { 1323 {
1314 rcon = 0.0; 1324 // Now calc the condition number for non-singular matrix.
1315 mattype.mark_as_unsymmetric (); 1325 char job = '1';
1316 typ = MatrixType::Full; 1326 Array<float> z (dim_vector (4 * nc, 1));
1317 }
1318 else
1319 {
1320 Array<float> z (dim_vector (3 * nc, 1));
1321 float *pz = z.fortran_vec (); 1327 float *pz = z.fortran_vec ();
1322 Array<octave_idx_type> iz (dim_vector (nc, 1)); 1328 Array<octave_idx_type> iz (dim_vector (nc, 1));
1323 octave_idx_type *piz = iz.fortran_vec (); 1329 octave_idx_type *piz = iz.fortran_vec ();
1324 1330
1325 F77_XFCN (spocon, SPOCON, (F77_CONST_CHAR_ARG2 (&job, 1), 1331 F77_XFCN (sgecon, SGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
1326 nr, tmp_data, nr, anorm, 1332 nc, tmp_data, nr, anorm,
1327 rcon, pz, piz, info 1333 rcon, pz, piz, info
1328 F77_CHAR_ARG_LEN (1))); 1334 F77_CHAR_ARG_LEN (1)));
1329
1330 if (info != 0)
1331 rcon = 0.0;
1332
1333 for (octave_idx_type i = 0; i < nc; i++)
1334 retval *= atmp (i,i);
1335
1336 retval = retval.square ();
1337 } 1335 }
1338 } 1336
1339 else if (typ != MatrixType::Full)
1340 (*current_liboctave_error_handler) ("det: invalid dense matrix type");
1341
1342 if (typ == MatrixType::Full)
1343 {
1344 Array<octave_idx_type> ipvt (dim_vector (nr, 1));
1345 octave_idx_type *pipvt = ipvt.fortran_vec ();
1346
1347 FloatMatrix atmp = *this;
1348 float *tmp_data = atmp.fortran_vec ();
1349
1350 info = 0;
1351
1352 // Calculate the norm of the matrix, for later use.
1353 float anorm = 0;
1354 if (calc_cond) anorm = xnorm (*this, 1);
1355
1356 F77_XFCN (sgetrf, SGETRF, (nr, nr, tmp_data, nr, pipvt, info));
1357
1358 // Throw-away extra info LAPACK gives so as to not change output.
1359 rcon = 0.0;
1360 if (info != 0) 1337 if (info != 0)
1361 { 1338 {
1362 info = -1; 1339 info = -1;
1363 retval = FloatDET (); 1340 retval = FloatDET ();
1364 } 1341 }
1365 else 1342 else
1366 { 1343 {
1367 if (calc_cond) 1344 for (octave_idx_type i = 0; i < nc; i++)
1368 { 1345 {
1369 // Now calc the condition number for non-singular matrix. 1346 float c = atmp(i,i);
1370 char job = '1'; 1347 retval *= (ipvt(i) != (i+1)) ? -c : c;
1371 Array<float> z (dim_vector (4 * nc, 1));
1372 float *pz = z.fortran_vec ();
1373 Array<octave_idx_type> iz (dim_vector (nc, 1));
1374 octave_idx_type *piz = iz.fortran_vec ();
1375
1376 F77_XFCN (sgecon, SGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
1377 nc, tmp_data, nr, anorm,
1378 rcon, pz, piz, info
1379 F77_CHAR_ARG_LEN (1)));
1380 }
1381
1382 if (info != 0)
1383 {
1384 info = -1;
1385 retval = FloatDET ();
1386 }
1387 else
1388 {
1389 for (octave_idx_type i = 0; i < nc; i++)
1390 {
1391 float c = atmp(i,i);
1392 retval *= (ipvt(i) != (i+1)) ? -c : c;
1393 }
1394 } 1348 }
1395 } 1349 }
1396 } 1350 }
1397 } 1351 }
1398 1352
1413 octave_idx_type nr = rows (); 1367 octave_idx_type nr = rows ();
1414 octave_idx_type nc = cols (); 1368 octave_idx_type nc = cols ();
1415 1369
1416 if (nr != nc) 1370 if (nr != nc)
1417 (*current_liboctave_error_handler) ("matrix must be square"); 1371 (*current_liboctave_error_handler) ("matrix must be square");
1418 else if (nr == 0 || nc == 0) 1372
1373 if (nr == 0 || nc == 0)
1419 rcon = octave_Inf; 1374 rcon = octave_Inf;
1420 else 1375 else
1421 { 1376 {
1422 volatile int typ = mattype.type (); 1377 volatile int typ = mattype.type ();
1423 1378
1448 F77_CHAR_ARG_LEN (1))); 1403 F77_CHAR_ARG_LEN (1)));
1449 1404
1450 if (info != 0) 1405 if (info != 0)
1451 rcon = 0.0; 1406 rcon = 0.0;
1452 } 1407 }
1453 else if (typ == MatrixType::Permuted_Upper) 1408 else if (typ == MatrixType::Permuted_Upper)
1454 (*current_liboctave_error_handler) 1409 (*current_liboctave_error_handler)
1455 ("permuted triangular matrix not implemented"); 1410 ("permuted triangular matrix not implemented");
1456 else if (typ == MatrixType::Lower) 1411 else if (typ == MatrixType::Lower)
1457 { 1412 {
1458 const float *tmp_data = fortran_vec (); 1413 const float *tmp_data = fortran_vec ();
1581 octave_idx_type nc = cols (); 1536 octave_idx_type nc = cols ();
1582 1537
1583 if (nr != b.rows ()) 1538 if (nr != b.rows ())
1584 (*current_liboctave_error_handler) 1539 (*current_liboctave_error_handler)
1585 ("matrix dimension mismatch solution of linear equations"); 1540 ("matrix dimension mismatch solution of linear equations");
1586 else if (nr == 0 || nc == 0 || b.cols () == 0) 1541
1542 if (nr == 0 || nc == 0 || b.cols () == 0)
1587 retval = FloatMatrix (nc, b.cols (), 0.0); 1543 retval = FloatMatrix (nc, b.cols (), 0.0);
1588 else 1544 else
1589 { 1545 {
1590 volatile int typ = mattype.type (); 1546 volatile int typ = mattype.type ();
1591 1547
1594 octave_idx_type b_nc = b.cols (); 1550 octave_idx_type b_nc = b.cols ();
1595 rcon = 1.; 1551 rcon = 1.;
1596 info = 0; 1552 info = 0;
1597 1553
1598 if (typ == MatrixType::Permuted_Upper) 1554 if (typ == MatrixType::Permuted_Upper)
1599 { 1555 (*current_liboctave_error_handler)
1600 (*current_liboctave_error_handler) 1556 ("permuted triangular matrix not implemented");
1601 ("permuted triangular matrix not implemented");
1602 }
1603 else 1557 else
1604 { 1558 {
1605 const float *tmp_data = fortran_vec (); 1559 const float *tmp_data = fortran_vec ();
1606 1560
1607 retval = b; 1561 retval = b;
1677 octave_idx_type nc = cols (); 1631 octave_idx_type nc = cols ();
1678 1632
1679 if (nr != b.rows ()) 1633 if (nr != b.rows ())
1680 (*current_liboctave_error_handler) 1634 (*current_liboctave_error_handler)
1681 ("matrix dimension mismatch solution of linear equations"); 1635 ("matrix dimension mismatch solution of linear equations");
1682 else if (nr == 0 || nc == 0 || b.cols () == 0) 1636
1637 if (nr == 0 || nc == 0 || b.cols () == 0)
1683 retval = FloatMatrix (nc, b.cols (), 0.0); 1638 retval = FloatMatrix (nc, b.cols (), 0.0);
1684 else 1639 else
1685 { 1640 {
1686 volatile int typ = mattype.type (); 1641 volatile int typ = mattype.type ();
1687 1642
1690 octave_idx_type b_nc = b.cols (); 1645 octave_idx_type b_nc = b.cols ();
1691 rcon = 1.; 1646 rcon = 1.;
1692 info = 0; 1647 info = 0;
1693 1648
1694 if (typ == MatrixType::Permuted_Lower) 1649 if (typ == MatrixType::Permuted_Lower)
1695 { 1650 (*current_liboctave_error_handler)
1696 (*current_liboctave_error_handler) 1651 ("permuted triangular matrix not implemented");
1697 ("permuted triangular matrix not implemented");
1698 }
1699 else 1652 else
1700 { 1653 {
1701 const float *tmp_data = fortran_vec (); 1654 const float *tmp_data = fortran_vec ();
1702 1655
1703 retval = b; 1656 retval = b;
1772 octave_idx_type nc = cols (); 1725 octave_idx_type nc = cols ();
1773 1726
1774 if (nr != nc || nr != b.rows ()) 1727 if (nr != nc || nr != b.rows ())
1775 (*current_liboctave_error_handler) 1728 (*current_liboctave_error_handler)
1776 ("matrix dimension mismatch solution of linear equations"); 1729 ("matrix dimension mismatch solution of linear equations");
1777 else if (nr == 0 || b.cols () == 0) 1730
1731 if (nr == 0 || b.cols () == 0)
1778 retval = FloatMatrix (nc, b.cols (), 0.0); 1732 retval = FloatMatrix (nc, b.cols (), 0.0);
1779 else 1733 else
1780 { 1734 {
1781 volatile int typ = mattype.type (); 1735 volatile int typ = mattype.type ();
1782 1736
1985 return transpose ().solve (mattype, b, info, rcon, sing_handler, 1939 return transpose ().solve (mattype, b, info, rcon, sing_handler,
1986 singular_fallback); 1940 singular_fallback);
1987 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian) 1941 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
1988 retval = fsolve (mattype, b, info, rcon, sing_handler, true); 1942 retval = fsolve (mattype, b, info, rcon, sing_handler, true);
1989 else if (typ != MatrixType::Rectangular) 1943 else if (typ != MatrixType::Rectangular)
1990 { 1944 (*current_liboctave_error_handler) ("unknown matrix type");
1991 (*current_liboctave_error_handler) ("unknown matrix type");
1992 return FloatMatrix ();
1993 }
1994 1945
1995 // Rectangular or one of the above solvers flags a singular matrix 1946 // Rectangular or one of the above solvers flags a singular matrix
1996 if (singular_fallback && mattype.type () == MatrixType::Rectangular) 1947 if (singular_fallback && mattype.type () == MatrixType::Rectangular)
1997 { 1948 {
1998 octave_idx_type rank; 1949 octave_idx_type rank;
2297 octave_idx_type n = cols (); 2248 octave_idx_type n = cols ();
2298 2249
2299 if (m != b.rows ()) 2250 if (m != b.rows ())
2300 (*current_liboctave_error_handler) 2251 (*current_liboctave_error_handler)
2301 ("matrix dimension mismatch solution of linear equations"); 2252 ("matrix dimension mismatch solution of linear equations");
2302 else if (m == 0 || n == 0 || b.cols () == 0) 2253
2254 if (m == 0 || n == 0 || b.cols () == 0)
2303 retval = FloatMatrix (n, b.cols (), 0.0); 2255 retval = FloatMatrix (n, b.cols (), 0.0);
2304 else 2256 else
2305 { 2257 {
2306 volatile octave_idx_type minmn = (m < n ? m : n); 2258 volatile octave_idx_type minmn = (m < n ? m : n);
2307 octave_idx_type maxmn = m > n ? m : n; 2259 octave_idx_type maxmn = m > n ? m : n;
2492 octave_idx_type n = cols (); 2444 octave_idx_type n = cols ();
2493 2445
2494 if (m != b.numel ()) 2446 if (m != b.numel ())
2495 (*current_liboctave_error_handler) 2447 (*current_liboctave_error_handler)
2496 ("matrix dimension mismatch solution of linear equations"); 2448 ("matrix dimension mismatch solution of linear equations");
2497 else if (m == 0 || n == 0) 2449
2450 if (m == 0 || n == 0)
2498 retval = FloatColumnVector (n, 0.0); 2451 retval = FloatColumnVector (n, 0.0);
2499 else 2452 else
2500 { 2453 {
2501 volatile octave_idx_type minmn = (m < n ? m : n); 2454 volatile octave_idx_type minmn = (m < n ? m : n);
2502 octave_idx_type maxmn = m > n ? m : n; 2455 octave_idx_type maxmn = m > n ? m : n;
2726 octave_idx_type nc = cols (); 2679 octave_idx_type nc = cols ();
2727 2680
2728 if (nr == 1 || nc == 1) 2681 if (nr == 1 || nc == 1)
2729 retval = FloatDiagMatrix (*this, m, n); 2682 retval = FloatDiagMatrix (*this, m, n);
2730 else 2683 else
2731 (*current_liboctave_error_handler) 2684 (*current_liboctave_error_handler) ("diag: expecting vector argument");
2732 ("diag: expecting vector argument");
2733 2685
2734 return retval; 2686 return retval;
2735 } 2687 }
2736 2688
2737 FloatColumnVector 2689 FloatColumnVector
3233 { 3185 {
3234 octave_idx_type nr = a.rows (); 3186 octave_idx_type nr = a.rows ();
3235 octave_idx_type nc = a.columns (); 3187 octave_idx_type nc = a.columns ();
3236 3188
3237 if (nr != b.rows () || nc != b.columns ()) 3189 if (nr != b.rows () || nc != b.columns ())
3238 { 3190 (*current_liboctave_error_handler)
3239 (*current_liboctave_error_handler) 3191 ("two-arg min requires same size arguments");
3240 ("two-arg min requires same size arguments");
3241 return FloatMatrix ();
3242 }
3243 3192
3244 EMPTY_RETURN_CHECK (FloatMatrix); 3193 EMPTY_RETURN_CHECK (FloatMatrix);
3245 3194
3246 FloatMatrix result (nr, nc); 3195 FloatMatrix result (nr, nc);
3247 3196
3300 { 3249 {
3301 octave_idx_type nr = a.rows (); 3250 octave_idx_type nr = a.rows ();
3302 octave_idx_type nc = a.columns (); 3251 octave_idx_type nc = a.columns ();
3303 3252
3304 if (nr != b.rows () || nc != b.columns ()) 3253 if (nr != b.rows () || nc != b.columns ())
3305 { 3254 (*current_liboctave_error_handler)
3306 (*current_liboctave_error_handler) 3255 ("two-arg max requires same size arguments");
3307 ("two-arg max requires same size arguments");
3308 return FloatMatrix ();
3309 }
3310 3256
3311 EMPTY_RETURN_CHECK (FloatMatrix); 3257 EMPTY_RETURN_CHECK (FloatMatrix);
3312 3258
3313 FloatMatrix result (nr, nc); 3259 FloatMatrix result (nr, nc);
3314 3260