Mercurial > octave-nkf
changeset 7482:29980c6b8604
don't check f77_exception_encountered
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Thu, 14 Feb 2008 21:57:50 -0500 |
| parents | 78f3811155f7 |
| children | fb66330b2608 |
| files | ChangeLog examples/fortdemo.cc liboctave/CColVector.cc liboctave/CMatrix.cc liboctave/CRowVector.cc liboctave/CSparse.cc liboctave/ChangeLog liboctave/CmplxAEPBAL.cc liboctave/CmplxCHOL.cc liboctave/CmplxHESS.cc liboctave/CmplxLU.cc liboctave/CmplxQR.cc liboctave/CmplxQRP.cc liboctave/CmplxSCHUR.cc liboctave/CmplxSVD.cc liboctave/DASPK.cc liboctave/DASRT.cc liboctave/DASSL.cc liboctave/EIG.cc liboctave/LSODE.cc liboctave/NLEqn.cc liboctave/Quad.cc liboctave/dColVector.cc liboctave/dMatrix.cc liboctave/dRowVector.cc liboctave/dSparse.cc liboctave/dbleAEPBAL.cc liboctave/dbleCHOL.cc liboctave/dbleHESS.cc liboctave/dbleLU.cc liboctave/dbleQR.cc liboctave/dbleQRP.cc liboctave/dbleSCHUR.cc liboctave/dbleSVD.cc src/ChangeLog src/DLD-FUNCTIONS/balance.cc src/DLD-FUNCTIONS/qz.cc |
| diffstat | 37 files changed, 2615 insertions(+), 3544 deletions(-) [+] |
line wrap: on
line diff
--- a/ChangeLog Thu Feb 14 17:14:23 2008 -0500 +++ b/ChangeLog Thu Feb 14 21:57:50 2008 -0500 @@ -1,3 +1,7 @@ +2008-02-14 John W. Eaton <jwe@octave.org> + + * examples/fortdemo.cc: Don't check f77_exception_encountered. + 2008-02-06 John W. Eaton <jwe@octave.org> * examples/Makefile.in (octave.desktop):
--- a/examples/fortdemo.cc Thu Feb 14 17:14:23 2008 -0500 +++ b/examples/fortdemo.cc Thu Feb 14 21:57:50 2008 -0500 @@ -50,13 +50,8 @@ F77_XFCN (fortsub, FORTSUB, (na, av, ctmp F77_CHAR_ARG_LEN (128))); - if (f77_exception_encountered) - error ("fortdemo: error in fortran"); - else - { - retval(1) = std::string (ctmp); - retval(0) = a; - } + retval(1) = std::string (ctmp); + retval(0) = a; } } return retval;
--- a/liboctave/CColVector.cc Thu Feb 14 17:14:23 2008 -0500 +++ b/liboctave/CColVector.cc Thu Feb 14 21:57:50 2008 -0500 @@ -353,10 +353,6 @@ nr, nc, 1.0, m.data (), ld, a.data (), 1, 0.0, y, 1 F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zgemv"); } }
--- a/liboctave/CMatrix.cc Thu Feb 14 17:14:23 2008 -0500 +++ b/liboctave/CMatrix.cc Thu Feb 14 21:57:50 2008 -0500 @@ -1043,38 +1043,29 @@ F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in ztrtri"); - else + // Throw-away extra info LAPACK gives so as to not change output. + rcond = 0.0; + if (info != 0) + info = -1; + else if (calc_cond) { - // Throw-away extra info LAPACK gives so as to not change output. - rcond = 0.0; - if (info != 0) + octave_idx_type ztrcon_info = 0; + char job = '1'; + + OCTAVE_LOCAL_BUFFER (Complex, cwork, 2*nr); + OCTAVE_LOCAL_BUFFER (double, rwork, nr); + + F77_XFCN (ztrcon, ZTRCON, (F77_CONST_CHAR_ARG2 (&job, 1), + F77_CONST_CHAR_ARG2 (&uplo, 1), + F77_CONST_CHAR_ARG2 (&udiag, 1), + nr, tmp_data, nr, rcond, + cwork, rwork, ztrcon_info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (ztrcon_info != 0) info = -1; - else if (calc_cond) - { - octave_idx_type ztrcon_info = 0; - char job = '1'; - - OCTAVE_LOCAL_BUFFER (Complex, cwork, 2*nr); - OCTAVE_LOCAL_BUFFER (double, rwork, nr); - - F77_XFCN (ztrcon, ZTRCON, (F77_CONST_CHAR_ARG2 (&job, 1), - F77_CONST_CHAR_ARG2 (&uplo, 1), - F77_CONST_CHAR_ARG2 (&udiag, 1), - nr, tmp_data, nr, rcond, - cwork, rwork, ztrcon_info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in ztrcon"); - - if (ztrcon_info != 0) - info = -1; - } } if (info == -1 && ! force) @@ -1111,13 +1102,6 @@ F77_XFCN (zgetri, ZGETRI, (nc, tmp_data, nr, pipvt, z.fortran_vec (), lwork, info)); - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("unrecoverable error in zgetri"); - return retval; - } - lwork = static_cast<octave_idx_type> (std::real(z(0))); lwork = (lwork < 2 *nc ? 2*nc : lwork); z.resize (lwork); @@ -1132,50 +1116,37 @@ F77_XFCN (zgetrf, ZGETRF, (nc, nc, tmp_data, nr, pipvt, info)); - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in zgetrf"); + // Throw-away extra info LAPACK gives so as to not change output. + rcond = 0.0; + if (info != 0) + info = -1; + else if (calc_cond) + { + // Now calculate the condition number for non-singular matrix. + octave_idx_type zgecon_info = 0; + char job = '1'; + Array<double> rz (2 * nc); + double *prz = rz.fortran_vec (); + F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1), + nc, tmp_data, nr, anorm, + rcond, pz, prz, zgecon_info + F77_CHAR_ARG_LEN (1))); + + if (zgecon_info != 0) + info = -1; + } + + if (info == -1 && ! force) + retval = *this; // Restore contents. else { - // Throw-away extra info LAPACK gives so as to not change output. - rcond = 0.0; - if (info != 0) + octave_idx_type zgetri_info = 0; + + F77_XFCN (zgetri, ZGETRI, (nc, tmp_data, nr, pipvt, + pz, lwork, zgetri_info)); + + if (zgetri_info != 0) info = -1; - else if (calc_cond) - { - // Now calculate the condition number for non-singular matrix. - octave_idx_type zgecon_info = 0; - char job = '1'; - Array<double> rz (2 * nc); - double *prz = rz.fortran_vec (); - F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1), - nc, tmp_data, nr, anorm, - rcond, pz, prz, zgecon_info - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zgecon"); - - if (zgecon_info != 0) - info = -1; - } - - if (info == -1 && ! force) - retval = *this; // Restore contents. - else - { - octave_idx_type zgetri_info = 0; - - F77_XFCN (zgetri, ZGETRI, (nc, tmp_data, nr, pipvt, - pz, lwork, zgetri_info)); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zgetri"); - - if (zgetri_info != 0) - info = -1; - } } if (info != 0) @@ -1621,12 +1592,30 @@ F77_XFCN (zgetrf, ZGETRF, (nr, nc, tmp_data, nr, pipvt, info)); - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in zgetrf"); - else + // Throw-away extra info LAPACK gives so as to not change output. + rcond = 0.0; + if (info != 0) + { + info = -1; + retval = ComplexDET (); + } + else { - // Throw-away extra info LAPACK gives so as to not change output. - rcond = 0.0; + if (calc_cond) + { + // Now calc the condition number for non-singular matrix. + char job = '1'; + Array<Complex> z (2*nr); + Complex *pz = z.fortran_vec (); + Array<double> rz (2*nr); + double *prz = rz.fortran_vec (); + + F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1), + nc, tmp_data, nr, anorm, + rcond, pz, prz, info + F77_CHAR_ARG_LEN (1))); + } + if (info != 0) { info = -1; @@ -1634,60 +1623,33 @@ } else { - if (calc_cond) - { - // Now calc the condition number for non-singular matrix. - char job = '1'; - Array<Complex> z (2*nr); - Complex *pz = z.fortran_vec (); - Array<double> rz (2*nr); - double *prz = rz.fortran_vec (); - - F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1), - nc, tmp_data, nr, anorm, - rcond, pz, prz, info - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zgecon"); - } - - if (info != 0) - { - info = -1; - retval = ComplexDET (); - } - else + Complex c = 1.0; + int e = 0; + + for (octave_idx_type i = 0; i < nc; i++) { - Complex c = 1.0; - int e = 0; - - for (octave_idx_type i = 0; i < nc; i++) + if (ipvt(i) != (i+1)) + c = -c; + + c *= atmp(i,i); + + if (c == 0.0) + break; + + while (std::abs(c) < 0.5) { - if (ipvt(i) != (i+1)) - c = -c; - - c *= atmp(i,i); - - if (c == 0.0) - break; - - while (std::abs(c) < 0.5) - { - c *= 2.0; - e--; - } - - while (std::abs(c) >= 2.0) - { - c /= 2.0; - e++; - } + c *= 2.0; + e--; } - retval = ComplexDET (c, e); + while (std::abs(c) >= 2.0) + { + c /= 2.0; + e++; + } } + + retval = ComplexDET (c, e); } } } @@ -1751,10 +1713,6 @@ F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in ztrcon"); - if (info != 0) info = -2; @@ -1790,10 +1748,6 @@ F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dtrtrs"); } } } @@ -1860,10 +1814,6 @@ F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in ztrcon"); - if (info != 0) info = -2; @@ -1899,10 +1849,6 @@ F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dtrtrs"); } } } @@ -1949,78 +1895,64 @@ tmp_data, nr, info F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zpotrf"); - else + // Throw-away extra info LAPACK gives so as to not change output. + rcond = 0.0; + if (info != 0) + { + info = -2; + + mattype.mark_as_unsymmetric (); + typ = MatrixType::Full; + } + else { - // Throw-away extra info LAPACK gives so as to not change output. - rcond = 0.0; - if (info != 0) + if (calc_cond) { - info = -2; - + Array<Complex> z (2 * nc); + Complex *pz = z.fortran_vec (); + Array<double> rz (nc); + double *prz = rz.fortran_vec (); + + F77_XFCN (zpocon, ZPOCON, (F77_CONST_CHAR_ARG2 (&job, 1), + nr, tmp_data, nr, anorm, + rcond, pz, prz, info + F77_CHAR_ARG_LEN (1))); + + if (info != 0) + info = -2; + + volatile double rcond_plus_one = rcond + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcond)) + { + info = -2; + + if (sing_handler) + sing_handler (rcond); + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision, rcond = %g", + rcond); + } + } + + if (info == 0) + { + retval = b; + Complex *result = retval.fortran_vec (); + + octave_idx_type b_nc = b.cols (); + + F77_XFCN (zpotrs, ZPOTRS, (F77_CONST_CHAR_ARG2 (&job, 1), + nr, b_nc, tmp_data, nr, + result, b.rows(), info + F77_CHAR_ARG_LEN (1))); + } + else + { mattype.mark_as_unsymmetric (); typ = MatrixType::Full; } - else - { - if (calc_cond) - { - Array<Complex> z (2 * nc); - Complex *pz = z.fortran_vec (); - Array<double> rz (nc); - double *prz = rz.fortran_vec (); - - F77_XFCN (zpocon, ZPOCON, (F77_CONST_CHAR_ARG2 (&job, 1), - nr, tmp_data, nr, anorm, - rcond, pz, prz, info - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zpocon"); - - if (info != 0) - info = -2; - - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) - { - info = -2; - - if (sing_handler) - sing_handler (rcond); - else - (*current_liboctave_error_handler) - ("matrix singular to machine precision, rcond = %g", - rcond); - } - } - - if (info == 0) - { - retval = b; - Complex *result = retval.fortran_vec (); - - octave_idx_type b_nc = b.cols (); - - F77_XFCN (zpotrs, ZPOTRS, (F77_CONST_CHAR_ARG2 (&job, 1), - nr, b_nc, tmp_data, nr, - result, b.rows(), info - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zpotrs"); - } - else - { - mattype.mark_as_unsymmetric (); - typ = MatrixType::Full; - } - } } } @@ -2045,79 +1977,65 @@ F77_XFCN (zgetrf, ZGETRF, (nr, nr, tmp_data, nr, pipvt, info)); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zgetrf"); - else + // Throw-away extra info LAPACK gives so as to not change output. + rcond = 0.0; + if (info != 0) + { + info = -2; + + if (sing_handler) + sing_handler (rcond); + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision"); + + mattype.mark_as_rectangular (); + } + else { - // Throw-away extra info LAPACK gives so as to not change output. - rcond = 0.0; - if (info != 0) - { - info = -2; - - if (sing_handler) - sing_handler (rcond); - else - (*current_liboctave_error_handler) - ("matrix singular to machine precision"); - - mattype.mark_as_rectangular (); - } - else + if (calc_cond) { - if (calc_cond) + // Now calculate the condition number for + // non-singular matrix. + char job = '1'; + F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1), + nc, tmp_data, nr, anorm, + rcond, pz, prz, info + F77_CHAR_ARG_LEN (1))); + + if (info != 0) + info = -2; + + volatile double rcond_plus_one = rcond + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcond)) { - // Now calculate the condition number for - // non-singular matrix. - char job = '1'; - F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1), - nc, tmp_data, nr, anorm, - rcond, pz, prz, info - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zgecon"); - - if (info != 0) - info = -2; - - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) - { - info = -2; - - if (sing_handler) - sing_handler (rcond); - else - (*current_liboctave_error_handler) - ("matrix singular to machine precision, rcond = %g", - rcond); - } + info = -2; + + if (sing_handler) + sing_handler (rcond); + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision, rcond = %g", + rcond); } - - if (info == 0) - { - retval = b; - Complex *result = retval.fortran_vec (); - - octave_idx_type b_nc = b.cols (); - - char job = 'N'; - F77_XFCN (zgetrs, ZGETRS, (F77_CONST_CHAR_ARG2 (&job, 1), - nr, b_nc, tmp_data, nr, - pipvt, result, b.rows(), info - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zgetrs"); - } - else - mattype.mark_as_rectangular (); } + + if (info == 0) + { + retval = b; + Complex *result = retval.fortran_vec (); + + octave_idx_type b_nc = b.cols (); + + char job = 'N'; + F77_XFCN (zgetrs, ZGETRS, (F77_CONST_CHAR_ARG2 (&job, 1), + nr, b_nc, tmp_data, nr, + pipvt, result, b.rows(), info + F77_CHAR_ARG_LEN (1))); + } + else + mattype.mark_as_rectangular (); } } } @@ -2561,37 +2479,25 @@ work(0) = lworkaround; } - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zgelsd"); + lwork = static_cast<octave_idx_type> (std::real (work(0))); + work.resize (lwork); + + F77_XFCN (zgelsd, ZGELSD, (m, n, nrhs, tmp_data, m, pretval, + maxmn, ps, rcond, rank, + work.fortran_vec (), lwork, + prwork, piwork, info)); + + if (rank < minmn) + (*current_liboctave_warning_handler) + ("zgelsd: rank deficient %dx%d matrix, rank = %d, tol = %e", + m, n, rank, rcond); + + if (s.elem (0) == 0.0) + rcond = 0.0; else - { - lwork = static_cast<octave_idx_type> (std::real (work(0))); - work.resize (lwork); - - F77_XFCN (zgelsd, ZGELSD, (m, n, nrhs, tmp_data, m, pretval, - maxmn, ps, rcond, rank, - work.fortran_vec (), lwork, - prwork, piwork, info)); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zgelsd"); - else - { - if (rank < minmn) - (*current_liboctave_warning_handler) - ("zgelsd: rank deficient %dx%d matrix, rank = %d, tol = %e", - m, n, rank, rcond); - - if (s.elem (0) == 0.0) - rcond = 0.0; - else - rcond = s.elem (minmn - 1) / s.elem (0); - - retval.resize (n, nrhs); - } - } + rcond = s.elem (minmn - 1) / s.elem (0); + + retval.resize (n, nrhs); } return retval; @@ -2733,38 +2639,29 @@ ps, rcond, rank, work.fortran_vec (), lwork, prwork, piwork, info)); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zgelsd"); - else + lwork = static_cast<octave_idx_type> (std::real (work(0))); + work.resize (lwork); + rwork.resize (static_cast<octave_idx_type> (rwork(0))); + iwork.resize (iwork(0)); + + F77_XFCN (zgelsd, ZGELSD, (m, n, nrhs, tmp_data, m, pretval, + maxmn, ps, rcond, rank, + work.fortran_vec (), lwork, + prwork, piwork, info)); + + if (rank < minmn) { - lwork = static_cast<octave_idx_type> (std::real (work(0))); - work.resize (lwork); - rwork.resize (static_cast<octave_idx_type> (rwork(0))); - iwork.resize (iwork(0)); - - F77_XFCN (zgelsd, ZGELSD, (m, n, nrhs, tmp_data, m, pretval, - maxmn, ps, rcond, rank, - work.fortran_vec (), lwork, - prwork, piwork, info)); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zgelsd"); - else if (rank < minmn) - { - if (rank < minmn) - (*current_liboctave_warning_handler) - ("zgelsd: rank deficient %dx%d matrix, rank = %d, tol = %e", - m, n, rank, rcond); - - if (s.elem (0) == 0.0) - rcond = 0.0; - else - rcond = s.elem (minmn - 1) / s.elem (0); - - retval.resize (n, nrhs); - } + if (rank < minmn) + (*current_liboctave_warning_handler) + ("zgelsd: rank deficient %dx%d matrix, rank = %d, tol = %e", + m, n, rank, rcond); + + if (s.elem (0) == 0.0) + rcond = 0.0; + else + rcond = s.elem (minmn - 1) / s.elem (0); + + retval.resize (n, nrhs); } } @@ -2841,12 +2738,6 @@ dpermute.fortran_vec (), info F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) ("unrecoverable error in zgebal"); - return retval; - } - // then scale job = 'S'; F77_XFCN (zgebal, ZGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1), @@ -2854,12 +2745,6 @@ dscale.fortran_vec (), info F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) ("unrecoverable error in zgebal"); - return retval; - } - // Preconditioning step 3: scaling. ColumnVector work (nc); @@ -2870,12 +2755,6 @@ work.fortran_vec (), inf_norm F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) ("unrecoverable error in zlange"); - return retval; - } - int sqpow = (inf_norm > 0.0 ? static_cast<int> (1.0 + log (inf_norm) / log (2.0)) : 0); @@ -3051,10 +2930,6 @@ a.data (), 1, 0.0, c, len F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zgemm"); } return retval; @@ -3939,14 +3814,9 @@ F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in ztrsyl"); - else - { - // FIXME -- check info? - - retval = -ua * cx * ub.hermitian (); - } + // FIXME -- check info? + + retval = -ua * cx * ub.hermitian (); return retval; } @@ -4016,10 +3886,6 @@ nr, nc, 1.0, m.data (), ld, a.data (), 1, 0.0, c, 1 F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zgemv"); } } else @@ -4030,10 +3896,6 @@ ld, a.data (), lda, 0.0, c, nr F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zgemm"); } } }
--- a/liboctave/CRowVector.cc Thu Feb 14 17:14:23 2008 -0500 +++ b/liboctave/CRowVector.cc Thu Feb 14 21:57:50 2008 -0500 @@ -351,10 +351,6 @@ a_nr, a_nc, 1.0, a.data (), ld, v.data (), 1, 0.0, y, 1 F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zgemv"); } }
--- a/liboctave/CSparse.cc Thu Feb 14 17:14:23 2008 -0500 +++ b/liboctave/CSparse.cc Thu Feb 14 21:57:50 2008 -0500 @@ -3749,10 +3749,7 @@ F77_XFCN (zptsv, ZPTSV, (nr, b_nc, D, DL, result, b.rows(), err)); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zptsv"); - else if (err != 0) + if (err != 0) { err = 0; mattype.mark_as_unsymmetric (); @@ -3809,10 +3806,7 @@ F77_XFCN (zgtsv, ZGTSV, (nr, b_nc, DL, D, DU, result, b.rows(), err)); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zgtsv"); - else if (err != 0) + if (err != 0) { rcond = 0.; err = -2; @@ -3910,84 +3904,71 @@ F77_XFCN (zgttrf, ZGTTRF, (nr, DL, D, DU, DU2, pipvt, err)); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zgttrf"); - else - { - if (err != 0) - { - err = -2; - rcond = 0.0; - - if (sing_handler) - { - sing_handler (rcond); - mattype.mark_as_rectangular (); - } - else - (*current_liboctave_error_handler) - ("matrix singular to machine precision"); - - } - else - { - char job = 'N'; - volatile octave_idx_type x_nz = b.nnz (); - octave_idx_type b_nc = b.cols (); - retval = SparseComplexMatrix (nr, b_nc, x_nz); - retval.xcidx(0) = 0; - volatile octave_idx_type ii = 0; - rcond = 1.0; - - OCTAVE_LOCAL_BUFFER (Complex, work, nr); - - for (volatile octave_idx_type j = 0; j < b_nc; j++) + if (err != 0) + { + err = -2; + rcond = 0.0; + + if (sing_handler) + { + sing_handler (rcond); + mattype.mark_as_rectangular (); + } + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision"); + + } + else + { + char job = 'N'; + volatile octave_idx_type x_nz = b.nnz (); + octave_idx_type b_nc = b.cols (); + retval = SparseComplexMatrix (nr, b_nc, x_nz); + retval.xcidx(0) = 0; + volatile octave_idx_type ii = 0; + rcond = 1.0; + + OCTAVE_LOCAL_BUFFER (Complex, work, nr); + + for (volatile octave_idx_type j = 0; j < b_nc; j++) + { + for (octave_idx_type i = 0; i < nr; i++) + work[i] = 0.; + for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++) + work[b.ridx(i)] = b.data(i); + + F77_XFCN (zgttrs, ZGTTRS, + (F77_CONST_CHAR_ARG2 (&job, 1), + nr, 1, DL, D, DU, DU2, pipvt, + work, b.rows (), err + F77_CHAR_ARG_LEN (1))); + + // Count non-zeros in work vector and adjust + // space in retval if needed + octave_idx_type new_nnz = 0; + for (octave_idx_type i = 0; i < nr; i++) + if (work[i] != 0.) + new_nnz++; + + if (ii + new_nnz > x_nz) { - for (octave_idx_type i = 0; i < nr; i++) - work[i] = 0.; - for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++) - work[b.ridx(i)] = b.data(i); - - F77_XFCN (zgttrs, ZGTTRS, - (F77_CONST_CHAR_ARG2 (&job, 1), - nr, 1, DL, D, DU, DU2, pipvt, - work, b.rows (), err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("unrecoverable error in zgttrs"); - break; - } - - // Count non-zeros in work vector and adjust - // space in retval if needed - octave_idx_type new_nnz = 0; - for (octave_idx_type i = 0; i < nr; i++) - if (work[i] != 0.) - new_nnz++; - - if (ii + new_nnz > x_nz) - { - // Resize the sparse matrix - octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; - retval.change_capacity (sz); - x_nz = sz; - } - - for (octave_idx_type i = 0; i < nr; i++) - if (work[i] != 0.) - { - retval.xridx(ii) = i; - retval.xdata(ii++) = work[i]; - } - retval.xcidx(j+1) = ii; + // Resize the sparse matrix + octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; + retval.change_capacity (sz); + x_nz = sz; } - retval.maybe_compress (); - } + for (octave_idx_type i = 0; i < nr; i++) + if (work[i] != 0.) + { + retval.xridx(ii) = i; + retval.xdata(ii++) = work[i]; + } + retval.xcidx(j+1) = ii; + } + + retval.maybe_compress (); } } else if (typ != MatrixType::Tridiagonal_Hermitian) @@ -4069,13 +4050,7 @@ F77_XFCN (zptsv, ZPTSV, (nr, b_nc, D, DL, result, b_nr, err)); - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("unrecoverable error in zptsv"); - err = -1; - } - else if (err != 0) + if (err != 0) { err = 0; mattype.mark_as_unsymmetric (); @@ -4133,13 +4108,7 @@ F77_XFCN (zgtsv, ZGTSV, (nr, b_nc, DL, D, DU, result, b_nr, err)); - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("unrecoverable error in zgtsv"); - err = -1; - } - else if (err != 0) + if (err != 0) { rcond = 0.; err = -2; @@ -4235,95 +4204,82 @@ F77_XFCN (zgttrf, ZGTTRF, (nr, DL, D, DU, DU2, pipvt, err)); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zgttrf"); - else - { - if (err != 0) - { - rcond = 0.0; - err = -2; - - if (sing_handler) - { - sing_handler (rcond); - mattype.mark_as_rectangular (); - } - else - (*current_liboctave_error_handler) - ("matrix singular to machine precision"); - } - else - { - rcond = 1.; - char job = 'N'; - octave_idx_type b_nr = b.rows (); - octave_idx_type b_nc = b.cols (); - OCTAVE_LOCAL_BUFFER (Complex, Bx, b_nr); - - // Take a first guess that the number of non-zero terms - // will be as many as in b - volatile octave_idx_type x_nz = b.nnz (); - volatile octave_idx_type ii = 0; - retval = SparseComplexMatrix (b_nr, b_nc, x_nz); - - retval.xcidx(0) = 0; - for (volatile octave_idx_type j = 0; j < b_nc; j++) + if (err != 0) + { + rcond = 0.0; + err = -2; + + if (sing_handler) + { + sing_handler (rcond); + mattype.mark_as_rectangular (); + } + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision"); + } + else + { + rcond = 1.; + char job = 'N'; + octave_idx_type b_nr = b.rows (); + octave_idx_type b_nc = b.cols (); + OCTAVE_LOCAL_BUFFER (Complex, Bx, b_nr); + + // Take a first guess that the number of non-zero terms + // will be as many as in b + volatile octave_idx_type x_nz = b.nnz (); + volatile octave_idx_type ii = 0; + retval = SparseComplexMatrix (b_nr, b_nc, x_nz); + + retval.xcidx(0) = 0; + for (volatile octave_idx_type j = 0; j < b_nc; j++) + { + + for (octave_idx_type i = 0; i < b_nr; i++) + Bx[i] = b (i,j); + + F77_XFCN (zgttrs, ZGTTRS, + (F77_CONST_CHAR_ARG2 (&job, 1), + nr, 1, DL, D, DU, DU2, pipvt, + Bx, b_nr, err + F77_CHAR_ARG_LEN (1))); + + if (err != 0) { - - for (octave_idx_type i = 0; i < b_nr; i++) - Bx[i] = b (i,j); - - F77_XFCN (zgttrs, ZGTTRS, - (F77_CONST_CHAR_ARG2 (&job, 1), - nr, 1, DL, D, DU, DU2, pipvt, - Bx, b_nr, err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("unrecoverable error in zgttrs"); - break; - } - - if (err != 0) - { - (*current_liboctave_error_handler) - ("SparseComplexMatrix::solve solve failed"); - - err = -1; - break; - } - - // Count non-zeros in work vector and adjust - // space in retval if needed - octave_idx_type new_nnz = 0; - for (octave_idx_type i = 0; i < nr; i++) - if (Bx[i] != 0.) - new_nnz++; - - if (ii + new_nnz > x_nz) - { - // Resize the sparse matrix - octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; - retval.change_capacity (sz); - x_nz = sz; - } - - for (octave_idx_type i = 0; i < nr; i++) - if (Bx[i] != 0.) - { - retval.xridx(ii) = i; - retval.xdata(ii++) = Bx[i]; - } - - retval.xcidx(j+1) = ii; + (*current_liboctave_error_handler) + ("SparseComplexMatrix::solve solve failed"); + + err = -1; + break; } - retval.maybe_compress (); - } + // Count non-zeros in work vector and adjust + // space in retval if needed + octave_idx_type new_nnz = 0; + for (octave_idx_type i = 0; i < nr; i++) + if (Bx[i] != 0.) + new_nnz++; + + if (ii + new_nnz > x_nz) + { + // Resize the sparse matrix + octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; + retval.change_capacity (sz); + x_nz = sz; + } + + for (octave_idx_type i = 0; i < nr; i++) + if (Bx[i] != 0.) + { + retval.xridx(ii) = i; + retval.xdata(ii++) = Bx[i]; + } + + retval.xcidx(j+1) = ii; + } + + retval.maybe_compress (); } } else if (typ != MatrixType::Tridiagonal_Hermitian) @@ -4390,85 +4346,71 @@ nr, n_lower, tmp_data, ldm, err F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zpbtrf"); - else - { - if (err != 0) - { - rcond = 0.0; - // Matrix is not positive definite!! Fall through to - // unsymmetric banded solver. - mattype.mark_as_unsymmetric (); - typ = MatrixType::Banded; - err = 0; - } - else - { - if (calc_cond) + if (err != 0) + { + rcond = 0.0; + // Matrix is not positive definite!! Fall through to + // unsymmetric banded solver. + mattype.mark_as_unsymmetric (); + typ = MatrixType::Banded; + err = 0; + } + else + { + if (calc_cond) + { + Array<Complex> z (2 * nr); + Complex *pz = z.fortran_vec (); + Array<double> iz (nr); + double *piz = iz.fortran_vec (); + + F77_XFCN (zpbcon, ZPBCON, + (F77_CONST_CHAR_ARG2 (&job, 1), + nr, n_lower, tmp_data, ldm, + anorm, rcond, pz, piz, err + F77_CHAR_ARG_LEN (1))); + + if (err != 0) + err = -2; + + volatile double rcond_plus_one = rcond + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcond)) { - Array<Complex> z (2 * nr); - Complex *pz = z.fortran_vec (); - Array<double> iz (nr); - double *piz = iz.fortran_vec (); - - F77_XFCN (zpbcon, ZPBCON, - (F77_CONST_CHAR_ARG2 (&job, 1), - nr, n_lower, tmp_data, ldm, - anorm, rcond, pz, piz, err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zpbcon"); - - if (err != 0) - err = -2; - - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) - { - err = -2; - - if (sing_handler) - { - sing_handler (rcond); - mattype.mark_as_rectangular (); - } - else - (*current_liboctave_error_handler) - ("matrix singular to machine precision, rcond = %g", - rcond); - } + err = -2; + + if (sing_handler) + { + sing_handler (rcond); + mattype.mark_as_rectangular (); + } + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision, rcond = %g", + rcond); } - else - rcond = 1.0; - - if (err == 0) + } + else + rcond = 1.0; + + if (err == 0) + { + retval = ComplexMatrix (b); + Complex *result = retval.fortran_vec (); + + octave_idx_type b_nc = b.cols (); + + F77_XFCN (zpbtrs, ZPBTRS, + (F77_CONST_CHAR_ARG2 (&job, 1), + nr, n_lower, b_nc, tmp_data, + ldm, result, b.rows(), err + F77_CHAR_ARG_LEN (1))); + + if (err != 0) { - retval = ComplexMatrix (b); - Complex *result = retval.fortran_vec (); - - octave_idx_type b_nc = b.cols (); - - F77_XFCN (zpbtrs, ZPBTRS, - (F77_CONST_CHAR_ARG2 (&job, 1), - nr, n_lower, b_nc, tmp_data, - ldm, result, b.rows(), err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zpbtrs"); - - if (err != 0) - { - (*current_liboctave_error_handler) - ("SparseMatrix::solve solve failed"); - err = -1; - } + (*current_liboctave_error_handler) + ("SparseMatrix::solve solve failed"); + err = -1; } } } @@ -4517,88 +4459,74 @@ F77_XFCN (zgbtrf, ZGBTRF, (nr, nc, n_lower, n_upper, tmp_data, ldm, pipvt, err)); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zgbtrf"); - else - { - // Throw-away extra info LAPACK gives so as to not - // change output. - if (err != 0) - { - rcond = 0.0; - err = -2; - - if (sing_handler) - { - sing_handler (rcond); - mattype.mark_as_rectangular (); - } - else - (*current_liboctave_error_handler) - ("matrix singular to machine precision"); - } - else - { - if (calc_cond) + // Throw-away extra info LAPACK gives so as to not + // change output. + if (err != 0) + { + rcond = 0.0; + err = -2; + + if (sing_handler) + { + sing_handler (rcond); + mattype.mark_as_rectangular (); + } + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision"); + } + else + { + if (calc_cond) + { + char job = '1'; + Array<Complex> z (2 * nr); + Complex *pz = z.fortran_vec (); + Array<double> iz (nr); + double *piz = iz.fortran_vec (); + + F77_XFCN (zgbcon, ZGBCON, + (F77_CONST_CHAR_ARG2 (&job, 1), + nc, n_lower, n_upper, tmp_data, ldm, pipvt, + anorm, rcond, pz, piz, err + F77_CHAR_ARG_LEN (1))); + + if (err != 0) + err = -2; + + volatile double rcond_plus_one = rcond + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcond)) { - char job = '1'; - Array<Complex> z (2 * nr); - Complex *pz = z.fortran_vec (); - Array<double> iz (nr); - double *piz = iz.fortran_vec (); - - F77_XFCN (zgbcon, ZGBCON, - (F77_CONST_CHAR_ARG2 (&job, 1), - nc, n_lower, n_upper, tmp_data, ldm, pipvt, - anorm, rcond, pz, piz, err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zgbcon"); - - if (err != 0) - err = -2; - - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) - { - err = -2; - - if (sing_handler) - { - sing_handler (rcond); - mattype.mark_as_rectangular (); - } - else - (*current_liboctave_error_handler) - ("matrix singular to machine precision, rcond = %g", - rcond); - } + err = -2; + + if (sing_handler) + { + sing_handler (rcond); + mattype.mark_as_rectangular (); + } + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision, rcond = %g", + rcond); } - else - rcond = 1.; - - if (err == 0) - { - retval = ComplexMatrix (b); - Complex *result = retval.fortran_vec (); - - octave_idx_type b_nc = b.cols (); - - char job = 'N'; - F77_XFCN (zgbtrs, ZGBTRS, - (F77_CONST_CHAR_ARG2 (&job, 1), - nr, n_lower, n_upper, b_nc, tmp_data, - ldm, pipvt, result, b.rows(), err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zgbtrs"); - } + } + else + rcond = 1.; + + if (err == 0) + { + retval = ComplexMatrix (b); + Complex *result = retval.fortran_vec (); + + octave_idx_type b_nc = b.cols (); + + char job = 'N'; + F77_XFCN (zgbtrs, ZGBTRS, + (F77_CONST_CHAR_ARG2 (&job, 1), + nr, n_lower, n_upper, b_nc, tmp_data, + ldm, pipvt, result, b.rows(), err + F77_CHAR_ARG_LEN (1))); } } } @@ -4667,123 +4595,105 @@ nr, n_lower, tmp_data, ldm, err F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zpbtrf"); - else - { - if (err != 0) - { - rcond = 0.0; - mattype.mark_as_unsymmetric (); - typ = MatrixType::Banded; - err = 0; - } - else - { - if (calc_cond) - { - Array<Complex> z (2 * nr); - Complex *pz = z.fortran_vec (); - Array<double> iz (nr); - double *piz = iz.fortran_vec (); - - F77_XFCN (zpbcon, ZPBCON, - (F77_CONST_CHAR_ARG2 (&job, 1), - nr, n_lower, tmp_data, ldm, - anorm, rcond, pz, piz, err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zpbcon"); - - if (err != 0) - err = -2; - - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) - { - err = -2; - - if (sing_handler) - { - sing_handler (rcond); - mattype.mark_as_rectangular (); - } - else - (*current_liboctave_error_handler) - ("matrix singular to machine precision, rcond = %g", - rcond); - } - } - else - rcond = 1.0; - - if (err == 0) + if (err != 0) + { + rcond = 0.0; + mattype.mark_as_unsymmetric (); + typ = MatrixType::Banded; + err = 0; + } + else + { + if (calc_cond) + { + Array<Complex> z (2 * nr); + Complex *pz = z.fortran_vec (); + Array<double> iz (nr); + double *piz = iz.fortran_vec (); + + F77_XFCN (zpbcon, ZPBCON, + (F77_CONST_CHAR_ARG2 (&job, 1), + nr, n_lower, tmp_data, ldm, + anorm, rcond, pz, piz, err + F77_CHAR_ARG_LEN (1))); + + if (err != 0) + err = -2; + + volatile double rcond_plus_one = rcond + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcond)) { - octave_idx_type b_nr = b.rows (); - octave_idx_type b_nc = b.cols (); - OCTAVE_LOCAL_BUFFER (Complex, Bx, b_nr); - - // Take a first guess that the number of non-zero terms - // will be as many as in b - volatile octave_idx_type x_nz = b.nnz (); - volatile octave_idx_type ii = 0; - retval = SparseComplexMatrix (b_nr, b_nc, x_nz); - - retval.xcidx(0) = 0; - for (volatile octave_idx_type j = 0; j < b_nc; j++) + err = -2; + + if (sing_handler) { - for (octave_idx_type i = 0; i < b_nr; i++) - Bx[i] = b.elem (i, j); - - F77_XFCN (zpbtrs, ZPBTRS, - (F77_CONST_CHAR_ARG2 (&job, 1), - nr, n_lower, 1, tmp_data, - ldm, Bx, b_nr, err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) + sing_handler (rcond); + mattype.mark_as_rectangular (); + } + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision, rcond = %g", + rcond); + } + } + else + rcond = 1.0; + + if (err == 0) + { + octave_idx_type b_nr = b.rows (); + octave_idx_type b_nc = b.cols (); + OCTAVE_LOCAL_BUFFER (Complex, Bx, b_nr); + + // Take a first guess that the number of non-zero terms + // will be as many as in b + volatile octave_idx_type x_nz = b.nnz (); + volatile octave_idx_type ii = 0; + retval = SparseComplexMatrix (b_nr, b_nc, x_nz); + + retval.xcidx(0) = 0; + for (volatile octave_idx_type j = 0; j < b_nc; j++) + { + for (octave_idx_type i = 0; i < b_nr; i++) + Bx[i] = b.elem (i, j); + + F77_XFCN (zpbtrs, ZPBTRS, + (F77_CONST_CHAR_ARG2 (&job, 1), + nr, n_lower, 1, tmp_data, + ldm, Bx, b_nr, err + F77_CHAR_ARG_LEN (1))); + + if (err != 0) + { + (*current_liboctave_error_handler) + ("SparseComplexMatrix::solve solve failed"); + err = -1; + break; + } + + for (octave_idx_type i = 0; i < b_nr; i++) + { + Complex tmp = Bx[i]; + if (tmp != 0.0) { - (*current_liboctave_error_handler) - ("unrecoverable error in dpbtrs"); - err = -1; - break; + if (ii == x_nz) + { + // Resize the sparse matrix + octave_idx_type sz = x_nz * + (b_nc - j) / b_nc; + sz = (sz > 10 ? sz : 10) + x_nz; + retval.change_capacity (sz); + x_nz = sz; + } + retval.xdata(ii) = tmp; + retval.xridx(ii++) = i; } - - if (err != 0) - { - (*current_liboctave_error_handler) - ("SparseComplexMatrix::solve solve failed"); - err = -1; - break; - } - - for (octave_idx_type i = 0; i < b_nr; i++) - { - Complex tmp = Bx[i]; - if (tmp != 0.0) - { - if (ii == x_nz) - { - // Resize the sparse matrix - octave_idx_type sz = x_nz * - (b_nc - j) / b_nc; - sz = (sz > 10 ? sz : 10) + x_nz; - retval.change_capacity (sz); - x_nz = sz; - } - retval.xdata(ii) = tmp; - retval.xridx(ii++) = i; - } - } - retval.xcidx(j+1) = ii; } - - retval.maybe_compress (); + retval.xcidx(j+1) = ii; } + + retval.maybe_compress (); } } } @@ -4831,127 +4741,110 @@ F77_XFCN (zgbtrf, ZGBTRF, (nr, nr, n_lower, n_upper, tmp_data, ldm, pipvt, err)); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zgbtrf"); - else - { - if (err != 0) - { - rcond = 0.0; - err = -2; - - if (sing_handler) - { - sing_handler (rcond); - mattype.mark_as_rectangular (); - } - else - (*current_liboctave_error_handler) - ("matrix singular to machine precision"); - - } - else - { - if (calc_cond) + if (err != 0) + { + rcond = 0.0; + err = -2; + + if (sing_handler) + { + sing_handler (rcond); + mattype.mark_as_rectangular (); + } + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision"); + + } + else + { + if (calc_cond) + { + char job = '1'; + Array<Complex> z (2 * nr); + Complex *pz = z.fortran_vec (); + Array<double> iz (nr); + double *piz = iz.fortran_vec (); + + F77_XFCN (zgbcon, ZGBCON, + (F77_CONST_CHAR_ARG2 (&job, 1), + nc, n_lower, n_upper, tmp_data, ldm, pipvt, + anorm, rcond, pz, piz, err + F77_CHAR_ARG_LEN (1))); + + if (err != 0) + err = -2; + + volatile double rcond_plus_one = rcond + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcond)) { - char job = '1'; - Array<Complex> z (2 * nr); - Complex *pz = z.fortran_vec (); - Array<double> iz (nr); - double *piz = iz.fortran_vec (); - - F77_XFCN (zgbcon, ZGBCON, - (F77_CONST_CHAR_ARG2 (&job, 1), - nc, n_lower, n_upper, tmp_data, ldm, pipvt, - anorm, rcond, pz, piz, err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zgbcon"); - - if (err != 0) - err = -2; - - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) - { - err = -2; - - if (sing_handler) - { - sing_handler (rcond); - mattype.mark_as_rectangular (); - } - else - (*current_liboctave_error_handler) - ("matrix singular to machine precision, rcond = %g", - rcond); - } + err = -2; + + if (sing_handler) + { + sing_handler (rcond); + mattype.mark_as_rectangular (); + } + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision, rcond = %g", + rcond); } - else - rcond = 1.; - - if (err == 0) + } + else + rcond = 1.; + + if (err == 0) + { + char job = 'N'; + volatile octave_idx_type x_nz = b.nnz (); + octave_idx_type b_nc = b.cols (); + retval = SparseComplexMatrix (nr, b_nc, x_nz); + retval.xcidx(0) = 0; + volatile octave_idx_type ii = 0; + + OCTAVE_LOCAL_BUFFER (Complex, work, nr); + + for (volatile octave_idx_type j = 0; j < b_nc; j++) { - char job = 'N'; - volatile octave_idx_type x_nz = b.nnz (); - octave_idx_type b_nc = b.cols (); - retval = SparseComplexMatrix (nr, b_nc, x_nz); - retval.xcidx(0) = 0; - volatile octave_idx_type ii = 0; - - OCTAVE_LOCAL_BUFFER (Complex, work, nr); - - for (volatile octave_idx_type j = 0; j < b_nc; j++) + for (octave_idx_type i = 0; i < nr; i++) + work[i] = 0.; + for (octave_idx_type i = b.cidx(j); + i < b.cidx(j+1); i++) + work[b.ridx(i)] = b.data(i); + + F77_XFCN (zgbtrs, ZGBTRS, + (F77_CONST_CHAR_ARG2 (&job, 1), + nr, n_lower, n_upper, 1, tmp_data, + ldm, pipvt, work, b.rows (), err + F77_CHAR_ARG_LEN (1))); + + // Count non-zeros in work vector and adjust + // space in retval if needed + octave_idx_type new_nnz = 0; + for (octave_idx_type i = 0; i < nr; i++) + if (work[i] != 0.) + new_nnz++; + + if (ii + new_nnz > x_nz) { - for (octave_idx_type i = 0; i < nr; i++) - work[i] = 0.; - for (octave_idx_type i = b.cidx(j); - i < b.cidx(j+1); i++) - work[b.ridx(i)] = b.data(i); - - F77_XFCN (zgbtrs, ZGBTRS, - (F77_CONST_CHAR_ARG2 (&job, 1), - nr, n_lower, n_upper, 1, tmp_data, - ldm, pipvt, work, b.rows (), err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("unrecoverable error in zgbtrs"); - break; - } - - // Count non-zeros in work vector and adjust - // space in retval if needed - octave_idx_type new_nnz = 0; - for (octave_idx_type i = 0; i < nr; i++) - if (work[i] != 0.) - new_nnz++; - - if (ii + new_nnz > x_nz) - { - // Resize the sparse matrix - octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; - retval.change_capacity (sz); - x_nz = sz; - } - - for (octave_idx_type i = 0; i < nr; i++) - if (work[i] != 0.) - { - retval.xridx(ii) = i; - retval.xdata(ii++) = work[i]; - } - retval.xcidx(j+1) = ii; + // Resize the sparse matrix + octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; + retval.change_capacity (sz); + x_nz = sz; } - retval.maybe_compress (); + for (octave_idx_type i = 0; i < nr; i++) + if (work[i] != 0.) + { + retval.xridx(ii) = i; + retval.xdata(ii++) = work[i]; + } + retval.xcidx(j+1) = ii; } + + retval.maybe_compress (); } } } @@ -5020,88 +4913,71 @@ nr, n_lower, tmp_data, ldm, err F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zpbtrf"); - else - { - if (err != 0) - { - // Matrix is not positive definite!! Fall through to - // unsymmetric banded solver. - rcond = 0.0; - mattype.mark_as_unsymmetric (); - typ = MatrixType::Banded; - err = 0; - } - else - { - if (calc_cond) + if (err != 0) + { + // Matrix is not positive definite!! Fall through to + // unsymmetric banded solver. + rcond = 0.0; + mattype.mark_as_unsymmetric (); + typ = MatrixType::Banded; + err = 0; + } + else + { + if (calc_cond) + { + Array<Complex> z (2 * nr); + Complex *pz = z.fortran_vec (); + Array<double> iz (nr); + double *piz = iz.fortran_vec (); + + F77_XFCN (zpbcon, ZPBCON, + (F77_CONST_CHAR_ARG2 (&job, 1), + nr, n_lower, tmp_data, ldm, + anorm, rcond, pz, piz, err + F77_CHAR_ARG_LEN (1))); + + if (err != 0) + err = -2; + + volatile double rcond_plus_one = rcond + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcond)) { - Array<Complex> z (2 * nr); - Complex *pz = z.fortran_vec (); - Array<double> iz (nr); - double *piz = iz.fortran_vec (); - - F77_XFCN (zpbcon, ZPBCON, - (F77_CONST_CHAR_ARG2 (&job, 1), - nr, n_lower, tmp_data, ldm, - anorm, rcond, pz, piz, err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zpbcon"); - - if (err != 0) - err = -2; - - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) - { - err = -2; - - if (sing_handler) - { - sing_handler (rcond); - mattype.mark_as_rectangular (); - } - else - (*current_liboctave_error_handler) - ("matrix singular to machine precision, rcond = %g", - rcond); - } + err = -2; + + if (sing_handler) + { + sing_handler (rcond); + mattype.mark_as_rectangular (); + } + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision, rcond = %g", + rcond); } - else - rcond = 1.0; - - if (err == 0) + } + else + rcond = 1.0; + + if (err == 0) + { + octave_idx_type b_nr = b.rows (); + octave_idx_type b_nc = b.cols (); + retval = ComplexMatrix (b); + Complex *result = retval.fortran_vec (); + + F77_XFCN (zpbtrs, ZPBTRS, + (F77_CONST_CHAR_ARG2 (&job, 1), + nr, n_lower, b_nc, tmp_data, + ldm, result, b_nr, err + F77_CHAR_ARG_LEN (1))); + + if (err != 0) { - octave_idx_type b_nr = b.rows (); - octave_idx_type b_nc = b.cols (); - retval = ComplexMatrix (b); - Complex *result = retval.fortran_vec (); - - F77_XFCN (zpbtrs, ZPBTRS, - (F77_CONST_CHAR_ARG2 (&job, 1), - nr, n_lower, b_nc, tmp_data, - ldm, result, b_nr, err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("unrecoverable error in zpbtrs"); - err = -1; - } - - if (err != 0) - { - (*current_liboctave_error_handler) - ("SparseComplexMatrix::solve solve failed"); - err = -1; - } + (*current_liboctave_error_handler) + ("SparseComplexMatrix::solve solve failed"); + err = -1; } } } @@ -5150,87 +5026,71 @@ F77_XFCN (zgbtrf, ZGBTRF, (nr, nr, n_lower, n_upper, tmp_data, ldm, pipvt, err)); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zgbtrf"); - else - { - if (err != 0) - { - err = -2; - rcond = 0.0; - - if (sing_handler) - { - sing_handler (rcond); - mattype.mark_as_rectangular (); - } - else - (*current_liboctave_error_handler) - ("matrix singular to machine precision"); - } - else - { - if (calc_cond) + if (err != 0) + { + err = -2; + rcond = 0.0; + + if (sing_handler) + { + sing_handler (rcond); + mattype.mark_as_rectangular (); + } + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision"); + } + else + { + if (calc_cond) + { + char job = '1'; + Array<Complex> z (2 * nr); + Complex *pz = z.fortran_vec (); + Array<double> iz (nr); + double *piz = iz.fortran_vec (); + + F77_XFCN (zgbcon, ZGBCON, + (F77_CONST_CHAR_ARG2 (&job, 1), + nc, n_lower, n_upper, tmp_data, ldm, pipvt, + anorm, rcond, pz, piz, err + F77_CHAR_ARG_LEN (1))); + + if (err != 0) + err = -2; + + volatile double rcond_plus_one = rcond + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcond)) { - char job = '1'; - Array<Complex> z (2 * nr); - Complex *pz = z.fortran_vec (); - Array<double> iz (nr); - double *piz = iz.fortran_vec (); - - F77_XFCN (zgbcon, ZGBCON, - (F77_CONST_CHAR_ARG2 (&job, 1), - nc, n_lower, n_upper, tmp_data, ldm, pipvt, - anorm, rcond, pz, piz, err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zgbcon"); - - if (err != 0) - err = -2; - - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) - { - err = -2; - - if (sing_handler) - { - sing_handler (rcond); - mattype.mark_as_rectangular (); - } - else - (*current_liboctave_error_handler) - ("matrix singular to machine precision, rcond = %g", - rcond); - } + err = -2; + + if (sing_handler) + { + sing_handler (rcond); + mattype.mark_as_rectangular (); + } + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision, rcond = %g", + rcond); } - else - rcond = 1.; - - if (err == 0) - { - char job = 'N'; - octave_idx_type b_nc = b.cols (); - retval = ComplexMatrix (b); - Complex *result = retval.fortran_vec (); - - F77_XFCN (zgbtrs, ZGBTRS, - (F77_CONST_CHAR_ARG2 (&job, 1), - nr, n_lower, n_upper, b_nc, tmp_data, - ldm, pipvt, result, b.rows (), err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("unrecoverable error in dgbtrs"); - } - } + } + else + rcond = 1.; + + if (err == 0) + { + char job = 'N'; + octave_idx_type b_nc = b.cols (); + retval = ComplexMatrix (b); + Complex *result = retval.fortran_vec (); + + F77_XFCN (zgbtrs, ZGBTRS, + (F77_CONST_CHAR_ARG2 (&job, 1), + nr, n_lower, n_upper, b_nc, tmp_data, + ldm, pipvt, result, b.rows (), err + F77_CHAR_ARG_LEN (1))); } } } @@ -5299,132 +5159,114 @@ nr, n_lower, tmp_data, ldm, err F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zpbtrf"); - else - { - if (err != 0) - { - // Matrix is not positive definite!! Fall through to - // unsymmetric banded solver. - mattype.mark_as_unsymmetric (); - typ = MatrixType::Banded; - - rcond = 0.0; - err = 0; - } - else - { - if (calc_cond) - { - Array<Complex> z (2 * nr); - Complex *pz = z.fortran_vec (); - Array<double> iz (nr); - double *piz = iz.fortran_vec (); - - F77_XFCN (zpbcon, ZPBCON, - (F77_CONST_CHAR_ARG2 (&job, 1), - nr, n_lower, tmp_data, ldm, - anorm, rcond, pz, piz, err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zpbcon"); - - if (err != 0) - err = -2; - - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) - { - err = -2; - - if (sing_handler) - { - sing_handler (rcond); - mattype.mark_as_rectangular (); - } - else - (*current_liboctave_error_handler) - ("matrix singular to machine precision, rcond = %g", - rcond); - } - } - else - rcond = 1.0; - - if (err == 0) + if (err != 0) + { + // Matrix is not positive definite!! Fall through to + // unsymmetric banded solver. + mattype.mark_as_unsymmetric (); + typ = MatrixType::Banded; + + rcond = 0.0; + err = 0; + } + else + { + if (calc_cond) + { + Array<Complex> z (2 * nr); + Complex *pz = z.fortran_vec (); + Array<double> iz (nr); + double *piz = iz.fortran_vec (); + + F77_XFCN (zpbcon, ZPBCON, + (F77_CONST_CHAR_ARG2 (&job, 1), + nr, n_lower, tmp_data, ldm, + anorm, rcond, pz, piz, err + F77_CHAR_ARG_LEN (1))); + + if (err != 0) + err = -2; + + volatile double rcond_plus_one = rcond + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcond)) { - octave_idx_type b_nr = b.rows (); - octave_idx_type b_nc = b.cols (); - OCTAVE_LOCAL_BUFFER (Complex, Bx, b_nr); - - // Take a first guess that the number of non-zero terms - // will be as many as in b - volatile octave_idx_type x_nz = b.nnz (); - volatile octave_idx_type ii = 0; - retval = SparseComplexMatrix (b_nr, b_nc, x_nz); - - retval.xcidx(0) = 0; - for (volatile octave_idx_type j = 0; j < b_nc; j++) + err = -2; + + if (sing_handler) { - - for (octave_idx_type i = 0; i < b_nr; i++) - Bx[i] = b (i,j); - - F77_XFCN (zpbtrs, ZPBTRS, - (F77_CONST_CHAR_ARG2 (&job, 1), - nr, n_lower, 1, tmp_data, - ldm, Bx, b_nr, err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("unrecoverable error in zpbtrs"); - err = -1; - break; - } - - if (err != 0) - { - (*current_liboctave_error_handler) - ("SparseMatrix::solve solve failed"); - err = -1; - break; - } - - // Count non-zeros in work vector and adjust - // space in retval if needed - octave_idx_type new_nnz = 0; - for (octave_idx_type i = 0; i < nr; i++) - if (Bx[i] != 0.) - new_nnz++; - - if (ii + new_nnz > x_nz) - { - // Resize the sparse matrix - octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; - retval.change_capacity (sz); - x_nz = sz; - } - - for (octave_idx_type i = 0; i < nr; i++) - if (Bx[i] != 0.) - { - retval.xridx(ii) = i; - retval.xdata(ii++) = Bx[i]; - } - - retval.xcidx(j+1) = ii; + sing_handler (rcond); + mattype.mark_as_rectangular (); } - - retval.maybe_compress (); + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision, rcond = %g", + rcond); } } + else + rcond = 1.0; + + if (err == 0) + { + octave_idx_type b_nr = b.rows (); + octave_idx_type b_nc = b.cols (); + OCTAVE_LOCAL_BUFFER (Complex, Bx, b_nr); + + // Take a first guess that the number of non-zero terms + // will be as many as in b + volatile octave_idx_type x_nz = b.nnz (); + volatile octave_idx_type ii = 0; + retval = SparseComplexMatrix (b_nr, b_nc, x_nz); + + retval.xcidx(0) = 0; + for (volatile octave_idx_type j = 0; j < b_nc; j++) + { + + for (octave_idx_type i = 0; i < b_nr; i++) + Bx[i] = b (i,j); + + F77_XFCN (zpbtrs, ZPBTRS, + (F77_CONST_CHAR_ARG2 (&job, 1), + nr, n_lower, 1, tmp_data, + ldm, Bx, b_nr, err + F77_CHAR_ARG_LEN (1))); + + if (err != 0) + { + (*current_liboctave_error_handler) + ("SparseMatrix::solve solve failed"); + err = -1; + break; + } + + // Count non-zeros in work vector and adjust + // space in retval if needed + octave_idx_type new_nnz = 0; + for (octave_idx_type i = 0; i < nr; i++) + if (Bx[i] != 0.) + new_nnz++; + + if (ii + new_nnz > x_nz) + { + // Resize the sparse matrix + octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; + retval.change_capacity (sz); + x_nz = sz; + } + + for (octave_idx_type i = 0; i < nr; i++) + if (Bx[i] != 0.) + { + retval.xridx(ii) = i; + retval.xdata(ii++) = Bx[i]; + } + + retval.xcidx(j+1) = ii; + } + + retval.maybe_compress (); + } } } @@ -5471,128 +5313,111 @@ F77_XFCN (zgbtrf, ZGBTRF, (nr, nr, n_lower, n_upper, tmp_data, ldm, pipvt, err)); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in xgbtrf"); - else - { - if (err != 0) - { - err = -2; - rcond = 0.0; - - if (sing_handler) - { - sing_handler (rcond); - mattype.mark_as_rectangular (); - } - else - (*current_liboctave_error_handler) - ("matrix singular to machine precision"); - - } - else - { - if (calc_cond) + if (err != 0) + { + err = -2; + rcond = 0.0; + + if (sing_handler) + { + sing_handler (rcond); + mattype.mark_as_rectangular (); + } + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision"); + + } + else + { + if (calc_cond) + { + char job = '1'; + Array<Complex> z (2 * nr); + Complex *pz = z.fortran_vec (); + Array<double> iz (nr); + double *piz = iz.fortran_vec (); + + F77_XFCN (zgbcon, ZGBCON, + (F77_CONST_CHAR_ARG2 (&job, 1), + nc, n_lower, n_upper, tmp_data, ldm, pipvt, + anorm, rcond, pz, piz, err + F77_CHAR_ARG_LEN (1))); + + if (err != 0) + err = -2; + + volatile double rcond_plus_one = rcond + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcond)) { - char job = '1'; - Array<Complex> z (2 * nr); - Complex *pz = z.fortran_vec (); - Array<double> iz (nr); - double *piz = iz.fortran_vec (); - - F77_XFCN (zgbcon, ZGBCON, - (F77_CONST_CHAR_ARG2 (&job, 1), - nc, n_lower, n_upper, tmp_data, ldm, pipvt, - anorm, rcond, pz, piz, err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zgbcon"); - - if (err != 0) - err = -2; - - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) - { - err = -2; - - if (sing_handler) - { - sing_handler (rcond); - mattype.mark_as_rectangular (); - } - else - (*current_liboctave_error_handler) - ("matrix singular to machine precision, rcond = %g", - rcond); - } + err = -2; + + if (sing_handler) + { + sing_handler (rcond); + mattype.mark_as_rectangular (); + } + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision, rcond = %g", + rcond); } - else - rcond = 1.; - - if (err == 0) + } + else + rcond = 1.; + + if (err == 0) + { + char job = 'N'; + volatile octave_idx_type x_nz = b.nnz (); + octave_idx_type b_nc = b.cols (); + retval = SparseComplexMatrix (nr, b_nc, x_nz); + retval.xcidx(0) = 0; + volatile octave_idx_type ii = 0; + + OCTAVE_LOCAL_BUFFER (Complex, Bx, nr); + + for (volatile octave_idx_type j = 0; j < b_nc; j++) { - char job = 'N'; - volatile octave_idx_type x_nz = b.nnz (); - octave_idx_type b_nc = b.cols (); - retval = SparseComplexMatrix (nr, b_nc, x_nz); - retval.xcidx(0) = 0; - volatile octave_idx_type ii = 0; - - OCTAVE_LOCAL_BUFFER (Complex, Bx, nr); - - for (volatile octave_idx_type j = 0; j < b_nc; j++) + for (octave_idx_type i = 0; i < nr; i++) + Bx[i] = 0.; + + for (octave_idx_type i = b.cidx(j); + i < b.cidx(j+1); i++) + Bx[b.ridx(i)] = b.data(i); + + F77_XFCN (zgbtrs, ZGBTRS, + (F77_CONST_CHAR_ARG2 (&job, 1), + nr, n_lower, n_upper, 1, tmp_data, + ldm, pipvt, Bx, b.rows (), err + F77_CHAR_ARG_LEN (1))); + + // Count non-zeros in work vector and adjust + // space in retval if needed + octave_idx_type new_nnz = 0; + for (octave_idx_type i = 0; i < nr; i++) + if (Bx[i] != 0.) + new_nnz++; + + if (ii + new_nnz > x_nz) { - for (octave_idx_type i = 0; i < nr; i++) - Bx[i] = 0.; - - for (octave_idx_type i = b.cidx(j); - i < b.cidx(j+1); i++) - Bx[b.ridx(i)] = b.data(i); - - F77_XFCN (zgbtrs, ZGBTRS, - (F77_CONST_CHAR_ARG2 (&job, 1), - nr, n_lower, n_upper, 1, tmp_data, - ldm, pipvt, Bx, b.rows (), err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("unrecoverable error in dgbtrs"); - break; - } - - // Count non-zeros in work vector and adjust - // space in retval if needed - octave_idx_type new_nnz = 0; - for (octave_idx_type i = 0; i < nr; i++) - if (Bx[i] != 0.) - new_nnz++; - - if (ii + new_nnz > x_nz) - { - // Resize the sparse matrix - octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; - retval.change_capacity (sz); - x_nz = sz; - } - - for (octave_idx_type i = 0; i < nr; i++) - if (Bx[i] != 0.) - { - retval.xridx(ii) = i; - retval.xdata(ii++) = Bx[i]; - } - retval.xcidx(j+1) = ii; + // Resize the sparse matrix + octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; + retval.change_capacity (sz); + x_nz = sz; } - retval.maybe_compress (); + for (octave_idx_type i = 0; i < nr; i++) + if (Bx[i] != 0.) + { + retval.xridx(ii) = i; + retval.xdata(ii++) = Bx[i]; + } + retval.xcidx(j+1) = ii; } + + retval.maybe_compress (); } } }
--- a/liboctave/ChangeLog Thu Feb 14 17:14:23 2008 -0500 +++ b/liboctave/ChangeLog Thu Feb 14 21:57:50 2008 -0500 @@ -1,3 +1,14 @@ +2008-02-14 John W. Eaton <jwe@octave.org> + + * CColVector.cc, CMatrix.cc, CRowVector.cc, CSparse.cc, + CmplxAEPBAL.cc, CmplxCHOL.cc, CmplxHESS.cc, CmplxLU.cc, + CmplxQR.cc, CmplxQRP.cc, CmplxSCHUR.cc, CmplxSVD.cc, DASPK.cc, + DASRT.cc, DASSL.cc, EIG.cc, LSODE.cc, NLEqn.cc, Quad.cc, + dColVector.cc, dMatrix.cc, dRowVector.cc, dSparse.cc, + dbleAEPBAL.cc, dbleCHOL.cc, dbleHESS.cc, dbleLU.cc, dbleQR.cc, + dbleQRP.cc, dbleSCHUR.cc, dbleSVD.cc: + Don't check f77_exception_encountered. + 2008-02-12 John W. Eaton <jwe@octave.org> * CMatrix.cc: Declare xilaenv instead of ilaenv.
--- a/liboctave/CmplxAEPBAL.cc Thu Feb 14 17:14:23 2008 -0500 +++ b/liboctave/CmplxAEPBAL.cc Thu Feb 14 21:57:50 2008 -0500 @@ -77,28 +77,20 @@ pscale, info F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in zgebal"); - else - { - balancing_mat = ComplexMatrix (n, n, 0.0); - for (octave_idx_type i = 0; i < n; i++) - balancing_mat.elem (i, i) = 1.0; + balancing_mat = ComplexMatrix (n, n, 0.0); + for (octave_idx_type i = 0; i < n; i++) + balancing_mat.elem (i, i) = 1.0; - Complex *p_balancing_mat = balancing_mat.fortran_vec (); + Complex *p_balancing_mat = balancing_mat.fortran_vec (); - char side = 'R'; + char side = 'R'; - F77_XFCN (zgebak, ZGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1), - F77_CONST_CHAR_ARG2 (&side, 1), - n, ilo, ihi, pscale, n, - p_balancing_mat, n, info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in zgebak"); - } + F77_XFCN (zgebak, ZGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1), + F77_CONST_CHAR_ARG2 (&side, 1), + n, ilo, ihi, pscale, n, + p_balancing_mat, n, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); return info; }
--- a/liboctave/CmplxCHOL.cc Thu Feb 14 17:14:23 2008 -0500 +++ b/liboctave/CmplxCHOL.cc Thu Feb 14 21:57:50 2008 -0500 @@ -76,43 +76,34 @@ F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, n, info F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in zpotrf"); + xrcond = 0.0; + if (info != 0) + info = -1; + else if (calc_cond) + { + octave_idx_type zpocon_info = 0; + + // Now calculate the condition number for non-singular matrix. + Array<Complex> z (2*n); + Complex *pz = z.fortran_vec (); + Array<double> rz (n); + double *prz = rz.fortran_vec (); + F77_XFCN (zpocon, ZPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, + n, anorm, xrcond, pz, prz, zpocon_info + F77_CHAR_ARG_LEN (1))); + + if (zpocon_info != 0) + info = -1; + } else { - xrcond = 0.0; - if (info != 0) - info = -1; - else if (calc_cond) - { - octave_idx_type zpocon_info = 0; - - // Now calculate the condition number for non-singular matrix. - Array<Complex> z (2*n); - Complex *pz = z.fortran_vec (); - Array<double> rz (n); - double *prz = rz.fortran_vec (); - F77_XFCN (zpocon, ZPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, - n, anorm, xrcond, pz, prz, zpocon_info - F77_CHAR_ARG_LEN (1))); + // If someone thinks of a more graceful way of doing this (or + // faster for that matter :-)), please let me know! - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zpocon"); - - if (zpocon_info != 0) - info = -1; - } - else - { - // If someone thinks of a more graceful way of doing this (or - // faster for that matter :-)), please let me know! - - if (n > 1) - for (octave_idx_type j = 0; j < a_nc; j++) - for (octave_idx_type i = j+1; i < a_nr; i++) - chol_mat.xelem (i, j) = 0.0; - } + if (n > 1) + for (octave_idx_type j = 0; j < a_nc; j++) + for (octave_idx_type i = j+1; i < a_nr; i++) + chol_mat.xelem (i, j) = 0.0; } return info; @@ -137,20 +128,15 @@ tmp.fortran_vec (), n, info F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in zpotri"); - else - { - // If someone thinks of a more graceful way of doing this (or - // faster for that matter :-)), please let me know! + // If someone thinks of a more graceful way of doing this (or + // faster for that matter :-)), please let me know! - if (n > 1) - for (octave_idx_type j = 0; j < r_nc; j++) - for (octave_idx_type i = j+1; i < r_nr; i++) - tmp.xelem (i, j) = std::conj (tmp.xelem (j, i)); + if (n > 1) + for (octave_idx_type j = 0; j < r_nc; j++) + for (octave_idx_type i = j+1; i < r_nr; i++) + tmp.xelem (i, j) = std::conj (tmp.xelem (j, i)); - retval = tmp; - } + retval = tmp; } else (*current_liboctave_error_handler) ("chol2inv requires square matrix");
--- a/liboctave/CmplxHESS.cc Thu Feb 14 17:14:23 2008 -0500 +++ b/liboctave/CmplxHESS.cc Thu Feb 14 21:57:50 2008 -0500 @@ -88,56 +88,34 @@ n, h, n, ilo, ihi, pscale, info F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in zgebal"); - else - { - Array<Complex> tau (n-1); - Complex *ptau = tau.fortran_vec (); + Array<Complex> tau (n-1); + Complex *ptau = tau.fortran_vec (); - Array<Complex> work (lwork); - Complex *pwork = work.fortran_vec (); + Array<Complex> work (lwork); + Complex *pwork = work.fortran_vec (); - F77_XFCN (zgehrd, ZGEHRD, (n, ilo, ihi, h, n, ptau, pwork, lwork, info)); + F77_XFCN (zgehrd, ZGEHRD, (n, ilo, ihi, h, n, ptau, pwork, lwork, info)); - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in zgehrd"); - else - { - unitary_hess_mat = hess_mat; - Complex *z = unitary_hess_mat.fortran_vec (); + unitary_hess_mat = hess_mat; + Complex *z = unitary_hess_mat.fortran_vec (); - F77_XFCN (zunghr, ZUNGHR, (n, ilo, ihi, z, n, ptau, pwork, - lwork, info)); + F77_XFCN (zunghr, ZUNGHR, (n, ilo, ihi, z, n, ptau, pwork, + lwork, info)); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zunghr"); - else - { - F77_XFCN (zgebak, ZGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1), - F77_CONST_CHAR_ARG2 (&side, 1), - n, ilo, ihi, pscale, n, z, n, info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); + F77_XFCN (zgebak, ZGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1), + F77_CONST_CHAR_ARG2 (&side, 1), + n, ilo, ihi, pscale, n, z, n, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zgebak"); - else - { - // If someone thinks of a more graceful way of - // doing this (or faster for that matter :-)), - // please let me know! + // If someone thinks of a more graceful way of + // doing this (or faster for that matter :-)), + // please let me know! - if (n > 2) - for (octave_idx_type j = 0; j < a_nc; j++) - for (octave_idx_type i = j+2; i < a_nr; i++) - hess_mat.elem (i, j) = 0; - } - } - } - } + if (n > 2) + for (octave_idx_type j = 0; j < a_nc; j++) + for (octave_idx_type i = j+2; i < a_nr; i++) + hess_mat.elem (i, j) = 0; return info; }
--- a/liboctave/CmplxLU.cc Thu Feb 14 17:14:23 2008 -0500 +++ b/liboctave/CmplxLU.cc Thu Feb 14 21:57:50 2008 -0500 @@ -61,10 +61,7 @@ F77_XFCN (zgetrf, ZGETRF, (a_nr, a_nc, tmp_data, a_nr, pipvt, info)); - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in zgetrf"); - else - ipvt -= static_cast<octave_idx_type> (1); + ipvt -= static_cast<octave_idx_type> (1); } /*
--- a/liboctave/CmplxQR.cc Thu Feb 14 17:14:23 2008 -0500 +++ b/liboctave/CmplxQR.cc Thu Feb 14 21:57:50 2008 -0500 @@ -85,56 +85,45 @@ F77_XFCN (zgeqrf, ZGEQRF, (m, n, tmp_data, m, ptau, pwork, lwork, info)); - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in zgeqrf"); + if (qr_type == QR::raw) + { + for (octave_idx_type j = 0; j < min_mn; j++) + { + octave_idx_type limit = j < min_mn - 1 ? j : min_mn - 1; + for (octave_idx_type i = limit + 1; i < m; i++) + A_fact.elem (i, j) *= tau.elem (j); + } + + r = A_fact; + + if (m > n) + r.resize (m, n); + } else { - if (qr_type == QR::raw) - { - for (octave_idx_type j = 0; j < min_mn; j++) - { - octave_idx_type limit = j < min_mn - 1 ? j : min_mn - 1; - for (octave_idx_type i = limit + 1; i < m; i++) - A_fact.elem (i, j) *= tau.elem (j); - } - - r = A_fact; + octave_idx_type n2 = (qr_type == QR::economy) ? min_mn : m; - if (m > n) - r.resize (m, n); - } + if (qr_type == QR::economy && m > n) + r.resize (n, n, 0.0); else - { - octave_idx_type n2 = (qr_type == QR::economy) ? min_mn : m; - - if (qr_type == QR::economy && m > n) - r.resize (n, n, 0.0); - else - r.resize (m, n, 0.0); + r.resize (m, n, 0.0); - for (octave_idx_type j = 0; j < n; j++) - { - octave_idx_type limit = j < min_mn-1 ? j : min_mn-1; - for (octave_idx_type i = 0; i <= limit; i++) - r.elem (i, j) = A_fact.elem (i, j); - } - - lwork = 32 * n2; - work.resize (lwork); - Complex *pwork2 = work.fortran_vec (); + for (octave_idx_type j = 0; j < n; j++) + { + octave_idx_type limit = j < min_mn-1 ? j : min_mn-1; + for (octave_idx_type i = 0; i <= limit; i++) + r.elem (i, j) = A_fact.elem (i, j); + } - F77_XFCN (zungqr, ZUNGQR, (m, n2, min_mn, tmp_data, m, ptau, - pwork2, lwork, info)); + lwork = 32 * n2; + work.resize (lwork); + Complex *pwork2 = work.fortran_vec (); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in zungqr"); - else - { - q = A_fact; - q.resize (m, n2); - } - } + F77_XFCN (zungqr, ZUNGQR, (m, n2, min_mn, tmp_data, m, ptau, + pwork2, lwork, info)); + + q = A_fact; + q.resize (m, n2); } }
--- a/liboctave/CmplxQRP.cc Thu Feb 14 17:14:23 2008 -0500 +++ b/liboctave/CmplxQRP.cc Thu Feb 14 21:57:50 2008 -0500 @@ -94,51 +94,41 @@ F77_XFCN (zgeqpf, ZGEQPF, (m, n, tmp_data, m, pjpvt, ptau, pwork, prwork, info)); - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in zgeqpf"); + // Form Permutation matrix (if economy is requested, return the + // indices only!) + + if (qr_type == QR::economy) + { + p.resize (1, n, 0.0); + for (octave_idx_type j = 0; j < n; j++) + p.elem (0, j) = jpvt.elem (j); + } else { - // Form Permutation matrix (if economy is requested, return the - // indices only!) + p.resize (n, n, 0.0); + for (octave_idx_type j = 0; j < n; j++) + p.elem (jpvt.elem (j) - 1, j) = 1.0; + } - if (qr_type == QR::economy) - { - p.resize (1, n, 0.0); - for (octave_idx_type j = 0; j < n; j++) - p.elem (0, j) = jpvt.elem (j); - } - else - { - p.resize (n, n, 0.0); - for (octave_idx_type j = 0; j < n; j++) - p.elem (jpvt.elem (j) - 1, j) = 1.0; - } + octave_idx_type n2 = (qr_type == QR::economy) ? min_mn : m; - octave_idx_type n2 = (qr_type == QR::economy) ? min_mn : m; + if (qr_type == QR::economy && m > n) + r.resize (n, n, 0.0); + else + r.resize (m, n, 0.0); - if (qr_type == QR::economy && m > n) - r.resize (n, n, 0.0); - else - r.resize (m, n, 0.0); - - for (octave_idx_type j = 0; j < n; j++) - { - octave_idx_type limit = j < min_mn-1 ? j : min_mn-1; - for (octave_idx_type i = 0; i <= limit; i++) - r.elem (i, j) = A_fact.elem (i, j); - } + for (octave_idx_type j = 0; j < n; j++) + { + octave_idx_type limit = j < min_mn-1 ? j : min_mn-1; + for (octave_idx_type i = 0; i <= limit; i++) + r.elem (i, j) = A_fact.elem (i, j); + } - F77_XFCN (zungqr, ZUNGQR, (m, n2, min_mn, tmp_data, m, ptau, - pwork, lwork, info)); + F77_XFCN (zungqr, ZUNGQR, (m, n2, min_mn, tmp_data, m, ptau, + pwork, lwork, info)); - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in zungqr"); - else - { - q = A_fact; - q.resize (m, n2); - } - } + q = A_fact; + q.resize (m, n2); } /*
--- a/liboctave/CmplxSCHUR.cc Thu Feb 14 17:14:23 2008 -0500 +++ b/liboctave/CmplxSCHUR.cc Thu Feb 14 21:57:50 2008 -0500 @@ -132,9 +132,6 @@ F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in zgeesx"); - return info; }
--- a/liboctave/CmplxSVD.cc Thu Feb 14 17:14:23 2008 -0500 +++ b/liboctave/CmplxSVD.cc Thu Feb 14 21:57:50 2008 -0500 @@ -149,29 +149,19 @@ F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in zgesvd"); - else - { - lwork = static_cast<octave_idx_type> (work(0).real ()); - work.resize (lwork); + lwork = static_cast<octave_idx_type> (work(0).real ()); + work.resize (lwork); - F77_XFCN (zgesvd, ZGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), - F77_CONST_CHAR_ARG2 (&jobv, 1), - m, n, tmp_data, m, s_vec, u, m, vt, - nrow_vt, work.fortran_vec (), lwork, - rwork.fortran_vec (), info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); + F77_XFCN (zgesvd, ZGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), + F77_CONST_CHAR_ARG2 (&jobv, 1), + m, n, tmp_data, m, s_vec, u, m, vt, + nrow_vt, work.fortran_vec (), lwork, + rwork.fortran_vec (), info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in zgesvd"); - else - { - if (! (jobv == 'N' || jobv == 'O')) - right_sm = right_sm.hermitian (); - } - } + if (! (jobv == 'N' || jobv == 'O')) + right_sm = right_sm.hermitian (); return info; }
--- a/liboctave/DASPK.cc Thu Feb 14 17:14:23 2008 -0500 +++ b/liboctave/DASPK.cc Thu Feb 14 21:57:50 2008 -0500 @@ -462,67 +462,59 @@ piwork, liw, dummy, idummy, ddaspk_j, ddaspk_psol)); - if (f77_exception_encountered) - { - integration_error = true; - (*current_liboctave_error_handler) ("unrecoverable error in daspk"); - } - else + switch (istate) { - switch (istate) - { - case 1: // A step was successfully taken in intermediate-output - // mode. The code has not yet reached TOUT. - case 2: // The integration to TSTOP was successfully completed - // (T=TSTOP) by stepping exactly to TSTOP. - case 3: // The integration to TOUT was successfully completed - // (T=TOUT) by stepping past TOUT. Y(*) is obtained by - // interpolation. YPRIME(*) is obtained by interpolation. - case 4: // The initial condition calculation, with - // INFO(11) > 0, was successful, and INFO(14) = 1. - // No integration steps were taken, and the solution - // is not considered to have been started. - retval = x; - t = tout; - break; + case 1: // A step was successfully taken in intermediate-output + // mode. The code has not yet reached TOUT. + case 2: // The integration to TSTOP was successfully completed + // (T=TSTOP) by stepping exactly to TSTOP. + case 3: // The integration to TOUT was successfully completed + // (T=TOUT) by stepping past TOUT. Y(*) is obtained by + // interpolation. YPRIME(*) is obtained by interpolation. + case 4: // The initial condition calculation, with + // INFO(11) > 0, was successful, and INFO(14) = 1. + // No integration steps were taken, and the solution + // is not considered to have been started. + retval = x; + t = tout; + break; - case -1: // A large amount of work has been expended. (~500 steps). - case -2: // The error tolerances are too stringent. - case -3: // The local error test cannot be satisfied because you - // specified a zero component in ATOL and the - // corresponding computed solution component is zero. - // Thus, a pure relative error test is impossible for - // this component. - case -6: // DDASPK had repeated error test failures on the last - // attempted step. - case -7: // The corrector could not converge. - case -8: // The matrix of partial derivatives is singular. - case -9: // The corrector could not converge. There were repeated - // error test failures in this step. - case -10: // The corrector could not converge because IRES was - // equal to minus one. - case -11: // IRES equal to -2 was encountered and control is being - // returned to the calling program. - case -12: // DDASPK failed to compute the initial YPRIME. - case -13: // Unrecoverable error encountered inside user's - // PSOL routine, and control is being returned to - // the calling program. - case -14: // The Krylov linear system solver could not - // achieve convergence. - case -33: // The code has encountered trouble from which it cannot - // recover. A message is printed explaining the trouble - // and control is returned to the calling program. For - // example, this occurs when invalid input is detected. - integration_error = true; - break; + case -1: // A large amount of work has been expended. (~500 steps). + case -2: // The error tolerances are too stringent. + case -3: // The local error test cannot be satisfied because you + // specified a zero component in ATOL and the + // corresponding computed solution component is zero. + // Thus, a pure relative error test is impossible for + // this component. + case -6: // DDASPK had repeated error test failures on the last + // attempted step. + case -7: // The corrector could not converge. + case -8: // The matrix of partial derivatives is singular. + case -9: // The corrector could not converge. There were repeated + // error test failures in this step. + case -10: // The corrector could not converge because IRES was + // equal to minus one. + case -11: // IRES equal to -2 was encountered and control is being + // returned to the calling program. + case -12: // DDASPK failed to compute the initial YPRIME. + case -13: // Unrecoverable error encountered inside user's + // PSOL routine, and control is being returned to + // the calling program. + case -14: // The Krylov linear system solver could not + // achieve convergence. + case -33: // The code has encountered trouble from which it cannot + // recover. A message is printed explaining the trouble + // and control is returned to the calling program. For + // example, this occurs when invalid input is detected. + integration_error = true; + break; - default: - integration_error = true; - (*current_liboctave_error_handler) - ("unrecognized value of istate (= %d) returned from ddaspk", - istate); - break; - } + default: + integration_error = true; + (*current_liboctave_error_handler) + ("unrecognized value of istate (= %d) returned from ddaspk", + istate); + break; } return retval;
--- a/liboctave/DASRT.cc Thu Feb 14 17:14:23 2008 -0500 +++ b/liboctave/DASRT.cc Thu Feb 14 21:57:50 2008 -0500 @@ -324,61 +324,53 @@ piwork, liw, dummy, idummy, ddasrt_j, ddasrt_g, ng, pjroot)); - if (f77_exception_encountered) - { - integration_error = true; - (*current_liboctave_error_handler) ("unrecoverable error in dasrt"); - } - else + switch (istate) { - switch (istate) - { - case 1: // A step was successfully taken in intermediate-output - // mode. The code has not yet reached TOUT. - case 2: // The integration to TOUT was successfully completed - // (T=TOUT) by stepping exactly to TOUT. - case 3: // The integration to TOUT was successfully completed - // (T=TOUT) by stepping past TOUT. Y(*) is obtained by - // interpolation. YPRIME(*) is obtained by interpolation. - t = tout; - break; + case 1: // A step was successfully taken in intermediate-output + // mode. The code has not yet reached TOUT. + case 2: // The integration to TOUT was successfully completed + // (T=TOUT) by stepping exactly to TOUT. + case 3: // The integration to TOUT was successfully completed + // (T=TOUT) by stepping past TOUT. Y(*) is obtained by + // interpolation. YPRIME(*) is obtained by interpolation. + t = tout; + break; - case 4: // The integration was successfully completed - // by finding one or more roots of G at T. - break; + case 4: // The integration was successfully completed + // by finding one or more roots of G at T. + break; - case -1: // A large amount of work has been expended. - case -2: // The error tolerances are too stringent. - case -3: // The local error test cannot be satisfied because you - // specified a zero component in ATOL and the - // corresponding computed solution component is zero. - // Thus, a pure relative error test is impossible for - // this component. - case -6: // DDASRT had repeated error test failures on the last - // attempted step. - case -7: // The corrector could not converge. - case -8: // The matrix of partial derivatives is singular. - case -9: // The corrector could not converge. There were repeated - // error test failures in this step. - case -10: // The corrector could not converge because IRES was - // equal to minus one. - case -11: // IRES equal to -2 was encountered and control is being - // returned to the calling program. - case -12: // DASSL failed to compute the initial YPRIME. - case -33: // The code has encountered trouble from which it cannot - // recover. A message is printed explaining the trouble - // and control is returned to the calling program. For - // example, this occurs when invalid input is detected. - integration_error = true; - break; + case -1: // A large amount of work has been expended. + case -2: // The error tolerances are too stringent. + case -3: // The local error test cannot be satisfied because you + // specified a zero component in ATOL and the + // corresponding computed solution component is zero. + // Thus, a pure relative error test is impossible for + // this component. + case -6: // DDASRT had repeated error test failures on the last + // attempted step. + case -7: // The corrector could not converge. + case -8: // The matrix of partial derivatives is singular. + case -9: // The corrector could not converge. There were repeated + // error test failures in this step. + case -10: // The corrector could not converge because IRES was + // equal to minus one. + case -11: // IRES equal to -2 was encountered and control is being + // returned to the calling program. + case -12: // DASSL failed to compute the initial YPRIME. + case -33: // The code has encountered trouble from which it cannot + // recover. A message is printed explaining the trouble + // and control is returned to the calling program. For + // example, this occurs when invalid input is detected. + integration_error = true; + break; - default: - integration_error = true; - (*current_liboctave_error_handler) - ("unrecognized value of istate (= %d) returned from ddasrt", - istate); - break; - } + default: + integration_error = true; + (*current_liboctave_error_handler) + ("unrecognized value of istate (= %d) returned from ddasrt", + istate); + break; } }
--- a/liboctave/DASSL.cc Thu Feb 14 17:14:23 2008 -0500 +++ b/liboctave/DASSL.cc Thu Feb 14 21:57:50 2008 -0500 @@ -287,58 +287,50 @@ prel_tol, pabs_tol, istate, prwork, lrw, piwork, liw, dummy, idummy, ddassl_j)); - if (f77_exception_encountered) - { - integration_error = true; - (*current_liboctave_error_handler) ("unrecoverable error in dassl"); - } - else + switch (istate) { - switch (istate) - { - case 1: // A step was successfully taken in intermediate-output - // mode. The code has not yet reached TOUT. - case 2: // The integration to TSTOP was successfully completed - // (T=TSTOP) by stepping exactly to TSTOP. - case 3: // The integration to TOUT was successfully completed - // (T=TOUT) by stepping past TOUT. Y(*) is obtained by - // interpolation. YPRIME(*) is obtained by interpolation. - retval = x; - t = tout; - break; + case 1: // A step was successfully taken in intermediate-output + // mode. The code has not yet reached TOUT. + case 2: // The integration to TSTOP was successfully completed + // (T=TSTOP) by stepping exactly to TSTOP. + case 3: // The integration to TOUT was successfully completed + // (T=TOUT) by stepping past TOUT. Y(*) is obtained by + // interpolation. YPRIME(*) is obtained by interpolation. + retval = x; + t = tout; + break; - case -1: // A large amount of work has been expended. (~500 steps). - case -2: // The error tolerances are too stringent. - case -3: // The local error test cannot be satisfied because you - // specified a zero component in ATOL and the - // corresponding computed solution component is zero. - // Thus, a pure relative error test is impossible for - // this component. - case -6: // DDASSL had repeated error test failures on the last - // attempted step. - case -7: // The corrector could not converge. - case -8: // The matrix of partial derivatives is singular. - case -9: // The corrector could not converge. There were repeated - // error test failures in this step. - case -10: // The corrector could not converge because IRES was - // equal to minus one. - case -11: // IRES equal to -2 was encountered and control is being - // returned to the calling program. - case -12: // DDASSL failed to compute the initial YPRIME. - case -33: // The code has encountered trouble from which it cannot - // recover. A message is printed explaining the trouble - // and control is returned to the calling program. For - // example, this occurs when invalid input is detected. - integration_error = true; - break; + case -1: // A large amount of work has been expended. (~500 steps). + case -2: // The error tolerances are too stringent. + case -3: // The local error test cannot be satisfied because you + // specified a zero component in ATOL and the + // corresponding computed solution component is zero. + // Thus, a pure relative error test is impossible for + // this component. + case -6: // DDASSL had repeated error test failures on the last + // attempted step. + case -7: // The corrector could not converge. + case -8: // The matrix of partial derivatives is singular. + case -9: // The corrector could not converge. There were repeated + // error test failures in this step. + case -10: // The corrector could not converge because IRES was + // equal to minus one. + case -11: // IRES equal to -2 was encountered and control is being + // returned to the calling program. + case -12: // DDASSL failed to compute the initial YPRIME. + case -33: // The code has encountered trouble from which it cannot + // recover. A message is printed explaining the trouble + // and control is returned to the calling program. For + // example, this occurs when invalid input is detected. + integration_error = true; + break; - default: - integration_error = true; - (*current_liboctave_error_handler) - ("unrecognized value of istate (= %d) returned from ddassl", - istate); - break; - } + default: + integration_error = true; + (*current_liboctave_error_handler) + ("unrecognized value of istate (= %d) returned from ddassl", + istate); + break; } return retval;
--- a/liboctave/EIG.cc Thu Feb 14 17:14:23 2008 -0500 +++ b/liboctave/EIG.cc Thu Feb 14 21:57:50 2008 -0500 @@ -116,7 +116,7 @@ F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - if (! f77_exception_encountered && info == 0) + if (info == 0) { lwork = static_cast<octave_idx_type> (dummy_work); Array<double> work (lwork); @@ -129,7 +129,7 @@ F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered || info < 0) + if (info < 0) { (*current_liboctave_error_handler) ("unrecoverable error in dgeev"); return info; @@ -208,7 +208,7 @@ F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - if (! f77_exception_encountered && info == 0) + if (info == 0) { lwork = static_cast<octave_idx_type> (dummy_work); Array<double> work (lwork); @@ -220,7 +220,7 @@ F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered || info < 0) + if (info < 0) { (*current_liboctave_error_handler) ("unrecoverable error in dsyev"); return info; @@ -291,7 +291,7 @@ F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - if (! f77_exception_encountered && info == 0) + if (info == 0) { lwork = static_cast<octave_idx_type> (dummy_work.real ()); Array<Complex> work (lwork); @@ -304,7 +304,7 @@ F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered || info < 0) + if (info < 0) { (*current_liboctave_error_handler) ("unrecoverable error in zgeev"); return info; @@ -358,7 +358,7 @@ F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - if (! f77_exception_encountered && info == 0) + if (info == 0) { lwork = static_cast<octave_idx_type> (dummy_work.real ()); Array<Complex> work (lwork); @@ -370,7 +370,7 @@ F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered || info < 0) + if (info < 0) { (*current_liboctave_error_handler) ("unrecoverable error in zheev"); return info;
--- a/liboctave/LSODE.cc Thu Feb 14 17:14:23 2008 -0500 +++ b/liboctave/LSODE.cc Thu Feb 14 21:57:50 2008 -0500 @@ -275,40 +275,32 @@ pabs_tol, itask, istate, iopt, prwork, lrw, piwork, liw, lsode_j, method_flag)); - if (f77_exception_encountered) - { - integration_error = true; - (*current_liboctave_error_handler) ("unrecoverable error in lsode"); - } - else + switch (istate) { - switch (istate) - { - case 1: // prior to initial integration step. - case 2: // lsode was successful. - retval = x; - t = tout; - break; - - case -1: // excess work done on this call (perhaps wrong mf). - case -2: // excess accuracy requested (tolerances too small). - case -3: // illegal input detected (see printed message). - case -4: // repeated error test failures (check all inputs). - case -5: // repeated convergence failures (perhaps bad jacobian - // supplied or wrong choice of mf or tolerances). - case -6: // error weight became zero during problem. (solution - // component i vanished, and atol or atol(i) = 0.) - case -13: // return requested in user-supplied function. - integration_error = true; - break; + case 1: // prior to initial integration step. + case 2: // lsode was successful. + retval = x; + t = tout; + break; - default: - integration_error = true; - (*current_liboctave_error_handler) - ("unrecognized value of istate (= %d) returned from lsode", - istate); - break; - } + case -1: // excess work done on this call (perhaps wrong mf). + case -2: // excess accuracy requested (tolerances too small). + case -3: // illegal input detected (see printed message). + case -4: // repeated error test failures (check all inputs). + case -5: // repeated convergence failures (perhaps bad jacobian + // supplied or wrong choice of mf or tolerances). + case -6: // error weight became zero during problem. (solution + // component i vanished, and atol or atol(i) = 0.) + case -13: // return requested in user-supplied function. + integration_error = true; + break; + + default: + integration_error = true; + (*current_liboctave_error_handler) + ("unrecognized value of istate (= %d) returned from lsode", + istate); + break; } return retval;
--- a/liboctave/NLEqn.cc Thu Feb 14 17:14:23 2008 -0500 +++ b/liboctave/NLEqn.cc Thu Feb 14 21:57:50 2008 -0500 @@ -176,10 +176,7 @@ solution_status = info; - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in hybrj1"); - else - fval = ColumnVector (fvec); + fval = ColumnVector (fvec); } else { @@ -195,10 +192,7 @@ solution_status = info; - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in hybrd1"); - else - fval = ColumnVector (fvec); + fval = ColumnVector (fvec); } return retval;
--- a/liboctave/Quad.cc Thu Feb 14 17:14:23 2008 -0500 +++ b/liboctave/Quad.cc Thu Feb 14 21:57:50 2008 -0500 @@ -112,9 +112,6 @@ abserr, neval, ier, leniw, lenw, last, piwork, pwork)); - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in dqagp"); - return result; } @@ -161,9 +158,6 @@ result, abserr, neval, ier, leniw, lenw, last, piwork, pwork)); - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in dqagi"); - return result; }
--- a/liboctave/dColVector.cc Thu Feb 14 17:14:23 2008 -0500 +++ b/liboctave/dColVector.cc Thu Feb 14 21:57:50 2008 -0500 @@ -221,10 +221,6 @@ nr, nc, 1.0, m.data (), ld, a.data (), 1, 0.0, y, 1 F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dgemv"); } }
--- a/liboctave/dMatrix.cc Thu Feb 14 17:14:23 2008 -0500 +++ b/liboctave/dMatrix.cc Thu Feb 14 21:57:50 2008 -0500 @@ -711,38 +711,29 @@ F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in dtrtri"); - else + // Throw-away extra info LAPACK gives so as to not change output. + rcond = 0.0; + if (info != 0) + info = -1; + else if (calc_cond) { - // Throw-away extra info LAPACK gives so as to not change output. - rcond = 0.0; - if (info != 0) + octave_idx_type dtrcon_info = 0; + char job = '1'; + + OCTAVE_LOCAL_BUFFER (double, work, 3 * nr); + OCTAVE_LOCAL_BUFFER (octave_idx_type, iwork, nr); + + F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&job, 1), + F77_CONST_CHAR_ARG2 (&uplo, 1), + F77_CONST_CHAR_ARG2 (&udiag, 1), + nr, tmp_data, nr, rcond, + work, iwork, dtrcon_info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (dtrcon_info != 0) info = -1; - else if (calc_cond) - { - octave_idx_type dtrcon_info = 0; - char job = '1'; - - OCTAVE_LOCAL_BUFFER (double, work, 3 * nr); - OCTAVE_LOCAL_BUFFER (octave_idx_type, iwork, nr); - - F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&job, 1), - F77_CONST_CHAR_ARG2 (&uplo, 1), - F77_CONST_CHAR_ARG2 (&udiag, 1), - nr, tmp_data, nr, rcond, - work, iwork, dtrcon_info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dtrcon"); - - if (dtrcon_info != 0) - info = -1; - } } if (info == -1 && ! force) @@ -779,13 +770,6 @@ F77_XFCN (dgetri, DGETRI, (nc, tmp_data, nr, pipvt, z.fortran_vec (), lwork, info)); - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("unrecoverable error in dgetri"); - return retval; - } - lwork = static_cast<octave_idx_type> (z(0)); lwork = (lwork < 2 *nc ? 2*nc : lwork); z.resize (lwork); @@ -800,51 +784,38 @@ F77_XFCN (dgetrf, DGETRF, (nc, nc, tmp_data, nr, pipvt, info)); - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in dgetrf"); + // Throw-away extra info LAPACK gives so as to not change output. + rcond = 0.0; + if (info != 0) + info = -1; + else if (calc_cond) + { + octave_idx_type dgecon_info = 0; + + // Now calculate the condition number for non-singular matrix. + char job = '1'; + Array<octave_idx_type> iz (nc); + octave_idx_type *piz = iz.fortran_vec (); + F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1), + nc, tmp_data, nr, anorm, + rcond, pz, piz, dgecon_info + F77_CHAR_ARG_LEN (1))); + + if (dgecon_info != 0) + info = -1; + } + + if (info == -1 && ! force) + retval = *this; // Restore matrix contents. else { - // Throw-away extra info LAPACK gives so as to not change output. - rcond = 0.0; - if (info != 0) + octave_idx_type dgetri_info = 0; + + F77_XFCN (dgetri, DGETRI, (nc, tmp_data, nr, pipvt, + pz, lwork, dgetri_info)); + + if (dgetri_info != 0) info = -1; - else if (calc_cond) - { - octave_idx_type dgecon_info = 0; - - // Now calculate the condition number for non-singular matrix. - char job = '1'; - Array<octave_idx_type> iz (nc); - octave_idx_type *piz = iz.fortran_vec (); - F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1), - nc, tmp_data, nr, anorm, - rcond, pz, piz, dgecon_info - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dgecon"); - - if (dgecon_info != 0) - info = -1; - } - - if (info == -1 && ! force) - retval = *this; // Restore matrix contents. - else - { - octave_idx_type dgetri_info = 0; - - F77_XFCN (dgetri, DGETRI, (nc, tmp_data, nr, pipvt, - pz, lwork, dgetri_info)); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dgetri"); - - if (dgetri_info != 0) - info = -1; - } } if (info != 0) @@ -1284,12 +1255,30 @@ F77_XFCN (dgetrf, DGETRF, (nr, nr, tmp_data, nr, pipvt, info)); - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in dgetrf"); - else + // Throw-away extra info LAPACK gives so as to not change output. + rcond = 0.0; + if (info != 0) + { + info = -1; + retval = DET (); + } + else { - // Throw-away extra info LAPACK gives so as to not change output. - rcond = 0.0; + if (calc_cond) + { + // Now calc the condition number for non-singular matrix. + char job = '1'; + Array<double> z (4 * nc); + double *pz = z.fortran_vec (); + Array<octave_idx_type> iz (nc); + octave_idx_type *piz = iz.fortran_vec (); + + F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1), + nc, tmp_data, nr, anorm, + rcond, pz, piz, info + F77_CHAR_ARG_LEN (1))); + } + if (info != 0) { info = -1; @@ -1297,60 +1286,33 @@ } else { - if (calc_cond) - { - // Now calc the condition number for non-singular matrix. - char job = '1'; - Array<double> z (4 * nc); - double *pz = z.fortran_vec (); - Array<octave_idx_type> iz (nc); - octave_idx_type *piz = iz.fortran_vec (); - - F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1), - nc, tmp_data, nr, anorm, - rcond, pz, piz, info - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dgecon"); - } - - if (info != 0) - { - info = -1; - retval = DET (); - } - else + double c = 1.0; + int e = 0; + + for (octave_idx_type i = 0; i < nc; i++) { - double c = 1.0; - int e = 0; - - for (octave_idx_type i = 0; i < nc; i++) + if (ipvt(i) != (i+1)) + c = -c; + + c *= atmp(i,i); + + if (c == 0.0) + break; + + while (fabs (c) < 0.5) { - if (ipvt(i) != (i+1)) - c = -c; - - c *= atmp(i,i); - - if (c == 0.0) - break; - - while (fabs (c) < 0.5) - { - c *= 2.0; - e--; - } - - while (fabs (c) >= 2.0) - { - c /= 2.0; - e++; - } + c *= 2.0; + e--; } - retval = DET (c, e); + while (fabs (c) >= 2.0) + { + c /= 2.0; + e++; + } } + + retval = DET (c, e); } } } @@ -1413,10 +1375,6 @@ F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dtrcon"); - if (info != 0) info = -2; @@ -1452,10 +1410,6 @@ F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dtrtrs"); } } } @@ -1521,10 +1475,6 @@ F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dtrcon"); - if (info != 0) info = -2; @@ -1560,10 +1510,6 @@ F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dtrtrs"); } } } @@ -1608,78 +1554,64 @@ tmp_data, nr, info F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dpotrf"); - else + // Throw-away extra info LAPACK gives so as to not change output. + rcond = 0.0; + if (info != 0) + { + info = -2; + + mattype.mark_as_unsymmetric (); + typ = MatrixType::Full; + } + else { - // Throw-away extra info LAPACK gives so as to not change output. - rcond = 0.0; - if (info != 0) + if (calc_cond) { - info = -2; - + Array<double> z (3 * nc); + double *pz = z.fortran_vec (); + Array<octave_idx_type> iz (nc); + octave_idx_type *piz = iz.fortran_vec (); + + F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 (&job, 1), + nr, tmp_data, nr, anorm, + rcond, pz, piz, info + F77_CHAR_ARG_LEN (1))); + + if (info != 0) + info = -2; + + volatile double rcond_plus_one = rcond + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcond)) + { + info = -2; + + if (sing_handler) + sing_handler (rcond); + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision, rcond = %g", + rcond); + } + } + + if (info == 0) + { + retval = b; + double *result = retval.fortran_vec (); + + octave_idx_type b_nc = b.cols (); + + F77_XFCN (dpotrs, DPOTRS, (F77_CONST_CHAR_ARG2 (&job, 1), + nr, b_nc, tmp_data, nr, + result, b.rows(), info + F77_CHAR_ARG_LEN (1))); + } + else + { mattype.mark_as_unsymmetric (); typ = MatrixType::Full; - } - else - { - if (calc_cond) - { - Array<double> z (3 * nc); - double *pz = z.fortran_vec (); - Array<octave_idx_type> iz (nc); - octave_idx_type *piz = iz.fortran_vec (); - - F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 (&job, 1), - nr, tmp_data, nr, anorm, - rcond, pz, piz, info - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dpocon"); - - if (info != 0) - info = -2; - - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) - { - info = -2; - - if (sing_handler) - sing_handler (rcond); - else - (*current_liboctave_error_handler) - ("matrix singular to machine precision, rcond = %g", - rcond); - } - } - - if (info == 0) - { - retval = b; - double *result = retval.fortran_vec (); - - octave_idx_type b_nc = b.cols (); - - F77_XFCN (dpotrs, DPOTRS, (F77_CONST_CHAR_ARG2 (&job, 1), - nr, b_nc, tmp_data, nr, - result, b.rows(), info - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dpotrs"); - } - else - { - mattype.mark_as_unsymmetric (); - typ = MatrixType::Full; - } - } + } } } @@ -1702,79 +1634,65 @@ F77_XFCN (dgetrf, DGETRF, (nr, nr, tmp_data, nr, pipvt, info)); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dgetrf"); - else + // Throw-away extra info LAPACK gives so as to not change output. + rcond = 0.0; + if (info != 0) { - // Throw-away extra info LAPACK gives so as to not change output. - rcond = 0.0; - if (info != 0) + info = -2; + + if (sing_handler) + sing_handler (rcond); + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision"); + + mattype.mark_as_rectangular (); + } + else + { + if (calc_cond) { - info = -2; - - if (sing_handler) - sing_handler (rcond); - else - (*current_liboctave_error_handler) - ("matrix singular to machine precision"); - - mattype.mark_as_rectangular (); - } - else - { - if (calc_cond) + // Now calculate the condition number for + // non-singular matrix. + char job = '1'; + F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1), + nc, tmp_data, nr, anorm, + rcond, pz, piz, info + F77_CHAR_ARG_LEN (1))); + + if (info != 0) + info = -2; + + volatile double rcond_plus_one = rcond + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcond)) { - // Now calculate the condition number for - // non-singular matrix. - char job = '1'; - F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1), - nc, tmp_data, nr, anorm, - rcond, pz, piz, info - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dgecon"); - - if (info != 0) - info = -2; - - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) - { - info = -2; - - if (sing_handler) - sing_handler (rcond); - else - (*current_liboctave_error_handler) - ("matrix singular to machine precision, rcond = %g", - rcond); - } + info = -2; + + if (sing_handler) + sing_handler (rcond); + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision, rcond = %g", + rcond); } - - if (info == 0) - { - retval = b; - double *result = retval.fortran_vec (); - - octave_idx_type b_nc = b.cols (); - - char job = 'N'; - F77_XFCN (dgetrs, DGETRS, (F77_CONST_CHAR_ARG2 (&job, 1), - nr, b_nc, tmp_data, nr, - pipvt, result, b.rows(), info - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dgetrs"); - } - else - mattype.mark_as_rectangular (); } + + if (info == 0) + { + retval = b; + double *result = retval.fortran_vec (); + + octave_idx_type b_nc = b.cols (); + + char job = 'N'; + F77_XFCN (dgetrs, DGETRS, (F77_CONST_CHAR_ARG2 (&job, 1), + nr, b_nc, tmp_data, nr, + pipvt, result, b.rows(), info + F77_CHAR_ARG_LEN (1))); + } + else + mattype.mark_as_rectangular (); } } else if (typ != MatrixType::Hermitian) @@ -2170,35 +2088,23 @@ work(0) = lworkaround; } - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dgelsd"); + lwork = static_cast<octave_idx_type> (work(0)); + work.resize (lwork); + + F77_XFCN (dgelsd, DGELSD, (m, n, nrhs, tmp_data, m, pretval, + maxmn, ps, rcond, rank, + work.fortran_vec (), lwork, + piwork, info)); + + if (rank < minmn) + (*current_liboctave_warning_handler) + ("dgelsd: rank deficient %dx%d matrix, rank = %d", m, n, rank); + if (s.elem (0) == 0.0) + rcond = 0.0; else - { - lwork = static_cast<octave_idx_type> (work(0)); - work.resize (lwork); - - F77_XFCN (dgelsd, DGELSD, (m, n, nrhs, tmp_data, m, pretval, - maxmn, ps, rcond, rank, - work.fortran_vec (), lwork, - piwork, info)); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dgelsd"); - else - { - if (rank < minmn) - (*current_liboctave_warning_handler) - ("dgelsd: rank deficient %dx%d matrix, rank = %d", m, n, rank); - if (s.elem (0) == 0.0) - rcond = 0.0; - else - rcond = s.elem (minmn - 1) / s.elem (0); - - retval.resize (n, nrhs); - } - } + rcond = s.elem (minmn - 1) / s.elem (0); + + retval.resize (n, nrhs); } return retval; @@ -2335,35 +2241,26 @@ ps, rcond, rank, work.fortran_vec (), lwork, piwork, info)); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dgelsd"); - else + lwork = static_cast<octave_idx_type> (work(0)); + work.resize (lwork); + + F77_XFCN (dgelsd, DGELSD, (m, n, nrhs, tmp_data, m, pretval, + maxmn, ps, rcond, rank, + work.fortran_vec (), lwork, + piwork, info)); + + if (rank < minmn) { - lwork = static_cast<octave_idx_type> (work(0)); - work.resize (lwork); - - F77_XFCN (dgelsd, DGELSD, (m, n, nrhs, tmp_data, m, pretval, - maxmn, ps, rcond, rank, - work.fortran_vec (), lwork, - piwork, info)); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dgelsd"); - else if (rank < minmn) - { - if (rank < minmn) - (*current_liboctave_warning_handler) - ("dgelsd: rank deficient %dx%d matrix, rank = %d", m, n, rank); - if (s.elem (0) == 0.0) - rcond = 0.0; - else - rcond = s.elem (minmn - 1) / s.elem (0); - } - - retval.resize (n, nrhs); + if (rank < minmn) + (*current_liboctave_warning_handler) + ("dgelsd: rank deficient %dx%d matrix, rank = %d", m, n, rank); + if (s.elem (0) == 0.0) + rcond = 0.0; + else + rcond = s.elem (minmn - 1) / s.elem (0); } + + retval.resize (n, nrhs); } return retval; @@ -2479,12 +2376,6 @@ dscale.fortran_vec (), info F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) ("unrecoverable error in dgebal"); - return retval; - } - // Preconditioning step 3: scaling. ColumnVector work(nc); @@ -2495,12 +2386,6 @@ work.fortran_vec (), inf_norm F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) ("unrecoverable error in dlange"); - return retval; - } - octave_idx_type sqpow = static_cast<octave_idx_type> (inf_norm > 0.0 ? (1.0 + log (inf_norm) / log (2.0)) : 0.0); @@ -2723,10 +2608,6 @@ a.data (), 1, 0.0, c, len F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dgemm"); } return retval; @@ -3327,14 +3208,9 @@ F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in dtrsyl"); - else - { - // FIXME -- check info? + // FIXME -- check info? - retval = -ua*cx*ub.transpose (); - } + retval = -ua*cx*ub.transpose (); return retval; } @@ -3393,10 +3269,6 @@ nr, nc, 1.0, m.data (), ld, a.data (), 1, 0.0, c, 1 F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dgemv"); } } else @@ -3407,10 +3279,6 @@ ld, a.data (), lda, 0.0, c, nr F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dgemm"); } } }
--- a/liboctave/dRowVector.cc Thu Feb 14 17:14:23 2008 -0500 +++ b/liboctave/dRowVector.cc Thu Feb 14 21:57:50 2008 -0500 @@ -223,10 +223,6 @@ a_nr, a_nc, 1.0, a.data (), ld, v.data (), 1, 0.0, y, 1 F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dgemv"); } }
--- a/liboctave/dSparse.cc Thu Feb 14 17:14:23 2008 -0500 +++ b/liboctave/dSparse.cc Thu Feb 14 21:57:50 2008 -0500 @@ -3821,10 +3821,7 @@ F77_XFCN (dptsv, DPTSV, (nr, b_nc, D, DL, result, b.rows(), err)); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dptsv"); - else if (err != 0) + if (err != 0) { err = 0; mattype.mark_as_unsymmetric (); @@ -3881,10 +3878,7 @@ F77_XFCN (dgtsv, DGTSV, (nr, b_nc, DL, D, DU, result, b.rows(), err)); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dgtsv"); - else if (err != 0) + if (err != 0) { rcond = 0.; err = -2; @@ -3982,84 +3976,71 @@ F77_XFCN (dgttrf, DGTTRF, (nr, DL, D, DU, DU2, pipvt, err)); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dgttrf"); - else - { - if (err != 0) - { - rcond = 0.0; - err = -2; - - if (sing_handler) - { - sing_handler (rcond); - mattype.mark_as_rectangular (); - } - else - (*current_liboctave_error_handler) - ("matrix singular to machine precision"); - - } - else - { - rcond = 1.0; - char job = 'N'; - volatile octave_idx_type x_nz = b.nnz (); - octave_idx_type b_nc = b.cols (); - retval = SparseMatrix (nr, b_nc, x_nz); - retval.xcidx(0) = 0; - volatile octave_idx_type ii = 0; - - OCTAVE_LOCAL_BUFFER (double, work, nr); - - for (volatile octave_idx_type j = 0; j < b_nc; j++) + if (err != 0) + { + rcond = 0.0; + err = -2; + + if (sing_handler) + { + sing_handler (rcond); + mattype.mark_as_rectangular (); + } + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision"); + + } + else + { + rcond = 1.0; + char job = 'N'; + volatile octave_idx_type x_nz = b.nnz (); + octave_idx_type b_nc = b.cols (); + retval = SparseMatrix (nr, b_nc, x_nz); + retval.xcidx(0) = 0; + volatile octave_idx_type ii = 0; + + OCTAVE_LOCAL_BUFFER (double, work, nr); + + for (volatile octave_idx_type j = 0; j < b_nc; j++) + { + for (octave_idx_type i = 0; i < nr; i++) + work[i] = 0.; + for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++) + work[b.ridx(i)] = b.data(i); + + F77_XFCN (dgttrs, DGTTRS, + (F77_CONST_CHAR_ARG2 (&job, 1), + nr, 1, DL, D, DU, DU2, pipvt, + work, b.rows (), err + F77_CHAR_ARG_LEN (1))); + + // Count non-zeros in work vector and adjust + // space in retval if needed + octave_idx_type new_nnz = 0; + for (octave_idx_type i = 0; i < nr; i++) + if (work[i] != 0.) + new_nnz++; + + if (ii + new_nnz > x_nz) { - for (octave_idx_type i = 0; i < nr; i++) - work[i] = 0.; - for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++) - work[b.ridx(i)] = b.data(i); - - F77_XFCN (dgttrs, DGTTRS, - (F77_CONST_CHAR_ARG2 (&job, 1), - nr, 1, DL, D, DU, DU2, pipvt, - work, b.rows (), err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("unrecoverable error in dgttrs"); - break; - } - - // Count non-zeros in work vector and adjust - // space in retval if needed - octave_idx_type new_nnz = 0; - for (octave_idx_type i = 0; i < nr; i++) - if (work[i] != 0.) - new_nnz++; - - if (ii + new_nnz > x_nz) - { - // Resize the sparse matrix - octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; - retval.change_capacity (sz); - x_nz = sz; - } - - for (octave_idx_type i = 0; i < nr; i++) - if (work[i] != 0.) - { - retval.xridx(ii) = i; - retval.xdata(ii++) = work[i]; - } - retval.xcidx(j+1) = ii; + // Resize the sparse matrix + octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; + retval.change_capacity (sz); + x_nz = sz; } - retval.maybe_compress (); - } + for (octave_idx_type i = 0; i < nr; i++) + if (work[i] != 0.) + { + retval.xridx(ii) = i; + retval.xdata(ii++) = work[i]; + } + retval.xcidx(j+1) = ii; + } + + retval.maybe_compress (); } } else if (typ != MatrixType::Tridiagonal_Hermitian) @@ -4141,13 +4122,7 @@ F77_XFCN (zptsv, ZPTSV, (nr, b_nc, D, DL, result, b_nr, err)); - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("unrecoverable error in zptsv"); - err = -1; - } - else if (err != 0) + if (err != 0) { err = 0; mattype.mark_as_unsymmetric (); @@ -4205,13 +4180,7 @@ F77_XFCN (zgtsv, ZGTSV, (nr, b_nc, DL, D, DU, result, b_nr, err)); - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("unrecoverable error in zgtsv"); - err = -1; - } - else if (err != 0) + if (err != 0) { rcond = 0.; err = -2; @@ -4306,123 +4275,103 @@ F77_XFCN (dgttrf, DGTTRF, (nr, DL, D, DU, DU2, pipvt, err)); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dgttrf"); - else - { - if (err != 0) - { - rcond = 0.0; - err = -2; - - if (sing_handler) + if (err != 0) + { + rcond = 0.0; + err = -2; + + if (sing_handler) + { + sing_handler (rcond); + mattype.mark_as_rectangular (); + } + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision"); + } + else + { + rcond = 1.; + char job = 'N'; + octave_idx_type b_nr = b.rows (); + octave_idx_type b_nc = b.cols (); + OCTAVE_LOCAL_BUFFER (double, Bx, b_nr); + OCTAVE_LOCAL_BUFFER (double, Bz, b_nr); + + // Take a first guess that the number of non-zero terms + // will be as many as in b + volatile octave_idx_type x_nz = b.nnz (); + volatile octave_idx_type ii = 0; + retval = SparseComplexMatrix (b_nr, b_nc, x_nz); + + retval.xcidx(0) = 0; + for (volatile octave_idx_type j = 0; j < b_nc; j++) + { + + for (octave_idx_type i = 0; i < b_nr; i++) { - sing_handler (rcond); - mattype.mark_as_rectangular (); + Complex c = b (i,j); + Bx[i] = std::real (c); + Bz[i] = std::imag (c); } - else - (*current_liboctave_error_handler) - ("matrix singular to machine precision"); - } - else - { - rcond = 1.; - char job = 'N'; - octave_idx_type b_nr = b.rows (); - octave_idx_type b_nc = b.cols (); - OCTAVE_LOCAL_BUFFER (double, Bx, b_nr); - OCTAVE_LOCAL_BUFFER (double, Bz, b_nr); - - // Take a first guess that the number of non-zero terms - // will be as many as in b - volatile octave_idx_type x_nz = b.nnz (); - volatile octave_idx_type ii = 0; - retval = SparseComplexMatrix (b_nr, b_nc, x_nz); - - retval.xcidx(0) = 0; - for (volatile octave_idx_type j = 0; j < b_nc; j++) + + F77_XFCN (dgttrs, DGTTRS, + (F77_CONST_CHAR_ARG2 (&job, 1), + nr, 1, DL, D, DU, DU2, pipvt, + Bx, b_nr, err + F77_CHAR_ARG_LEN (1))); + + if (err != 0) { - - for (octave_idx_type i = 0; i < b_nr; i++) - { - Complex c = b (i,j); - Bx[i] = std::real (c); - Bz[i] = std::imag (c); - } - - F77_XFCN (dgttrs, DGTTRS, - (F77_CONST_CHAR_ARG2 (&job, 1), - nr, 1, DL, D, DU, DU2, pipvt, - Bx, b_nr, err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("unrecoverable error in dgttrs"); - break; - } - - if (err != 0) - { - (*current_liboctave_error_handler) - ("SparseMatrix::solve solve failed"); - - err = -1; - break; - } - - F77_XFCN (dgttrs, DGTTRS, - (F77_CONST_CHAR_ARG2 (&job, 1), - nr, 1, DL, D, DU, DU2, pipvt, - Bz, b_nr, err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("unrecoverable error in dgttrs"); - break; - } - - if (err != 0) - { - (*current_liboctave_error_handler) - ("SparseMatrix::solve solve failed"); - - err = -1; - break; - } - - // Count non-zeros in work vector and adjust - // space in retval if needed - octave_idx_type new_nnz = 0; - for (octave_idx_type i = 0; i < nr; i++) - if (Bx[i] != 0. || Bz[i] != 0.) - new_nnz++; - - if (ii + new_nnz > x_nz) - { - // Resize the sparse matrix - octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; - retval.change_capacity (sz); - x_nz = sz; - } - - for (octave_idx_type i = 0; i < nr; i++) - if (Bx[i] != 0. || Bz[i] != 0.) - { - retval.xridx(ii) = i; - retval.xdata(ii++) = - Complex (Bx[i], Bz[i]); - } - - retval.xcidx(j+1) = ii; + (*current_liboctave_error_handler) + ("SparseMatrix::solve solve failed"); + + err = -1; + break; + } + + F77_XFCN (dgttrs, DGTTRS, + (F77_CONST_CHAR_ARG2 (&job, 1), + nr, 1, DL, D, DU, DU2, pipvt, + Bz, b_nr, err + F77_CHAR_ARG_LEN (1))); + + if (err != 0) + { + (*current_liboctave_error_handler) + ("SparseMatrix::solve solve failed"); + + err = -1; + break; } - retval.maybe_compress (); - } + // Count non-zeros in work vector and adjust + // space in retval if needed + octave_idx_type new_nnz = 0; + for (octave_idx_type i = 0; i < nr; i++) + if (Bx[i] != 0. || Bz[i] != 0.) + new_nnz++; + + if (ii + new_nnz > x_nz) + { + // Resize the sparse matrix + octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; + retval.change_capacity (sz); + x_nz = sz; + } + + for (octave_idx_type i = 0; i < nr; i++) + if (Bx[i] != 0. || Bz[i] != 0.) + { + retval.xridx(ii) = i; + retval.xdata(ii++) = + Complex (Bx[i], Bz[i]); + } + + retval.xcidx(j+1) = ii; + } + + retval.maybe_compress (); } } else if (typ != MatrixType::Tridiagonal_Hermitian) @@ -4489,85 +4438,71 @@ nr, n_lower, tmp_data, ldm, err F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dpbtrf"); - else - { - if (err != 0) - { - // Matrix is not positive definite!! Fall through to - // unsymmetric banded solver. - mattype.mark_as_unsymmetric (); - typ = MatrixType::Banded; - rcond = 0.0; - err = 0; - } - else - { - if (calc_cond) + if (err != 0) + { + // Matrix is not positive definite!! Fall through to + // unsymmetric banded solver. + mattype.mark_as_unsymmetric (); + typ = MatrixType::Banded; + rcond = 0.0; + err = 0; + } + else + { + if (calc_cond) + { + Array<double> z (3 * nr); + double *pz = z.fortran_vec (); + Array<octave_idx_type> iz (nr); + octave_idx_type *piz = iz.fortran_vec (); + + F77_XFCN (dpbcon, DGBCON, + (F77_CONST_CHAR_ARG2 (&job, 1), + nr, n_lower, tmp_data, ldm, + anorm, rcond, pz, piz, err + F77_CHAR_ARG_LEN (1))); + + if (err != 0) + err = -2; + + volatile double rcond_plus_one = rcond + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcond)) { - Array<double> z (3 * nr); - double *pz = z.fortran_vec (); - Array<octave_idx_type> iz (nr); - octave_idx_type *piz = iz.fortran_vec (); - - F77_XFCN (dpbcon, DGBCON, - (F77_CONST_CHAR_ARG2 (&job, 1), - nr, n_lower, tmp_data, ldm, - anorm, rcond, pz, piz, err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dpbcon"); - - if (err != 0) - err = -2; - - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) - { - err = -2; - - if (sing_handler) - { - sing_handler (rcond); - mattype.mark_as_rectangular (); - } - else - (*current_liboctave_error_handler) - ("matrix singular to machine precision, rcond = %g", - rcond); - } + err = -2; + + if (sing_handler) + { + sing_handler (rcond); + mattype.mark_as_rectangular (); + } + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision, rcond = %g", + rcond); } - else - rcond = 1.; - - if (err == 0) + } + else + rcond = 1.; + + if (err == 0) + { + retval = b; + double *result = retval.fortran_vec (); + + octave_idx_type b_nc = b.cols (); + + F77_XFCN (dpbtrs, DPBTRS, + (F77_CONST_CHAR_ARG2 (&job, 1), + nr, n_lower, b_nc, tmp_data, + ldm, result, b.rows(), err + F77_CHAR_ARG_LEN (1))); + + if (err != 0) { - retval = b; - double *result = retval.fortran_vec (); - - octave_idx_type b_nc = b.cols (); - - F77_XFCN (dpbtrs, DPBTRS, - (F77_CONST_CHAR_ARG2 (&job, 1), - nr, n_lower, b_nc, tmp_data, - ldm, result, b.rows(), err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dpbtrs"); - - if (err != 0) - { - (*current_liboctave_error_handler) - ("SparseMatrix::solve solve failed"); - err = -1; - } + (*current_liboctave_error_handler) + ("SparseMatrix::solve solve failed"); + err = -1; } } } @@ -4616,89 +4551,75 @@ F77_XFCN (dgbtrf, DGBTRF, (nr, nr, n_lower, n_upper, tmp_data, ldm, pipvt, err)); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dgbtrf"); - else - { - // Throw-away extra info LAPACK gives so as to not - // change output. - if (err != 0) - { - err = -2; - rcond = 0.0; - - if (sing_handler) - { - sing_handler (rcond); - mattype.mark_as_rectangular (); - } - else - (*current_liboctave_error_handler) - ("matrix singular to machine precision"); - - } - else - { - if (calc_cond) + // Throw-away extra info LAPACK gives so as to not + // change output. + if (err != 0) + { + err = -2; + rcond = 0.0; + + if (sing_handler) + { + sing_handler (rcond); + mattype.mark_as_rectangular (); + } + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision"); + + } + else + { + if (calc_cond) + { + char job = '1'; + Array<double> z (3 * nr); + double *pz = z.fortran_vec (); + Array<octave_idx_type> iz (nr); + octave_idx_type *piz = iz.fortran_vec (); + + F77_XFCN (dgbcon, DGBCON, + (F77_CONST_CHAR_ARG2 (&job, 1), + nc, n_lower, n_upper, tmp_data, ldm, pipvt, + anorm, rcond, pz, piz, err + F77_CHAR_ARG_LEN (1))); + + if (err != 0) + err = -2; + + volatile double rcond_plus_one = rcond + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcond)) { - char job = '1'; - Array<double> z (3 * nr); - double *pz = z.fortran_vec (); - Array<octave_idx_type> iz (nr); - octave_idx_type *piz = iz.fortran_vec (); - - F77_XFCN (dgbcon, DGBCON, - (F77_CONST_CHAR_ARG2 (&job, 1), - nc, n_lower, n_upper, tmp_data, ldm, pipvt, - anorm, rcond, pz, piz, err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dgbcon"); - - if (err != 0) - err = -2; - - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) - { - err = -2; - - if (sing_handler) - { - sing_handler (rcond); - mattype.mark_as_rectangular (); - } - else - (*current_liboctave_error_handler) - ("matrix singular to machine precision, rcond = %g", - rcond); - } + err = -2; + + if (sing_handler) + { + sing_handler (rcond); + mattype.mark_as_rectangular (); + } + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision, rcond = %g", + rcond); } - else - rcond = 1.; - - if (err == 0) - { - retval = b; - double *result = retval.fortran_vec (); - - octave_idx_type b_nc = b.cols (); - - char job = 'N'; - F77_XFCN (dgbtrs, DGBTRS, - (F77_CONST_CHAR_ARG2 (&job, 1), - nr, n_lower, n_upper, b_nc, tmp_data, - ldm, pipvt, result, b.rows(), err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dgbtrs"); - } + } + else + rcond = 1.; + + if (err == 0) + { + retval = b; + double *result = retval.fortran_vec (); + + octave_idx_type b_nc = b.cols (); + + char job = 'N'; + F77_XFCN (dgbtrs, DGBTRS, + (F77_CONST_CHAR_ARG2 (&job, 1), + nr, n_lower, n_upper, b_nc, tmp_data, + ldm, pipvt, result, b.rows(), err + F77_CHAR_ARG_LEN (1))); } } } @@ -4767,123 +4688,105 @@ nr, n_lower, tmp_data, ldm, err F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dpbtrf"); - else - { - if (err != 0) - { - mattype.mark_as_unsymmetric (); - typ = MatrixType::Banded; - rcond = 0.0; - err = 0; - } - else - { - if (calc_cond) - { - Array<double> z (3 * nr); - double *pz = z.fortran_vec (); - Array<octave_idx_type> iz (nr); - octave_idx_type *piz = iz.fortran_vec (); - - F77_XFCN (dpbcon, DGBCON, - (F77_CONST_CHAR_ARG2 (&job, 1), - nr, n_lower, tmp_data, ldm, - anorm, rcond, pz, piz, err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dpbcon"); - - if (err != 0) - err = -2; - - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) - { - err = -2; - - if (sing_handler) - { - sing_handler (rcond); - mattype.mark_as_rectangular (); - } - else - (*current_liboctave_error_handler) - ("matrix singular to machine precision, rcond = %g", - rcond); - } - } - else - rcond = 1.; - - if (err == 0) + if (err != 0) + { + mattype.mark_as_unsymmetric (); + typ = MatrixType::Banded; + rcond = 0.0; + err = 0; + } + else + { + if (calc_cond) + { + Array<double> z (3 * nr); + double *pz = z.fortran_vec (); + Array<octave_idx_type> iz (nr); + octave_idx_type *piz = iz.fortran_vec (); + + F77_XFCN (dpbcon, DGBCON, + (F77_CONST_CHAR_ARG2 (&job, 1), + nr, n_lower, tmp_data, ldm, + anorm, rcond, pz, piz, err + F77_CHAR_ARG_LEN (1))); + + if (err != 0) + err = -2; + + volatile double rcond_plus_one = rcond + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcond)) { - octave_idx_type b_nr = b.rows (); - octave_idx_type b_nc = b.cols (); - OCTAVE_LOCAL_BUFFER (double, Bx, b_nr); - - // Take a first guess that the number of non-zero terms - // will be as many as in b - volatile octave_idx_type x_nz = b.nnz (); - volatile octave_idx_type ii = 0; - retval = SparseMatrix (b_nr, b_nc, x_nz); - - retval.xcidx(0) = 0; - for (volatile octave_idx_type j = 0; j < b_nc; j++) + err = -2; + + if (sing_handler) { - for (octave_idx_type i = 0; i < b_nr; i++) - Bx[i] = b.elem (i, j); - - F77_XFCN (dpbtrs, DPBTRS, - (F77_CONST_CHAR_ARG2 (&job, 1), - nr, n_lower, 1, tmp_data, - ldm, Bx, b_nr, err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) + sing_handler (rcond); + mattype.mark_as_rectangular (); + } + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision, rcond = %g", + rcond); + } + } + else + rcond = 1.; + + if (err == 0) + { + octave_idx_type b_nr = b.rows (); + octave_idx_type b_nc = b.cols (); + OCTAVE_LOCAL_BUFFER (double, Bx, b_nr); + + // Take a first guess that the number of non-zero terms + // will be as many as in b + volatile octave_idx_type x_nz = b.nnz (); + volatile octave_idx_type ii = 0; + retval = SparseMatrix (b_nr, b_nc, x_nz); + + retval.xcidx(0) = 0; + for (volatile octave_idx_type j = 0; j < b_nc; j++) + { + for (octave_idx_type i = 0; i < b_nr; i++) + Bx[i] = b.elem (i, j); + + F77_XFCN (dpbtrs, DPBTRS, + (F77_CONST_CHAR_ARG2 (&job, 1), + nr, n_lower, 1, tmp_data, + ldm, Bx, b_nr, err + F77_CHAR_ARG_LEN (1))); + + if (err != 0) + { + (*current_liboctave_error_handler) + ("SparseMatrix::solve solve failed"); + err = -1; + break; + } + + for (octave_idx_type i = 0; i < b_nr; i++) + { + double tmp = Bx[i]; + if (tmp != 0.0) { - (*current_liboctave_error_handler) - ("unrecoverable error in dpbtrs"); - err = -1; - break; + if (ii == x_nz) + { + // Resize the sparse matrix + octave_idx_type sz = x_nz * + (b_nc - j) / b_nc; + sz = (sz > 10 ? sz : 10) + x_nz; + retval.change_capacity (sz); + x_nz = sz; + } + retval.xdata(ii) = tmp; + retval.xridx(ii++) = i; } - - if (err != 0) - { - (*current_liboctave_error_handler) - ("SparseMatrix::solve solve failed"); - err = -1; - break; - } - - for (octave_idx_type i = 0; i < b_nr; i++) - { - double tmp = Bx[i]; - if (tmp != 0.0) - { - if (ii == x_nz) - { - // Resize the sparse matrix - octave_idx_type sz = x_nz * - (b_nc - j) / b_nc; - sz = (sz > 10 ? sz : 10) + x_nz; - retval.change_capacity (sz); - x_nz = sz; - } - retval.xdata(ii) = tmp; - retval.xridx(ii++) = i; - } - } - retval.xcidx(j+1) = ii; } - - retval.maybe_compress (); + retval.xcidx(j+1) = ii; } + + retval.maybe_compress (); } } } @@ -4931,127 +4834,110 @@ F77_XFCN (dgbtrf, DGBTRF, (nr, nr, n_lower, n_upper, tmp_data, ldm, pipvt, err)); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dgbtrf"); - else - { - if (err != 0) - { - err = -2; - rcond = 0.0; - - if (sing_handler) - { - sing_handler (rcond); - mattype.mark_as_rectangular (); - } - else - (*current_liboctave_error_handler) - ("matrix singular to machine precision"); - - } - else - { - if (calc_cond) + if (err != 0) + { + err = -2; + rcond = 0.0; + + if (sing_handler) + { + sing_handler (rcond); + mattype.mark_as_rectangular (); + } + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision"); + + } + else + { + if (calc_cond) + { + char job = '1'; + Array<double> z (3 * nr); + double *pz = z.fortran_vec (); + Array<octave_idx_type> iz (nr); + octave_idx_type *piz = iz.fortran_vec (); + + F77_XFCN (dgbcon, DGBCON, + (F77_CONST_CHAR_ARG2 (&job, 1), + nc, n_lower, n_upper, tmp_data, ldm, pipvt, + anorm, rcond, pz, piz, err + F77_CHAR_ARG_LEN (1))); + + if (err != 0) + err = -2; + + volatile double rcond_plus_one = rcond + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcond)) { - char job = '1'; - Array<double> z (3 * nr); - double *pz = z.fortran_vec (); - Array<octave_idx_type> iz (nr); - octave_idx_type *piz = iz.fortran_vec (); - - F77_XFCN (dgbcon, DGBCON, - (F77_CONST_CHAR_ARG2 (&job, 1), - nc, n_lower, n_upper, tmp_data, ldm, pipvt, - anorm, rcond, pz, piz, err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dgbcon"); - - if (err != 0) - err = -2; - - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) - { - err = -2; - - if (sing_handler) - { - sing_handler (rcond); - mattype.mark_as_rectangular (); - } - else - (*current_liboctave_error_handler) - ("matrix singular to machine precision, rcond = %g", - rcond); - } + err = -2; + + if (sing_handler) + { + sing_handler (rcond); + mattype.mark_as_rectangular (); + } + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision, rcond = %g", + rcond); } - else - rcond = 1.; - - if (err == 0) + } + else + rcond = 1.; + + if (err == 0) + { + char job = 'N'; + volatile octave_idx_type x_nz = b.nnz (); + octave_idx_type b_nc = b.cols (); + retval = SparseMatrix (nr, b_nc, x_nz); + retval.xcidx(0) = 0; + volatile octave_idx_type ii = 0; + + OCTAVE_LOCAL_BUFFER (double, work, nr); + + for (volatile octave_idx_type j = 0; j < b_nc; j++) { - char job = 'N'; - volatile octave_idx_type x_nz = b.nnz (); - octave_idx_type b_nc = b.cols (); - retval = SparseMatrix (nr, b_nc, x_nz); - retval.xcidx(0) = 0; - volatile octave_idx_type ii = 0; - - OCTAVE_LOCAL_BUFFER (double, work, nr); - - for (volatile octave_idx_type j = 0; j < b_nc; j++) + for (octave_idx_type i = 0; i < nr; i++) + work[i] = 0.; + for (octave_idx_type i = b.cidx(j); + i < b.cidx(j+1); i++) + work[b.ridx(i)] = b.data(i); + + F77_XFCN (dgbtrs, DGBTRS, + (F77_CONST_CHAR_ARG2 (&job, 1), + nr, n_lower, n_upper, 1, tmp_data, + ldm, pipvt, work, b.rows (), err + F77_CHAR_ARG_LEN (1))); + + // Count non-zeros in work vector and adjust + // space in retval if needed + octave_idx_type new_nnz = 0; + for (octave_idx_type i = 0; i < nr; i++) + if (work[i] != 0.) + new_nnz++; + + if (ii + new_nnz > x_nz) { - for (octave_idx_type i = 0; i < nr; i++) - work[i] = 0.; - for (octave_idx_type i = b.cidx(j); - i < b.cidx(j+1); i++) - work[b.ridx(i)] = b.data(i); - - F77_XFCN (dgbtrs, DGBTRS, - (F77_CONST_CHAR_ARG2 (&job, 1), - nr, n_lower, n_upper, 1, tmp_data, - ldm, pipvt, work, b.rows (), err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("unrecoverable error in dgbtrs"); - break; - } - - // Count non-zeros in work vector and adjust - // space in retval if needed - octave_idx_type new_nnz = 0; - for (octave_idx_type i = 0; i < nr; i++) - if (work[i] != 0.) - new_nnz++; - - if (ii + new_nnz > x_nz) - { - // Resize the sparse matrix - octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; - retval.change_capacity (sz); - x_nz = sz; - } - - for (octave_idx_type i = 0; i < nr; i++) - if (work[i] != 0.) - { - retval.xridx(ii) = i; - retval.xdata(ii++) = work[i]; - } - retval.xcidx(j+1) = ii; + // Resize the sparse matrix + octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; + retval.change_capacity (sz); + x_nz = sz; } - retval.maybe_compress (); + for (octave_idx_type i = 0; i < nr; i++) + if (work[i] != 0.) + { + retval.xridx(ii) = i; + retval.xdata(ii++) = work[i]; + } + retval.xcidx(j+1) = ii; } + + retval.maybe_compress (); } } } @@ -5120,128 +5006,102 @@ nr, n_lower, tmp_data, ldm, err F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dpbtrf"); - else - { - if (err != 0) - { - // Matrix is not positive definite!! Fall through to - // unsymmetric banded solver. - mattype.mark_as_unsymmetric (); - typ = MatrixType::Banded; - rcond = 0.0; - err = 0; - } - else - { - if (calc_cond) - { - Array<double> z (3 * nr); - double *pz = z.fortran_vec (); - Array<octave_idx_type> iz (nr); - octave_idx_type *piz = iz.fortran_vec (); - - F77_XFCN (dpbcon, DGBCON, - (F77_CONST_CHAR_ARG2 (&job, 1), - nr, n_lower, tmp_data, ldm, - anorm, rcond, pz, piz, err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dpbcon"); - - if (err != 0) - err = -2; - - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) - { - err = -2; - - if (sing_handler) - { - sing_handler (rcond); - mattype.mark_as_rectangular (); - } - else - (*current_liboctave_error_handler) - ("matrix singular to machine precision, rcond = %g", - rcond); - } - } - else - rcond = 1.; - - if (err == 0) + if (err != 0) + { + // Matrix is not positive definite!! Fall through to + // unsymmetric banded solver. + mattype.mark_as_unsymmetric (); + typ = MatrixType::Banded; + rcond = 0.0; + err = 0; + } + else + { + if (calc_cond) + { + Array<double> z (3 * nr); + double *pz = z.fortran_vec (); + Array<octave_idx_type> iz (nr); + octave_idx_type *piz = iz.fortran_vec (); + + F77_XFCN (dpbcon, DGBCON, + (F77_CONST_CHAR_ARG2 (&job, 1), + nr, n_lower, tmp_data, ldm, + anorm, rcond, pz, piz, err + F77_CHAR_ARG_LEN (1))); + + if (err != 0) + err = -2; + + volatile double rcond_plus_one = rcond + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcond)) { - octave_idx_type b_nr = b.rows (); - octave_idx_type b_nc = b.cols (); - - OCTAVE_LOCAL_BUFFER (double, Bx, b_nr); - OCTAVE_LOCAL_BUFFER (double, Bz, b_nr); - - retval.resize (b_nr, b_nc); - - for (volatile octave_idx_type j = 0; j < b_nc; j++) + err = -2; + + if (sing_handler) + { + sing_handler (rcond); + mattype.mark_as_rectangular (); + } + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision, rcond = %g", + rcond); + } + } + else + rcond = 1.; + + if (err == 0) + { + octave_idx_type b_nr = b.rows (); + octave_idx_type b_nc = b.cols (); + + OCTAVE_LOCAL_BUFFER (double, Bx, b_nr); + OCTAVE_LOCAL_BUFFER (double, Bz, b_nr); + + retval.resize (b_nr, b_nc); + + for (volatile octave_idx_type j = 0; j < b_nc; j++) + { + for (octave_idx_type i = 0; i < b_nr; i++) { - for (octave_idx_type i = 0; i < b_nr; i++) - { - Complex c = b (i,j); - Bx[i] = std::real (c); - Bz[i] = std::imag (c); - } - - F77_XFCN (dpbtrs, DPBTRS, - (F77_CONST_CHAR_ARG2 (&job, 1), - nr, n_lower, 1, tmp_data, - ldm, Bx, b_nr, err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("unrecoverable error in dpbtrs"); - err = -1; - break; - } - - if (err != 0) - { - (*current_liboctave_error_handler) - ("SparseMatrix::solve solve failed"); - err = -1; - break; - } - - F77_XFCN (dpbtrs, DPBTRS, - (F77_CONST_CHAR_ARG2 (&job, 1), - nr, n_lower, 1, tmp_data, - ldm, Bz, b.rows(), err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("unrecoverable error in dpbtrs"); - err = -1; - break; - } - - if (err != 0) - { - (*current_liboctave_error_handler) - ("SparseMatrix::solve solve failed"); - err = -1; - break; - } - - for (octave_idx_type i = 0; i < b_nr; i++) - retval (i, j) = Complex (Bx[i], Bz[i]); + Complex c = b (i,j); + Bx[i] = std::real (c); + Bz[i] = std::imag (c); + } + + F77_XFCN (dpbtrs, DPBTRS, + (F77_CONST_CHAR_ARG2 (&job, 1), + nr, n_lower, 1, tmp_data, + ldm, Bx, b_nr, err + F77_CHAR_ARG_LEN (1))); + + if (err != 0) + { + (*current_liboctave_error_handler) + ("SparseMatrix::solve solve failed"); + err = -1; + break; } + + F77_XFCN (dpbtrs, DPBTRS, + (F77_CONST_CHAR_ARG2 (&job, 1), + nr, n_lower, 1, tmp_data, + ldm, Bz, b.rows(), err + F77_CHAR_ARG_LEN (1))); + + if (err != 0) + { + (*current_liboctave_error_handler) + ("SparseMatrix::solve solve failed"); + err = -1; + break; + } + + for (octave_idx_type i = 0; i < b_nr; i++) + retval (i, j) = Complex (Bx[i], Bz[i]); } } } @@ -5290,116 +5150,92 @@ F77_XFCN (dgbtrf, DGBTRF, (nr, nr, n_lower, n_upper, tmp_data, ldm, pipvt, err)); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dgbtrf"); - else - { - if (err != 0) - { - err = -2; - rcond = 0.0; - - if (sing_handler) - { - sing_handler (rcond); - mattype.mark_as_rectangular (); - } - else - (*current_liboctave_error_handler) - ("matrix singular to machine precision"); - - } - else - { - if (calc_cond) + if (err != 0) + { + err = -2; + rcond = 0.0; + + if (sing_handler) + { + sing_handler (rcond); + mattype.mark_as_rectangular (); + } + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision"); + + } + else + { + if (calc_cond) + { + char job = '1'; + Array<double> z (3 * nr); + double *pz = z.fortran_vec (); + Array<octave_idx_type> iz (nr); + octave_idx_type *piz = iz.fortran_vec (); + + F77_XFCN (dpbcon, DGBCON, + (F77_CONST_CHAR_ARG2 (&job, 1), + nr, n_lower, tmp_data, ldm, + anorm, rcond, pz, piz, err + F77_CHAR_ARG_LEN (1))); + + if (err != 0) + err = -2; + + volatile double rcond_plus_one = rcond + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcond)) { - char job = '1'; - Array<double> z (3 * nr); - double *pz = z.fortran_vec (); - Array<octave_idx_type> iz (nr); - octave_idx_type *piz = iz.fortran_vec (); - - F77_XFCN (dpbcon, DGBCON, - (F77_CONST_CHAR_ARG2 (&job, 1), - nr, n_lower, tmp_data, ldm, - anorm, rcond, pz, piz, err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dpbcon"); - - if (err != 0) - err = -2; - - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) - { - err = -2; - - if (sing_handler) - { - sing_handler (rcond); - mattype.mark_as_rectangular (); - } - else - (*current_liboctave_error_handler) - ("matrix singular to machine precision, rcond = %g", - rcond); - } + err = -2; + + if (sing_handler) + { + sing_handler (rcond); + mattype.mark_as_rectangular (); + } + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision, rcond = %g", + rcond); } - else - rcond = 1.; - - if (err == 0) + } + else + rcond = 1.; + + if (err == 0) + { + char job = 'N'; + octave_idx_type b_nc = b.cols (); + retval.resize (nr,b_nc); + + OCTAVE_LOCAL_BUFFER (double, Bz, nr); + OCTAVE_LOCAL_BUFFER (double, Bx, nr); + + for (volatile octave_idx_type j = 0; j < b_nc; j++) { - char job = 'N'; - octave_idx_type b_nc = b.cols (); - retval.resize (nr,b_nc); - - OCTAVE_LOCAL_BUFFER (double, Bz, nr); - OCTAVE_LOCAL_BUFFER (double, Bx, nr); - - for (volatile octave_idx_type j = 0; j < b_nc; j++) + for (octave_idx_type i = 0; i < nr; i++) { - for (octave_idx_type i = 0; i < nr; i++) - { - Complex c = b (i, j); - Bx[i] = std::real (c); - Bz[i] = std::imag (c); - } - - F77_XFCN (dgbtrs, DGBTRS, - (F77_CONST_CHAR_ARG2 (&job, 1), - nr, n_lower, n_upper, 1, tmp_data, - ldm, pipvt, Bx, b.rows (), err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("unrecoverable error in dgbtrs"); - break; - } - - F77_XFCN (dgbtrs, DGBTRS, - (F77_CONST_CHAR_ARG2 (&job, 1), - nr, n_lower, n_upper, 1, tmp_data, - ldm, pipvt, Bz, b.rows (), err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("unrecoverable error in dgbtrs"); - break; - } - - for (octave_idx_type i = 0; i < nr; i++) - retval (i, j) = Complex (Bx[i], Bz[i]); + Complex c = b (i, j); + Bx[i] = std::real (c); + Bz[i] = std::imag (c); } + + F77_XFCN (dgbtrs, DGBTRS, + (F77_CONST_CHAR_ARG2 (&job, 1), + nr, n_lower, n_upper, 1, tmp_data, + ldm, pipvt, Bx, b.rows (), err + F77_CHAR_ARG_LEN (1))); + + F77_XFCN (dgbtrs, DGBTRS, + (F77_CONST_CHAR_ARG2 (&job, 1), + nr, n_lower, n_upper, 1, tmp_data, + ldm, pipvt, Bz, b.rows (), err + F77_CHAR_ARG_LEN (1))); + + for (octave_idx_type i = 0; i < nr; i++) + retval (i, j) = Complex (Bx[i], Bz[i]); } } } @@ -5469,160 +5305,134 @@ nr, n_lower, tmp_data, ldm, err F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dpbtrf"); - else - { - if (err != 0) - { - // Matrix is not positive definite!! Fall through to - // unsymmetric banded solver. - mattype.mark_as_unsymmetric (); - typ = MatrixType::Banded; - - rcond = 0.0; - err = 0; - } - else - { - if (calc_cond) + if (err != 0) + { + // Matrix is not positive definite!! Fall through to + // unsymmetric banded solver. + mattype.mark_as_unsymmetric (); + typ = MatrixType::Banded; + + rcond = 0.0; + err = 0; + } + else + { + if (calc_cond) + { + Array<double> z (3 * nr); + double *pz = z.fortran_vec (); + Array<octave_idx_type> iz (nr); + octave_idx_type *piz = iz.fortran_vec (); + + F77_XFCN (dpbcon, DGBCON, + (F77_CONST_CHAR_ARG2 (&job, 1), + nr, n_lower, tmp_data, ldm, + anorm, rcond, pz, piz, err + F77_CHAR_ARG_LEN (1))); + + if (err != 0) + err = -2; + + volatile double rcond_plus_one = rcond + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcond)) { - Array<double> z (3 * nr); - double *pz = z.fortran_vec (); - Array<octave_idx_type> iz (nr); - octave_idx_type *piz = iz.fortran_vec (); - - F77_XFCN (dpbcon, DGBCON, - (F77_CONST_CHAR_ARG2 (&job, 1), - nr, n_lower, tmp_data, ldm, - anorm, rcond, pz, piz, err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dpbcon"); - - if (err != 0) - err = -2; - - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) - { - err = -2; - - if (sing_handler) - { - sing_handler (rcond); - mattype.mark_as_rectangular (); - } - else - (*current_liboctave_error_handler) - ("matrix singular to machine precision, rcond = %g", - rcond); - } - } - else - rcond = 1.; - - if (err == 0) - { - octave_idx_type b_nr = b.rows (); - octave_idx_type b_nc = b.cols (); - OCTAVE_LOCAL_BUFFER (double, Bx, b_nr); - OCTAVE_LOCAL_BUFFER (double, Bz, b_nr); - - // Take a first guess that the number of non-zero terms - // will be as many as in b - volatile octave_idx_type x_nz = b.nnz (); - volatile octave_idx_type ii = 0; - retval = SparseComplexMatrix (b_nr, b_nc, x_nz); - - retval.xcidx(0) = 0; - for (volatile octave_idx_type j = 0; j < b_nc; j++) + err = -2; + + if (sing_handler) { - - for (octave_idx_type i = 0; i < b_nr; i++) - { - Complex c = b (i,j); - Bx[i] = std::real (c); - Bz[i] = std::imag (c); - } - - F77_XFCN (dpbtrs, DPBTRS, - (F77_CONST_CHAR_ARG2 (&job, 1), - nr, n_lower, 1, tmp_data, - ldm, Bx, b_nr, err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("unrecoverable error in dpbtrs"); - err = -1; - break; - } - - if (err != 0) - { - (*current_liboctave_error_handler) - ("SparseMatrix::solve solve failed"); - err = -1; - break; - } - - F77_XFCN (dpbtrs, DPBTRS, - (F77_CONST_CHAR_ARG2 (&job, 1), - nr, n_lower, 1, tmp_data, - ldm, Bz, b_nr, err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("unrecoverable error in dpbtrs"); - err = -1; - break; - } - - if (err != 0) - { - (*current_liboctave_error_handler) - ("SparseMatrix::solve solve failed"); - - err = -1; - break; - } - - // Count non-zeros in work vector and adjust - // space in retval if needed - octave_idx_type new_nnz = 0; - for (octave_idx_type i = 0; i < nr; i++) - if (Bx[i] != 0. || Bz[i] != 0.) - new_nnz++; - - if (ii + new_nnz > x_nz) - { - // Resize the sparse matrix - octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; - retval.change_capacity (sz); - x_nz = sz; - } - - for (octave_idx_type i = 0; i < nr; i++) - if (Bx[i] != 0. || Bz[i] != 0.) - { - retval.xridx(ii) = i; - retval.xdata(ii++) = - Complex (Bx[i], Bz[i]); - } - - retval.xcidx(j+1) = ii; + sing_handler (rcond); + mattype.mark_as_rectangular (); + } + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision, rcond = %g", + rcond); + } + } + else + rcond = 1.; + + if (err == 0) + { + octave_idx_type b_nr = b.rows (); + octave_idx_type b_nc = b.cols (); + OCTAVE_LOCAL_BUFFER (double, Bx, b_nr); + OCTAVE_LOCAL_BUFFER (double, Bz, b_nr); + + // Take a first guess that the number of non-zero terms + // will be as many as in b + volatile octave_idx_type x_nz = b.nnz (); + volatile octave_idx_type ii = 0; + retval = SparseComplexMatrix (b_nr, b_nc, x_nz); + + retval.xcidx(0) = 0; + for (volatile octave_idx_type j = 0; j < b_nc; j++) + { + + for (octave_idx_type i = 0; i < b_nr; i++) + { + Complex c = b (i,j); + Bx[i] = std::real (c); + Bz[i] = std::imag (c); } - retval.maybe_compress (); + F77_XFCN (dpbtrs, DPBTRS, + (F77_CONST_CHAR_ARG2 (&job, 1), + nr, n_lower, 1, tmp_data, + ldm, Bx, b_nr, err + F77_CHAR_ARG_LEN (1))); + + if (err != 0) + { + (*current_liboctave_error_handler) + ("SparseMatrix::solve solve failed"); + err = -1; + break; + } + + F77_XFCN (dpbtrs, DPBTRS, + (F77_CONST_CHAR_ARG2 (&job, 1), + nr, n_lower, 1, tmp_data, + ldm, Bz, b_nr, err + F77_CHAR_ARG_LEN (1))); + + if (err != 0) + { + (*current_liboctave_error_handler) + ("SparseMatrix::solve solve failed"); + + err = -1; + break; + } + + // Count non-zeros in work vector and adjust + // space in retval if needed + octave_idx_type new_nnz = 0; + for (octave_idx_type i = 0; i < nr; i++) + if (Bx[i] != 0. || Bz[i] != 0.) + new_nnz++; + + if (ii + new_nnz > x_nz) + { + // Resize the sparse matrix + octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; + retval.change_capacity (sz); + x_nz = sz; + } + + for (octave_idx_type i = 0; i < nr; i++) + if (Bx[i] != 0. || Bz[i] != 0.) + { + retval.xridx(ii) = i; + retval.xdata(ii++) = + Complex (Bx[i], Bz[i]); + } + + retval.xcidx(j+1) = ii; } + + retval.maybe_compress (); } } } @@ -5670,149 +5480,125 @@ F77_XFCN (dgbtrf, DGBTRF, (nr, nr, n_lower, n_upper, tmp_data, ldm, pipvt, err)); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dgbtrf"); - else - { - if (err != 0) - { - err = -2; - rcond = 0.0; - - if (sing_handler) - { - sing_handler (rcond); - mattype.mark_as_rectangular (); - } - else - (*current_liboctave_error_handler) - ("matrix singular to machine precision"); - - } - else - { - if (calc_cond) - { - char job = '1'; - Array<double> z (3 * nr); - double *pz = z.fortran_vec (); - Array<octave_idx_type> iz (nr); - octave_idx_type *piz = iz.fortran_vec (); - - F77_XFCN (dgbcon, DGBCON, - (F77_CONST_CHAR_ARG2 (&job, 1), - nc, n_lower, n_upper, tmp_data, ldm, pipvt, - anorm, rcond, pz, piz, err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dgbcon"); - - if (err != 0) - err = -2; - - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) - { - err = -2; - - if (sing_handler) - { - sing_handler (rcond); - mattype.mark_as_rectangular (); - } - else - (*current_liboctave_error_handler) - ("matrix singular to machine precision, rcond = %g", - rcond); - } - } - else - rcond = 1.; - - if (err == 0) + if (err != 0) + { + err = -2; + rcond = 0.0; + + if (sing_handler) + { + sing_handler (rcond); + mattype.mark_as_rectangular (); + } + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision"); + + } + else + { + if (calc_cond) + { + char job = '1'; + Array<double> z (3 * nr); + double *pz = z.fortran_vec (); + Array<octave_idx_type> iz (nr); + octave_idx_type *piz = iz.fortran_vec (); + + F77_XFCN (dgbcon, DGBCON, + (F77_CONST_CHAR_ARG2 (&job, 1), + nc, n_lower, n_upper, tmp_data, ldm, pipvt, + anorm, rcond, pz, piz, err + F77_CHAR_ARG_LEN (1))); + + if (err != 0) + err = -2; + + volatile double rcond_plus_one = rcond + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcond)) { - char job = 'N'; - volatile octave_idx_type x_nz = b.nnz (); - octave_idx_type b_nc = b.cols (); - retval = SparseComplexMatrix (nr, b_nc, x_nz); - retval.xcidx(0) = 0; - volatile octave_idx_type ii = 0; - - OCTAVE_LOCAL_BUFFER (double, Bx, nr); - OCTAVE_LOCAL_BUFFER (double, Bz, nr); - - for (volatile octave_idx_type j = 0; j < b_nc; j++) + err = -2; + + if (sing_handler) + { + sing_handler (rcond); + mattype.mark_as_rectangular (); + } + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision, rcond = %g", + rcond); + } + } + else + rcond = 1.; + + if (err == 0) + { + char job = 'N'; + volatile octave_idx_type x_nz = b.nnz (); + octave_idx_type b_nc = b.cols (); + retval = SparseComplexMatrix (nr, b_nc, x_nz); + retval.xcidx(0) = 0; + volatile octave_idx_type ii = 0; + + OCTAVE_LOCAL_BUFFER (double, Bx, nr); + OCTAVE_LOCAL_BUFFER (double, Bz, nr); + + for (volatile octave_idx_type j = 0; j < b_nc; j++) + { + for (octave_idx_type i = 0; i < nr; i++) + { + Bx[i] = 0.; + Bz[i] = 0.; + } + for (octave_idx_type i = b.cidx(j); + i < b.cidx(j+1); i++) { - for (octave_idx_type i = 0; i < nr; i++) - { - Bx[i] = 0.; - Bz[i] = 0.; - } - for (octave_idx_type i = b.cidx(j); - i < b.cidx(j+1); i++) - { - Complex c = b.data(i); - Bx[b.ridx(i)] = std::real (c); - Bz[b.ridx(i)] = std::imag (c); - } - - F77_XFCN (dgbtrs, DGBTRS, - (F77_CONST_CHAR_ARG2 (&job, 1), - nr, n_lower, n_upper, 1, tmp_data, - ldm, pipvt, Bx, b.rows (), err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("unrecoverable error in dgbtrs"); - break; - } - - F77_XFCN (dgbtrs, DGBTRS, - (F77_CONST_CHAR_ARG2 (&job, 1), - nr, n_lower, n_upper, 1, tmp_data, - ldm, pipvt, Bz, b.rows (), err - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("unrecoverable error in dgbtrs"); - break; - } - - // Count non-zeros in work vector and adjust - // space in retval if needed - octave_idx_type new_nnz = 0; - for (octave_idx_type i = 0; i < nr; i++) - if (Bx[i] != 0. || Bz[i] != 0.) - new_nnz++; - - if (ii + new_nnz > x_nz) - { - // Resize the sparse matrix - octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; - retval.change_capacity (sz); - x_nz = sz; - } - - for (octave_idx_type i = 0; i < nr; i++) - if (Bx[i] != 0. || Bz[i] != 0.) - { - retval.xridx(ii) = i; - retval.xdata(ii++) = - Complex (Bx[i], Bz[i]); - } - retval.xcidx(j+1) = ii; + Complex c = b.data(i); + Bx[b.ridx(i)] = std::real (c); + Bz[b.ridx(i)] = std::imag (c); } - retval.maybe_compress (); + F77_XFCN (dgbtrs, DGBTRS, + (F77_CONST_CHAR_ARG2 (&job, 1), + nr, n_lower, n_upper, 1, tmp_data, + ldm, pipvt, Bx, b.rows (), err + F77_CHAR_ARG_LEN (1))); + + F77_XFCN (dgbtrs, DGBTRS, + (F77_CONST_CHAR_ARG2 (&job, 1), + nr, n_lower, n_upper, 1, tmp_data, + ldm, pipvt, Bz, b.rows (), err + F77_CHAR_ARG_LEN (1))); + + // Count non-zeros in work vector and adjust + // space in retval if needed + octave_idx_type new_nnz = 0; + for (octave_idx_type i = 0; i < nr; i++) + if (Bx[i] != 0. || Bz[i] != 0.) + new_nnz++; + + if (ii + new_nnz > x_nz) + { + // Resize the sparse matrix + octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; + retval.change_capacity (sz); + x_nz = sz; + } + + for (octave_idx_type i = 0; i < nr; i++) + if (Bx[i] != 0. || Bz[i] != 0.) + { + retval.xridx(ii) = i; + retval.xdata(ii++) = + Complex (Bx[i], Bz[i]); + } + retval.xcidx(j+1) = ii; } + + retval.maybe_compress (); } } }
--- a/liboctave/dbleAEPBAL.cc Thu Feb 14 17:14:23 2008 -0500 +++ b/liboctave/dbleAEPBAL.cc Thu Feb 14 21:57:50 2008 -0500 @@ -74,28 +74,20 @@ n, p_balanced_mat, n, ilo, ihi, pscale, info F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in dgebal"); - else - { - balancing_mat = Matrix (n, n, 0.0); - for (octave_idx_type i = 0; i < n; i++) - balancing_mat.elem (i ,i) = 1.0; + balancing_mat = Matrix (n, n, 0.0); + for (octave_idx_type i = 0; i < n; i++) + balancing_mat.elem (i ,i) = 1.0; - double *p_balancing_mat = balancing_mat.fortran_vec (); + double *p_balancing_mat = balancing_mat.fortran_vec (); - char side = 'R'; + char side = 'R'; - F77_XFCN (dgebak, DGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1), - F77_CONST_CHAR_ARG2 (&side, 1), - n, ilo, ihi, pscale, n, - p_balancing_mat, n, info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in dgebak"); - } + F77_XFCN (dgebak, DGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1), + F77_CONST_CHAR_ARG2 (&side, 1), + n, ilo, ihi, pscale, n, + p_balancing_mat, n, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); return info; }
--- a/liboctave/dbleCHOL.cc Thu Feb 14 17:14:23 2008 -0500 +++ b/liboctave/dbleCHOL.cc Thu Feb 14 21:57:50 2008 -0500 @@ -76,43 +76,34 @@ n, h, n, info F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in dpotrf"); + xrcond = 0.0; + if (info != 0) + info = -1; + else if (calc_cond) + { + octave_idx_type dpocon_info = 0; + + // Now calculate the condition number for non-singular matrix. + Array<double> z (3*n); + double *pz = z.fortran_vec (); + Array<octave_idx_type> iz (n); + octave_idx_type *piz = iz.fortran_vec (); + F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, + n, anorm, xrcond, pz, piz, dpocon_info + F77_CHAR_ARG_LEN (1))); + + if (dpocon_info != 0) + info = -1; + } else { - xrcond = 0.0; - if (info != 0) - info = -1; - else if (calc_cond) - { - octave_idx_type dpocon_info = 0; - - // Now calculate the condition number for non-singular matrix. - Array<double> z (3*n); - double *pz = z.fortran_vec (); - Array<octave_idx_type> iz (n); - octave_idx_type *piz = iz.fortran_vec (); - F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, - n, anorm, xrcond, pz, piz, dpocon_info - F77_CHAR_ARG_LEN (1))); + // If someone thinks of a more graceful way of doing this (or + // faster for that matter :-)), please let me know! - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dpocon"); - - if (dpocon_info != 0) - info = -1; - } - else - { - // If someone thinks of a more graceful way of doing this (or - // faster for that matter :-)), please let me know! - - if (n > 1) - for (octave_idx_type j = 0; j < a_nc; j++) - for (octave_idx_type i = j+1; i < a_nr; i++) - chol_mat.xelem (i, j) = 0.0; - } + if (n > 1) + for (octave_idx_type j = 0; j < a_nc; j++) + for (octave_idx_type i = j+1; i < a_nr; i++) + chol_mat.xelem (i, j) = 0.0; } return info; @@ -140,21 +131,15 @@ v, n, info F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dpotri"); - else - { - // If someone thinks of a more graceful way of doing this (or - // faster for that matter :-)), please let me know! + // If someone thinks of a more graceful way of doing this (or + // faster for that matter :-)), please let me know! - if (n > 1) - for (octave_idx_type j = 0; j < r_nc; j++) - for (octave_idx_type i = j+1; i < r_nr; i++) - tmp.xelem (i, j) = tmp.xelem (j, i); + if (n > 1) + for (octave_idx_type j = 0; j < r_nc; j++) + for (octave_idx_type i = j+1; i < r_nr; i++) + tmp.xelem (i, j) = tmp.xelem (j, i); - retval = tmp; - } + retval = tmp; } } else
--- a/liboctave/dbleHESS.cc Thu Feb 14 17:14:23 2008 -0500 +++ b/liboctave/dbleHESS.cc Thu Feb 14 21:57:50 2008 -0500 @@ -87,58 +87,36 @@ n, h, n, ilo, ihi, pscale, info F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in dgebal"); - else - { - Array<double> tau (n-1); - double *ptau = tau.fortran_vec (); + Array<double> tau (n-1); + double *ptau = tau.fortran_vec (); - Array<double> work (lwork); - double *pwork = work.fortran_vec (); + Array<double> work (lwork); + double *pwork = work.fortran_vec (); - F77_XFCN (dgehrd, DGEHRD, (n, ilo, ihi, h, n, ptau, pwork, - lwork, info)); + F77_XFCN (dgehrd, DGEHRD, (n, ilo, ihi, h, n, ptau, pwork, + lwork, info)); - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in dgehrd"); - else - { - unitary_hess_mat = hess_mat; - double *z = unitary_hess_mat.fortran_vec (); + unitary_hess_mat = hess_mat; + double *z = unitary_hess_mat.fortran_vec (); - F77_XFCN (dorghr, DORGHR, (n, ilo, ihi, z, n, ptau, pwork, - lwork, info)); + F77_XFCN (dorghr, DORGHR, (n, ilo, ihi, z, n, ptau, pwork, + lwork, info)); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dorghr"); - else - { - F77_XFCN (dgebak, DGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1), - F77_CONST_CHAR_ARG2 (&side, 1), - n, ilo, ihi, pscale, n, z, - n, info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); + F77_XFCN (dgebak, DGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1), + F77_CONST_CHAR_ARG2 (&side, 1), + n, ilo, ihi, pscale, n, z, + n, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dgebak"); - else - { - // If someone thinks of a more graceful way of doing - // this (or faster for that matter :-)), please let - // me know! + // If someone thinks of a more graceful way of doing + // this (or faster for that matter :-)), please let + // me know! - if (n > 2) - for (octave_idx_type j = 0; j < a_nc; j++) - for (octave_idx_type i = j+2; i < a_nr; i++) - hess_mat.elem (i, j) = 0; - } - } - } - } + if (n > 2) + for (octave_idx_type j = 0; j < a_nc; j++) + for (octave_idx_type i = j+2; i < a_nr; i++) + hess_mat.elem (i, j) = 0; return info; }
--- a/liboctave/dbleLU.cc Thu Feb 14 17:14:23 2008 -0500 +++ b/liboctave/dbleLU.cc Thu Feb 14 21:57:50 2008 -0500 @@ -61,10 +61,7 @@ F77_XFCN (dgetrf, DGETRF, (a_nr, a_nc, tmp_data, a_nr, pipvt, info)); - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in dgetrf"); - else - ipvt -= static_cast<octave_idx_type> (1); + ipvt -= static_cast<octave_idx_type> (1); } /*
--- a/liboctave/dbleQR.cc Thu Feb 14 17:14:23 2008 -0500 +++ b/liboctave/dbleQR.cc Thu Feb 14 21:57:50 2008 -0500 @@ -76,56 +76,45 @@ F77_XFCN (dgeqrf, DGEQRF, (m, n, tmp_data, m, ptau, pwork, lwork, info)); - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in dgeqrf"); + if (qr_type == QR::raw) + { + for (octave_idx_type j = 0; j < min_mn; j++) + { + octave_idx_type limit = j < min_mn - 1 ? j : min_mn - 1; + for (octave_idx_type i = limit + 1; i < m; i++) + A_fact.elem (i, j) *= tau.elem (j); + } + + r = A_fact; + + if (m > n) + r.resize (m, n); + } else { - if (qr_type == QR::raw) - { - for (octave_idx_type j = 0; j < min_mn; j++) - { - octave_idx_type limit = j < min_mn - 1 ? j : min_mn - 1; - for (octave_idx_type i = limit + 1; i < m; i++) - A_fact.elem (i, j) *= tau.elem (j); - } - - r = A_fact; + octave_idx_type n2 = (qr_type == QR::economy) ? min_mn : m; - if (m > n) - r.resize (m, n); - } + if (qr_type == QR::economy && m > n) + r.resize (n, n, 0.0); else - { - octave_idx_type n2 = (qr_type == QR::economy) ? min_mn : m; - - if (qr_type == QR::economy && m > n) - r.resize (n, n, 0.0); - else - r.resize (m, n, 0.0); + r.resize (m, n, 0.0); - for (octave_idx_type j = 0; j < n; j++) - { - octave_idx_type limit = j < min_mn-1 ? j : min_mn-1; - for (octave_idx_type i = 0; i <= limit; i++) - r.elem (i, j) = tmp_data[m*j+i]; - } - - lwork = 32 * n2; - work.resize (lwork); - double *pwork2 = work.fortran_vec (); + for (octave_idx_type j = 0; j < n; j++) + { + octave_idx_type limit = j < min_mn-1 ? j : min_mn-1; + for (octave_idx_type i = 0; i <= limit; i++) + r.elem (i, j) = tmp_data[m*j+i]; + } - F77_XFCN (dorgqr, DORGQR, (m, n2, min_mn, tmp_data, m, ptau, - pwork2, lwork, info)); + lwork = 32 * n2; + work.resize (lwork); + double *pwork2 = work.fortran_vec (); - if (f77_exception_encountered) - (*current_liboctave_error_handler) - ("unrecoverable error in dorgqr"); - else - { - q = A_fact; - q.resize (m, n2); - } - } + F77_XFCN (dorgqr, DORGQR, (m, n2, min_mn, tmp_data, m, ptau, + pwork2, lwork, info)); + + q = A_fact; + q.resize (m, n2); } }
--- a/liboctave/dbleQRP.cc Thu Feb 14 17:14:23 2008 -0500 +++ b/liboctave/dbleQRP.cc Thu Feb 14 21:57:50 2008 -0500 @@ -88,51 +88,41 @@ F77_XFCN (dgeqpf, DGEQPF, (m, n, tmp_data, m, pjpvt, ptau, pwork, info)); - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in dgeqpf"); + // Form Permutation matrix (if economy is requested, return the + // indices only!) + + if (qr_type == QR::economy) + { + p.resize (1, n, 0.0); + for (octave_idx_type j = 0; j < n; j++) + p.elem (0, j) = jpvt.elem (j); + } else { - // Form Permutation matrix (if economy is requested, return the - // indices only!) + p.resize (n, n, 0.0); + for (octave_idx_type j = 0; j < n; j++) + p.elem (jpvt.elem (j) - 1, j) = 1.0; + } - if (qr_type == QR::economy) - { - p.resize (1, n, 0.0); - for (octave_idx_type j = 0; j < n; j++) - p.elem (0, j) = jpvt.elem (j); - } - else - { - p.resize (n, n, 0.0); - for (octave_idx_type j = 0; j < n; j++) - p.elem (jpvt.elem (j) - 1, j) = 1.0; - } + octave_idx_type n2 = (qr_type == QR::economy) ? min_mn : m; - octave_idx_type n2 = (qr_type == QR::economy) ? min_mn : m; + if (qr_type == QR::economy && m > n) + r.resize (n, n, 0.0); + else + r.resize (m, n, 0.0); - if (qr_type == QR::economy && m > n) - r.resize (n, n, 0.0); - else - r.resize (m, n, 0.0); - - for (octave_idx_type j = 0; j < n; j++) - { - octave_idx_type limit = j < min_mn-1 ? j : min_mn-1; - for (octave_idx_type i = 0; i <= limit; i++) - r.elem (i, j) = A_fact.elem (i, j); - } + for (octave_idx_type j = 0; j < n; j++) + { + octave_idx_type limit = j < min_mn-1 ? j : min_mn-1; + for (octave_idx_type i = 0; i <= limit; i++) + r.elem (i, j) = A_fact.elem (i, j); + } - F77_XFCN (dorgqr, DORGQR, (m, n2, min_mn, tmp_data, m, ptau, - pwork, lwork, info)); + F77_XFCN (dorgqr, DORGQR, (m, n2, min_mn, tmp_data, m, ptau, + pwork, lwork, info)); - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in dorgqr"); - else - { - q = A_fact; - q.resize (m, n2); - } - } + q = A_fact; + q.resize (m, n2); } /*
--- a/liboctave/dbleSCHUR.cc Thu Feb 14 17:14:23 2008 -0500 +++ b/liboctave/dbleSCHUR.cc Thu Feb 14 21:57:50 2008 -0500 @@ -137,9 +137,6 @@ F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in dgeesx"); - return info; }
--- a/liboctave/dbleSVD.cc Thu Feb 14 17:14:23 2008 -0500 +++ b/liboctave/dbleSVD.cc Thu Feb 14 21:57:50 2008 -0500 @@ -144,28 +144,18 @@ F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in dgesvd"); - else - { - lwork = static_cast<octave_idx_type> (work(0)); - work.resize (lwork); + lwork = static_cast<octave_idx_type> (work(0)); + work.resize (lwork); - F77_XFCN (dgesvd, DGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), - F77_CONST_CHAR_ARG2 (&jobv, 1), - m, n, tmp_data, m, s_vec, u, m, vt, - nrow_vt, work.fortran_vec (), lwork, info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); + F77_XFCN (dgesvd, DGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), + F77_CONST_CHAR_ARG2 (&jobv, 1), + m, n, tmp_data, m, s_vec, u, m, vt, + nrow_vt, work.fortran_vec (), lwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in dgesvd"); - else - { - if (! (jobv == 'N' || jobv == 'O')) - right_sm = right_sm.transpose (); - } - } + if (! (jobv == 'N' || jobv == 'O')) + right_sm = right_sm.transpose (); return info; }
--- a/src/ChangeLog Thu Feb 14 17:14:23 2008 -0500 +++ b/src/ChangeLog Thu Feb 14 21:57:50 2008 -0500 @@ -1,5 +1,8 @@ 2008-02-14 John W. Eaton <jwe@octave.org> + * DLD-FUNCTIONS/balance.cc, DLD-FUNCTIONS/qz.cc: + Don't check f77_exception_encountered. + * sighandlers.cc (user_abort): If interrupting immediately, set octave_interrupt_state if it is not already set.
--- a/src/DLD-FUNCTIONS/balance.cc Thu Feb 14 17:14:23 2008 -0500 +++ b/src/DLD-FUNCTIONS/balance.cc Thu Feb 14 21:57:50 2008 -0500 @@ -253,12 +253,6 @@ nn, ilo, ihi, lscale.fortran_vec (), rscale.fortran_vec (), work.fortran_vec (), info F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - { - error ("unrecoverable error in balance GEP"); - return retval; - } } else { @@ -270,12 +264,6 @@ nn, ilo, ihi, lscale.fortran_vec (), rscale.fortran_vec (), work.fortran_vec (), info F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - { - error ("unrecoverable error in balance GEP"); - return retval; - } } // Since we just want the balancing matrices, we can use dggbal @@ -300,12 +288,6 @@ F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - { - error ("unrecoverable error in balance GEP(L)"); - return retval; - } - // then right F77_XFCN (dggbak, DGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1), @@ -315,12 +297,6 @@ F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - { - error ("unrecoverable error in balance GEP(R)"); - return retval; - } - switch (nargout) { case 0:
--- a/src/DLD-FUNCTIONS/qz.cc Thu Feb 14 17:14:23 2008 -0500 +++ b/src/DLD-FUNCTIONS/qz.cc Thu Feb 14 21:57:50 2008 -0500 @@ -475,12 +475,6 @@ nn, ilo, ihi, lscale.fortran_vec (), rscale.fortran_vec (), work.fortran_vec (), info F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - { - error ("unrecoverable error in qz (bal)"); - return retval; - } } // Since we just want the balancing matrices, we can use dggbal @@ -501,12 +495,6 @@ if (compq == 'V') std::cout << "qz: balancing done; QQ=" << std::endl << QQ << std::endl; #endif - - if (f77_exception_encountered) - { - error ("unrecoverable error in qz (bal-L)"); - return retval; - } } // then right @@ -524,12 +512,6 @@ if (compz == 'V') std::cout << "qz: balancing done; ZZ=" << std::endl << ZZ << std::endl; #endif - - if (f77_exception_encountered) - { - error ("unrecoverable error in qz (bal-R)"); - return retval; - } } static char qz_job; @@ -603,12 +585,6 @@ F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - { - error ("unrecoverable error in qz (dgghrd)"); - return retval; - } - // check if just computing generalized eigenvalues or if we're // actually computing the decomposition @@ -624,12 +600,6 @@ F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - { - error ("unrecoverable error in qz (dhgeqz)"); - return retval; - } } // order the QZ decomposition? @@ -904,12 +874,6 @@ F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); - if (f77_exception_encountered) - { - error ("unrecoverable error in qz (dtgevc)"); - return retval; - } - // now construct the complex form of VV, WW int jj = 0;