Mercurial > hg > octave-nkf

comparison liboctave/fMatrix.cc @ 15382:197774b411ec stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
rcond: use new copy of data for full factorization if positive definite cholesky factorization fails (bug #37336) * dMatrix.cc (Matrix::rcond): Don't reuse modified matrix data if positive definite cholesky factorization was attempted but fails. * CMatrix.cc (ComplexMatrix::rcond): Likewise. * fMatrix.cc (FloatMatrix::rcond): Likewise. * fCMatrix.cc (FloatComplexMatrix::rcond): Likewise. * rcond.cc: New tests.
author John W. Eaton <jwe@octave.org>
date Thu, 13 Sep 2012 15:14:47 -0400
parents 72c96de7a403
children
comparison
equal deleted inserted replaced
15304:4663cc835c65 15382:197774b411ec
1452 (*current_liboctave_error_handler) 1452 (*current_liboctave_error_handler)
1453 ("permuted triangular matrix not implemented"); 1453 ("permuted triangular matrix not implemented");
1454 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian) 1454 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
1455 { 1455 {
1456 float anorm = -1.0; 1456 float anorm = -1.0;
1457 FloatMatrix atmp = *this;
1458 float *tmp_data = atmp.fortran_vec ();
1459 1457
1460 if (typ == MatrixType::Hermitian) 1458 if (typ == MatrixType::Hermitian)
1461 { 1459 {
1462 octave_idx_type info = 0; 1460 octave_idx_type info = 0;
1463 char job = 'L'; 1461 char job = 'L';
1462
1463 FloatMatrix atmp = *this;
1464 float *tmp_data = atmp.fortran_vec ();
1465
1464 anorm = atmp.abs().sum(). 1466 anorm = atmp.abs().sum().
1465 row(static_cast<octave_idx_type>(0)).max(); 1467 row(static_cast<octave_idx_type>(0)).max();
1466 1468
1467 F77_XFCN (spotrf, SPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr, 1469 F77_XFCN (spotrf, SPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr,
1468 tmp_data, nr, info 1470 tmp_data, nr, info
1492 } 1494 }
1493 1495
1494 if (typ == MatrixType::Full) 1496 if (typ == MatrixType::Full)
1495 { 1497 {
1496 octave_idx_type info = 0; 1498 octave_idx_type info = 0;
1499
1500 FloatMatrix atmp = *this;
1501 float *tmp_data = atmp.fortran_vec ();
1497 1502
1498 Array<octave_idx_type> ipvt (dim_vector (nr, 1)); 1503 Array<octave_idx_type> ipvt (dim_vector (nr, 1));
1499 octave_idx_type *pipvt = ipvt.fortran_vec (); 1504 octave_idx_type *pipvt = ipvt.fortran_vec ();
1500 1505
1501 if(anorm < 0.) 1506 if(anorm < 0.)
1758 1763
1759 if (typ == MatrixType::Hermitian) 1764 if (typ == MatrixType::Hermitian)
1760 { 1765 {
1761 info = 0; 1766 info = 0;
1762 char job = 'L'; 1767 char job = 'L';
1768
1763 FloatMatrix atmp = *this; 1769 FloatMatrix atmp = *this;
1764 float *tmp_data = atmp.fortran_vec (); 1770 float *tmp_data = atmp.fortran_vec ();
1771
1765 anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max(); 1772 anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max();
1766 1773
1767 F77_XFCN (spotrf, SPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr, 1774 F77_XFCN (spotrf, SPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr,
1768 tmp_data, nr, info 1775 tmp_data, nr, info
1769 F77_CHAR_ARG_LEN (1))); 1776 F77_CHAR_ARG_LEN (1)));
1836 Array<octave_idx_type> ipvt (dim_vector (nr, 1)); 1843 Array<octave_idx_type> ipvt (dim_vector (nr, 1));
1837 octave_idx_type *pipvt = ipvt.fortran_vec (); 1844 octave_idx_type *pipvt = ipvt.fortran_vec ();
1838 1845
1839 FloatMatrix atmp = *this; 1846 FloatMatrix atmp = *this;
1840 float *tmp_data = atmp.fortran_vec (); 1847 float *tmp_data = atmp.fortran_vec ();
1848
1841 if(anorm < 0.) 1849 if(anorm < 0.)
1842 anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max(); 1850 anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max();
1843 1851
1844 Array<float> z (dim_vector (4 * nc, 1)); 1852 Array<float> z (dim_vector (4 * nc, 1));
1845 float *pz = z.fortran_vec (); 1853 float *pz = z.fortran_vec ();