Mercurial > hg > octave-nkf
diff liboctave/fCMatrix.cc @ 10314:07ebe522dac2
untabify liboctave C++ sources
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Thu, 11 Feb 2010 12:23:32 -0500 |
| parents | f7ba6cfe7fb7 |
| children | 12884915a8e4 |
line wrap: on
line diff
--- a/liboctave/fCMatrix.cc +++ b/liboctave/fCMatrix.cc @@ -65,35 +65,35 @@ { F77_RET_T F77_FUNC (xilaenv, XILAENV) (const octave_idx_type&, F77_CONST_CHAR_ARG_DECL, - F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, const octave_idx_type&, - const octave_idx_type&, const octave_idx_type&, - octave_idx_type& - F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL); + F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, const octave_idx_type&, + const octave_idx_type&, const octave_idx_type&, + octave_idx_type& + F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL); F77_RET_T F77_FUNC (cgebal, CGEBAL) (F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, FloatComplex*, const octave_idx_type&, octave_idx_type&, - octave_idx_type&, float*, octave_idx_type& - F77_CHAR_ARG_LEN_DECL); + const octave_idx_type&, FloatComplex*, const octave_idx_type&, octave_idx_type&, + octave_idx_type&, float*, octave_idx_type& + F77_CHAR_ARG_LEN_DECL); F77_RET_T F77_FUNC (sgebak, SGEBAK) (F77_CONST_CHAR_ARG_DECL, - F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, float*, - const octave_idx_type&, float*, const octave_idx_type&, octave_idx_type& - F77_CHAR_ARG_LEN_DECL - F77_CHAR_ARG_LEN_DECL); + F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, float*, + const octave_idx_type&, float*, const octave_idx_type&, octave_idx_type& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); F77_RET_T F77_FUNC (cgemm, CGEMM) (F77_CONST_CHAR_ARG_DECL, - F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, - const FloatComplex&, const FloatComplex*, const octave_idx_type&, - const FloatComplex*, const octave_idx_type&, const FloatComplex&, - FloatComplex*, const octave_idx_type& - F77_CHAR_ARG_LEN_DECL - F77_CHAR_ARG_LEN_DECL); + F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, + const FloatComplex&, const FloatComplex*, const octave_idx_type&, + const FloatComplex*, const octave_idx_type&, const FloatComplex&, + FloatComplex*, const octave_idx_type& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); F77_RET_T F77_FUNC (cgemv, CGEMV) (F77_CONST_CHAR_ARG_DECL, @@ -104,127 +104,127 @@ F77_RET_T F77_FUNC (xcdotu, XCDOTU) (const octave_idx_type&, const FloatComplex*, const octave_idx_type&, - const FloatComplex*, const octave_idx_type&, FloatComplex&); + const FloatComplex*, const octave_idx_type&, FloatComplex&); F77_RET_T F77_FUNC (xcdotc, XCDOTC) (const octave_idx_type&, const FloatComplex*, const octave_idx_type&, - const FloatComplex*, const octave_idx_type&, FloatComplex&); + const FloatComplex*, const octave_idx_type&, FloatComplex&); F77_RET_T F77_FUNC (csyrk, CSYRK) (F77_CONST_CHAR_ARG_DECL, - F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, const octave_idx_type&, - const FloatComplex&, const FloatComplex*, const octave_idx_type&, - const FloatComplex&, FloatComplex*, const octave_idx_type& - F77_CHAR_ARG_LEN_DECL - F77_CHAR_ARG_LEN_DECL); + F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, const octave_idx_type&, + const FloatComplex&, const FloatComplex*, const octave_idx_type&, + const FloatComplex&, FloatComplex*, const octave_idx_type& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); F77_RET_T F77_FUNC (cherk, CHERK) (F77_CONST_CHAR_ARG_DECL, - F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, const octave_idx_type&, - const float&, const FloatComplex*, const octave_idx_type&, - const float&, FloatComplex*, const octave_idx_type& - F77_CHAR_ARG_LEN_DECL - F77_CHAR_ARG_LEN_DECL); + F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, const octave_idx_type&, + const float&, const FloatComplex*, const octave_idx_type&, + const float&, FloatComplex*, const octave_idx_type& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); F77_RET_T F77_FUNC (cgetrf, CGETRF) (const octave_idx_type&, const octave_idx_type&, FloatComplex*, const octave_idx_type&, - octave_idx_type*, octave_idx_type&); + octave_idx_type*, octave_idx_type&); F77_RET_T F77_FUNC (cgetrs, CGETRS) (F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, const octave_idx_type&, FloatComplex*, const octave_idx_type&, - const octave_idx_type*, FloatComplex*, const octave_idx_type&, octave_idx_type& - F77_CHAR_ARG_LEN_DECL); + const octave_idx_type&, const octave_idx_type&, FloatComplex*, const octave_idx_type&, + const octave_idx_type*, FloatComplex*, const octave_idx_type&, octave_idx_type& + F77_CHAR_ARG_LEN_DECL); F77_RET_T F77_FUNC (cgetri, CGETRI) (const octave_idx_type&, FloatComplex*, const octave_idx_type&, const octave_idx_type*, - FloatComplex*, const octave_idx_type&, octave_idx_type&); + FloatComplex*, const octave_idx_type&, octave_idx_type&); F77_RET_T F77_FUNC (cgecon, CGECON) (F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, FloatComplex*, - const octave_idx_type&, const float&, float&, - FloatComplex*, float*, octave_idx_type& - F77_CHAR_ARG_LEN_DECL); + const octave_idx_type&, FloatComplex*, + const octave_idx_type&, const float&, float&, + FloatComplex*, float*, octave_idx_type& + F77_CHAR_ARG_LEN_DECL); F77_RET_T F77_FUNC (cgelsy, CGELSY) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, - FloatComplex*, const octave_idx_type&, FloatComplex*, - const octave_idx_type&, octave_idx_type*, float&, octave_idx_type&, - FloatComplex*, const octave_idx_type&, float*, octave_idx_type&); + FloatComplex*, const octave_idx_type&, FloatComplex*, + const octave_idx_type&, octave_idx_type*, float&, octave_idx_type&, + FloatComplex*, const octave_idx_type&, float*, octave_idx_type&); F77_RET_T F77_FUNC (cgelsd, CGELSD) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, - FloatComplex*, const octave_idx_type&, FloatComplex*, - const octave_idx_type&, float*, float&, octave_idx_type&, - FloatComplex*, const octave_idx_type&, float*, - octave_idx_type*, octave_idx_type&); + FloatComplex*, const octave_idx_type&, FloatComplex*, + const octave_idx_type&, float*, float&, octave_idx_type&, + FloatComplex*, const octave_idx_type&, float*, + octave_idx_type*, octave_idx_type&); F77_RET_T F77_FUNC (cpotrf, CPOTRF) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, - FloatComplex*, const octave_idx_type&, - octave_idx_type& F77_CHAR_ARG_LEN_DECL); + FloatComplex*, const octave_idx_type&, + octave_idx_type& F77_CHAR_ARG_LEN_DECL); F77_RET_T F77_FUNC (cpocon, CPOCON) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, - FloatComplex*, const octave_idx_type&, const float&, - float&, FloatComplex*, float*, - octave_idx_type& F77_CHAR_ARG_LEN_DECL); + FloatComplex*, const octave_idx_type&, const float&, + float&, FloatComplex*, float*, + octave_idx_type& F77_CHAR_ARG_LEN_DECL); F77_RET_T F77_FUNC (cpotrs, CPOTRS) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, - const octave_idx_type&, const FloatComplex*, - const octave_idx_type&, FloatComplex*, - const octave_idx_type&, octave_idx_type& - F77_CHAR_ARG_LEN_DECL); + const octave_idx_type&, const FloatComplex*, + const octave_idx_type&, FloatComplex*, + const octave_idx_type&, octave_idx_type& + F77_CHAR_ARG_LEN_DECL); F77_RET_T F77_FUNC (ctrtri, CTRTRI) (F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, const FloatComplex*, - const octave_idx_type&, octave_idx_type& - F77_CHAR_ARG_LEN_DECL - F77_CHAR_ARG_LEN_DECL); + const octave_idx_type&, const FloatComplex*, + const octave_idx_type&, octave_idx_type& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); F77_RET_T F77_FUNC (ctrcon, CTRCON) (F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, - F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, - const FloatComplex*, const octave_idx_type&, float&, - FloatComplex*, float*, octave_idx_type& - F77_CHAR_ARG_LEN_DECL - F77_CHAR_ARG_LEN_DECL - F77_CHAR_ARG_LEN_DECL); + F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, + const FloatComplex*, const octave_idx_type&, float&, + FloatComplex*, float*, octave_idx_type& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); F77_RET_T F77_FUNC (ctrtrs, CTRTRS) (F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, - F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, - const octave_idx_type&, const FloatComplex*, - const octave_idx_type&, FloatComplex*, - const octave_idx_type&, octave_idx_type& - F77_CHAR_ARG_LEN_DECL - F77_CHAR_ARG_LEN_DECL - F77_CHAR_ARG_LEN_DECL); + F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, + const octave_idx_type&, const FloatComplex*, + const octave_idx_type&, FloatComplex*, + const octave_idx_type&, octave_idx_type& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); F77_RET_T F77_FUNC (clartg, CLARTG) (const FloatComplex&, const FloatComplex&, - float&, FloatComplex&, FloatComplex&); + float&, FloatComplex&, FloatComplex&); F77_RET_T F77_FUNC (ctrsyl, CTRSYL) (F77_CONST_CHAR_ARG_DECL, - F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, - const FloatComplex*, const octave_idx_type&, - const FloatComplex*, const octave_idx_type&, - const FloatComplex*, const octave_idx_type&, float&, octave_idx_type& - F77_CHAR_ARG_LEN_DECL - F77_CHAR_ARG_LEN_DECL); + F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, + const FloatComplex*, const octave_idx_type&, + const FloatComplex*, const octave_idx_type&, + const FloatComplex*, const octave_idx_type&, float&, octave_idx_type& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); F77_RET_T F77_FUNC (xclange, XCLANGE) (F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, const octave_idx_type&, const FloatComplex*, - const octave_idx_type&, float*, float& - F77_CHAR_ARG_LEN_DECL); + const octave_idx_type&, const octave_idx_type&, const FloatComplex*, + const octave_idx_type&, float*, float& + F77_CHAR_ARG_LEN_DECL); } static const FloatComplex FloatComplex_NaN_result (octave_Float_NaN, octave_Float_NaN); @@ -331,9 +331,9 @@ if (is_square () && nr > 0) { for (octave_idx_type i = 0; i < nr; i++) - for (octave_idx_type j = i; j < nc; j++) - if (elem (i, j) != conj (elem (j, i))) - return false; + for (octave_idx_type j = i; j < nc; j++) + if (elem (i, j) != conj (elem (j, i))) + return false; return true; } @@ -360,8 +360,8 @@ make_unique (); for (octave_idx_type j = 0; j < a_nc; j++) - for (octave_idx_type i = 0; i < a_nr; i++) - xelem (r+i, c+j) = a.elem (i, j); + for (octave_idx_type i = 0; i < a_nr; i++) + xelem (r+i, c+j) = a.elem (i, j); } return *this; @@ -383,7 +383,7 @@ make_unique (); for (octave_idx_type i = 0; i < a_len; i++) - xelem (r, c+i) = a.elem (i); + xelem (r, c+i) = a.elem (i); } return *this; @@ -405,7 +405,7 @@ make_unique (); for (octave_idx_type i = 0; i < a_len; i++) - xelem (r+i, c) = a.elem (i); + xelem (r+i, c) = a.elem (i); } return *this; @@ -432,7 +432,7 @@ make_unique (); for (octave_idx_type i = 0; i < a_len; i++) - xelem (r+i, c+i) = a.elem (i, i); + xelem (r+i, c+i) = a.elem (i, i); } return *this; @@ -477,7 +477,7 @@ make_unique (); for (octave_idx_type i = 0; i < a_len; i++) - xelem (r+i, c) = a.elem (i); + xelem (r+i, c) = a.elem (i); } return *this; @@ -504,7 +504,7 @@ make_unique (); for (octave_idx_type i = 0; i < a_len; i++) - xelem (r+i, c+i) = a.elem (i, i); + xelem (r+i, c+i) = a.elem (i, i); } return *this; @@ -521,8 +521,8 @@ make_unique (); for (octave_idx_type j = 0; j < nc; j++) - for (octave_idx_type i = 0; i < nr; i++) - xelem (i, j) = val; + for (octave_idx_type i = 0; i < nr; i++) + xelem (i, j) = val; } return *this; @@ -539,8 +539,8 @@ make_unique (); for (octave_idx_type j = 0; j < nc; j++) - for (octave_idx_type i = 0; i < nr; i++) - xelem (i, j) = val; + for (octave_idx_type i = 0; i < nr; i++) + xelem (i, j) = val; } return *this; @@ -567,8 +567,8 @@ make_unique (); for (octave_idx_type j = c1; j <= c2; j++) - for (octave_idx_type i = r1; i <= r2; i++) - xelem (i, j) = val; + for (octave_idx_type i = r1; i <= r2; i++) + xelem (i, j) = val; } return *this; @@ -595,8 +595,8 @@ make_unique (); for (octave_idx_type j = c1; j <= c2; j++) - for (octave_idx_type i = r1; i <= r2; i++) - xelem (i, j) = val; + for (octave_idx_type i = r1; i <= r2; i++) + xelem (i, j) = val; } return *this; @@ -754,7 +754,7 @@ if (nc != a.cols ()) { (*current_liboctave_error_handler) - ("column dimension mismatch for stack"); + ("column dimension mismatch for stack"); return *this; } @@ -773,7 +773,7 @@ if (nc != a.length ()) { (*current_liboctave_error_handler) - ("column dimension mismatch for stack"); + ("column dimension mismatch for stack"); return *this; } @@ -792,7 +792,7 @@ if (nc != 1) { (*current_liboctave_error_handler) - ("column dimension mismatch for stack"); + ("column dimension mismatch for stack"); return *this; } @@ -811,7 +811,7 @@ if (nc != a.cols ()) { (*current_liboctave_error_handler) - ("column dimension mismatch for stack"); + ("column dimension mismatch for stack"); return *this; } @@ -830,7 +830,7 @@ if (nc != a.cols ()) { (*current_liboctave_error_handler) - ("column dimension mismatch for stack"); + ("column dimension mismatch for stack"); return *this; } @@ -849,7 +849,7 @@ if (nc != a.length ()) { (*current_liboctave_error_handler) - ("column dimension mismatch for stack"); + ("column dimension mismatch for stack"); return *this; } @@ -868,7 +868,7 @@ if (nc != 1) { (*current_liboctave_error_handler) - ("column dimension mismatch for stack"); + ("column dimension mismatch for stack"); return *this; } @@ -887,7 +887,7 @@ if (nc != a.cols ()) { (*current_liboctave_error_handler) - ("column dimension mismatch for stack"); + ("column dimension mismatch for stack"); return *this; } @@ -970,7 +970,7 @@ FloatComplexMatrix FloatComplexMatrix::inverse (octave_idx_type& info, float& rcon, int force, - int calc_cond) const + int calc_cond) const { MatrixType mattype (*this); return inverse (mattype, info, rcon, force, calc_cond); @@ -993,7 +993,7 @@ FloatComplexMatrix FloatComplexMatrix::tinverse (MatrixType &mattype, octave_idx_type& info, - float& rcon, int force, int calc_cond) const + float& rcon, int force, int calc_cond) const { FloatComplexMatrix retval; @@ -1011,38 +1011,38 @@ FloatComplex *tmp_data = retval.fortran_vec (); F77_XFCN (ctrtri, CTRTRI, (F77_CONST_CHAR_ARG2 (&uplo, 1), - F77_CONST_CHAR_ARG2 (&udiag, 1), - nr, tmp_data, nr, info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); + F77_CONST_CHAR_ARG2 (&udiag, 1), + nr, tmp_data, nr, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); // Throw-away extra info LAPACK gives so as to not change output. rcon = 0.0; if (info != 0) - info = -1; + info = -1; else if (calc_cond) - { - octave_idx_type ztrcon_info = 0; - char job = '1'; - - OCTAVE_LOCAL_BUFFER (FloatComplex, cwork, 2*nr); - OCTAVE_LOCAL_BUFFER (float, rwork, nr); - - F77_XFCN (ctrcon, CTRCON, (F77_CONST_CHAR_ARG2 (&job, 1), - F77_CONST_CHAR_ARG2 (&uplo, 1), - F77_CONST_CHAR_ARG2 (&udiag, 1), - nr, tmp_data, nr, rcon, - 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; - } + { + octave_idx_type ztrcon_info = 0; + char job = '1'; + + OCTAVE_LOCAL_BUFFER (FloatComplex, cwork, 2*nr); + OCTAVE_LOCAL_BUFFER (float, rwork, nr); + + F77_XFCN (ctrcon, CTRCON, (F77_CONST_CHAR_ARG2 (&job, 1), + F77_CONST_CHAR_ARG2 (&uplo, 1), + F77_CONST_CHAR_ARG2 (&udiag, 1), + nr, tmp_data, nr, rcon, + 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; + } if (info == -1 && ! force) - retval = *this; // Restore matrix contents. + retval = *this; // Restore matrix contents. } return retval; @@ -1050,7 +1050,7 @@ FloatComplexMatrix FloatComplexMatrix::finverse (MatrixType &mattype, octave_idx_type& info, - float& rcon, int force, int calc_cond) const + float& rcon, int force, int calc_cond) const { FloatComplexMatrix retval; @@ -1073,7 +1073,7 @@ // Query the optimum work array size. F77_XFCN (cgetri, CGETRI, (nc, tmp_data, nr, pipvt, - z.fortran_vec (), lwork, info)); + z.fortran_vec (), lwork, info)); lwork = static_cast<octave_idx_type> (std::real(z(0))); lwork = (lwork < 2 *nc ? 2*nc : lwork); @@ -1085,45 +1085,45 @@ // Calculate the norm of the matrix, for later use. float anorm; if (calc_cond) - anorm = retval.abs().sum().row(static_cast<octave_idx_type>(0)).max(); + anorm = retval.abs().sum().row(static_cast<octave_idx_type>(0)).max(); F77_XFCN (cgetrf, CGETRF, (nc, nc, tmp_data, nr, pipvt, info)); // Throw-away extra info LAPACK gives so as to not change output. rcon = 0.0; if (info != 0) - info = -1; + 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<float> rz (2 * nc); - float *prz = rz.fortran_vec (); - F77_XFCN (cgecon, CGECON, (F77_CONST_CHAR_ARG2 (&job, 1), - nc, tmp_data, nr, anorm, - rcon, pz, prz, zgecon_info - F77_CHAR_ARG_LEN (1))); - - if (zgecon_info != 0) - info = -1; - } + { + // Now calculate the condition number for non-singular matrix. + octave_idx_type zgecon_info = 0; + char job = '1'; + Array<float> rz (2 * nc); + float *prz = rz.fortran_vec (); + F77_XFCN (cgecon, CGECON, (F77_CONST_CHAR_ARG2 (&job, 1), + nc, tmp_data, nr, anorm, + rcon, pz, prz, zgecon_info + F77_CHAR_ARG_LEN (1))); + + if (zgecon_info != 0) + info = -1; + } if (info == -1 && ! force) - retval = *this; // Restore contents. + retval = *this; // Restore contents. else - { - octave_idx_type zgetri_info = 0; - - F77_XFCN (cgetri, CGETRI, (nc, tmp_data, nr, pipvt, - pz, lwork, zgetri_info)); - - if (zgetri_info != 0) - info = -1; - } + { + octave_idx_type zgetri_info = 0; + + F77_XFCN (cgetri, CGETRI, (nc, tmp_data, nr, pipvt, + pz, lwork, zgetri_info)); + + if (zgetri_info != 0) + info = -1; + } if (info != 0) - mattype.mark_as_rectangular(); + mattype.mark_as_rectangular(); } return retval; @@ -1131,7 +1131,7 @@ FloatComplexMatrix FloatComplexMatrix::inverse (MatrixType &mattype, octave_idx_type& info, - float& rcon, int force, int calc_cond) const + float& rcon, int force, int calc_cond) const { int typ = mattype.type (false); FloatComplexMatrix ret; @@ -1144,25 +1144,25 @@ else { if (mattype.is_hermitian ()) - { - FloatComplexCHOL chol (*this, info, calc_cond); - if (info == 0) - { - if (calc_cond) - rcon = chol.rcond(); - else - rcon = 1.0; - ret = chol.inverse (); - } - else - mattype.mark_as_unsymmetric (); - } + { + FloatComplexCHOL chol (*this, info, calc_cond); + if (info == 0) + { + if (calc_cond) + rcon = chol.rcond(); + else + rcon = 1.0; + ret = chol.inverse (); + } + else + mattype.mark_as_unsymmetric (); + } if (!mattype.is_hermitian ()) - ret = finverse(mattype, info, rcon, force, calc_cond); + ret = finverse(mattype, info, rcon, force, calc_cond); if ((mattype.is_hermitian () || calc_cond) && rcon == 0.) - ret = FloatComplexMatrix (rows (), columns (), FloatComplex (octave_Float_Inf, 0.)); + ret = FloatComplexMatrix (rows (), columns (), FloatComplex (octave_Float_Inf, 0.)); } return ret; @@ -1188,9 +1188,9 @@ if (tol <= 0.0) { if (nr > nc) - tol = nr * sigma.elem (0) * DBL_EPSILON; + tol = nr * sigma.elem (0) * DBL_EPSILON; else - tol = nc * sigma.elem (0) * DBL_EPSILON; + tol = nc * sigma.elem (0) * DBL_EPSILON; } while (r >= 0 && sigma.elem (r) < tol) @@ -1451,12 +1451,12 @@ octave_quit (); for (octave_idx_type i = 0; i < npts; i++) - prow[i] = tmp_data[i*nr + j]; + prow[i] = tmp_data[i*nr + j]; F77_FUNC (cfftf, CFFTF) (npts, prow, pwsave); for (octave_idx_type i = 0; i < npts; i++) - tmp_data[i*nr + j] = prow[i]; + tmp_data[i*nr + j] = prow[i]; } return retval; @@ -1520,12 +1520,12 @@ octave_quit (); for (octave_idx_type i = 0; i < npts; i++) - prow[i] = tmp_data[i*nr + j]; + prow[i] = tmp_data[i*nr + j]; F77_FUNC (cfftb, CFFTB) (npts, prow, pwsave); for (octave_idx_type i = 0; i < npts; i++) - tmp_data[i*nr + j] = prow[i] / static_cast<float> (npts); + tmp_data[i*nr + j] = prow[i] / static_cast<float> (npts); } return retval; @@ -1712,145 +1712,145 @@ int typ = mattype.type (); if (typ == MatrixType::Unknown) - typ = mattype.type (*this); + typ = mattype.type (*this); // Only calculate the condition number for LU/Cholesky if (typ == MatrixType::Upper) - { - const FloatComplex *tmp_data = fortran_vec (); - octave_idx_type info = 0; - char norm = '1'; - char uplo = 'U'; - char dia = 'N'; - - Array<FloatComplex> z (2 * nc); - FloatComplex *pz = z.fortran_vec (); - Array<float> rz (nc); - float *prz = rz.fortran_vec (); - - F77_XFCN (ctrcon, CTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1), - F77_CONST_CHAR_ARG2 (&uplo, 1), - F77_CONST_CHAR_ARG2 (&dia, 1), - nr, tmp_data, nr, rcon, - pz, prz, info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); - - if (info != 0) - rcon = 0; - } + { + const FloatComplex *tmp_data = fortran_vec (); + octave_idx_type info = 0; + char norm = '1'; + char uplo = 'U'; + char dia = 'N'; + + Array<FloatComplex> z (2 * nc); + FloatComplex *pz = z.fortran_vec (); + Array<float> rz (nc); + float *prz = rz.fortran_vec (); + + F77_XFCN (ctrcon, CTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1), + F77_CONST_CHAR_ARG2 (&uplo, 1), + F77_CONST_CHAR_ARG2 (&dia, 1), + nr, tmp_data, nr, rcon, + pz, prz, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (info != 0) + rcon = 0; + } else if (typ == MatrixType::Permuted_Upper) - (*current_liboctave_error_handler) - ("permuted triangular matrix not implemented"); + (*current_liboctave_error_handler) + ("permuted triangular matrix not implemented"); else if (typ == MatrixType::Lower) - { - const FloatComplex *tmp_data = fortran_vec (); - octave_idx_type info = 0; - char norm = '1'; - char uplo = 'L'; - char dia = 'N'; - - Array<FloatComplex> z (2 * nc); - FloatComplex *pz = z.fortran_vec (); - Array<float> rz (nc); - float *prz = rz.fortran_vec (); - - F77_XFCN (ctrcon, CTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1), - F77_CONST_CHAR_ARG2 (&uplo, 1), - F77_CONST_CHAR_ARG2 (&dia, 1), - nr, tmp_data, nr, rcon, - pz, prz, info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); - - if (info != 0) - rcon = 0.0; - } + { + const FloatComplex *tmp_data = fortran_vec (); + octave_idx_type info = 0; + char norm = '1'; + char uplo = 'L'; + char dia = 'N'; + + Array<FloatComplex> z (2 * nc); + FloatComplex *pz = z.fortran_vec (); + Array<float> rz (nc); + float *prz = rz.fortran_vec (); + + F77_XFCN (ctrcon, CTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1), + F77_CONST_CHAR_ARG2 (&uplo, 1), + F77_CONST_CHAR_ARG2 (&dia, 1), + nr, tmp_data, nr, rcon, + pz, prz, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (info != 0) + rcon = 0.0; + } else if (typ == MatrixType::Permuted_Lower) - (*current_liboctave_error_handler) - ("permuted triangular matrix not implemented"); + (*current_liboctave_error_handler) + ("permuted triangular matrix not implemented"); else if (typ == MatrixType::Full || typ == MatrixType::Hermitian) - { - float anorm = -1.0; - FloatComplexMatrix atmp = *this; - FloatComplex *tmp_data = atmp.fortran_vec (); - - if (typ == MatrixType::Hermitian) - { - octave_idx_type info = 0; - char job = 'L'; - anorm = atmp.abs().sum(). - row(static_cast<octave_idx_type>(0)).max(); - - F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr, - tmp_data, nr, info - F77_CHAR_ARG_LEN (1))); - - if (info != 0) - { - rcon = 0.0; - - mattype.mark_as_unsymmetric (); - typ = MatrixType::Full; - } - else - { - Array<FloatComplex> z (2 * nc); - FloatComplex *pz = z.fortran_vec (); - Array<float> rz (nc); - float *prz = rz.fortran_vec (); - - F77_XFCN (cpocon, CPOCON, (F77_CONST_CHAR_ARG2 (&job, 1), - nr, tmp_data, nr, anorm, - rcon, pz, prz, info - F77_CHAR_ARG_LEN (1))); - - if (info != 0) - rcon = 0.0; - } - } - - - if (typ == MatrixType::Full) - { - octave_idx_type info = 0; - - Array<octave_idx_type> ipvt (nr); - octave_idx_type *pipvt = ipvt.fortran_vec (); - - if(anorm < 0.) - anorm = atmp.abs().sum(). - row(static_cast<octave_idx_type>(0)).max(); - - Array<FloatComplex> z (2 * nc); - FloatComplex *pz = z.fortran_vec (); - Array<float> rz (2 * nc); - float *prz = rz.fortran_vec (); - - F77_XFCN (cgetrf, CGETRF, (nr, nr, tmp_data, nr, pipvt, info)); - - if (info != 0) - { - rcon = 0.0; - mattype.mark_as_rectangular (); - } - else - { - char job = '1'; - F77_XFCN (cgecon, CGECON, (F77_CONST_CHAR_ARG2 (&job, 1), - nc, tmp_data, nr, anorm, - rcon, pz, prz, info - F77_CHAR_ARG_LEN (1))); - - if (info != 0) - rcon = 0.0; - } - } - } + { + float anorm = -1.0; + FloatComplexMatrix atmp = *this; + FloatComplex *tmp_data = atmp.fortran_vec (); + + if (typ == MatrixType::Hermitian) + { + octave_idx_type info = 0; + char job = 'L'; + anorm = atmp.abs().sum(). + row(static_cast<octave_idx_type>(0)).max(); + + F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr, + tmp_data, nr, info + F77_CHAR_ARG_LEN (1))); + + if (info != 0) + { + rcon = 0.0; + + mattype.mark_as_unsymmetric (); + typ = MatrixType::Full; + } + else + { + Array<FloatComplex> z (2 * nc); + FloatComplex *pz = z.fortran_vec (); + Array<float> rz (nc); + float *prz = rz.fortran_vec (); + + F77_XFCN (cpocon, CPOCON, (F77_CONST_CHAR_ARG2 (&job, 1), + nr, tmp_data, nr, anorm, + rcon, pz, prz, info + F77_CHAR_ARG_LEN (1))); + + if (info != 0) + rcon = 0.0; + } + } + + + if (typ == MatrixType::Full) + { + octave_idx_type info = 0; + + Array<octave_idx_type> ipvt (nr); + octave_idx_type *pipvt = ipvt.fortran_vec (); + + if(anorm < 0.) + anorm = atmp.abs().sum(). + row(static_cast<octave_idx_type>(0)).max(); + + Array<FloatComplex> z (2 * nc); + FloatComplex *pz = z.fortran_vec (); + Array<float> rz (2 * nc); + float *prz = rz.fortran_vec (); + + F77_XFCN (cgetrf, CGETRF, (nr, nr, tmp_data, nr, pipvt, info)); + + if (info != 0) + { + rcon = 0.0; + mattype.mark_as_rectangular (); + } + else + { + char job = '1'; + F77_XFCN (cgecon, CGECON, (F77_CONST_CHAR_ARG2 (&job, 1), + nc, tmp_data, nr, anorm, + rcon, pz, prz, info + F77_CHAR_ARG_LEN (1))); + + if (info != 0) + rcon = 0.0; + } + } + } else - rcon = 0.0; + rcon = 0.0; } return rcon; @@ -1858,9 +1858,9 @@ FloatComplexMatrix FloatComplexMatrix::utsolve (MatrixType &mattype, const FloatComplexMatrix& b, - octave_idx_type& info, float& rcon, - solve_singularity_handler sing_handler, - bool calc_cond, blas_trans_type transt) const + octave_idx_type& info, float& rcon, + solve_singularity_handler sing_handler, + bool calc_cond, blas_trans_type transt) const { FloatComplexMatrix retval; @@ -1877,81 +1877,81 @@ volatile int typ = mattype.type (); if (typ == MatrixType::Permuted_Upper || - typ == MatrixType::Upper) - { - octave_idx_type b_nc = b.cols (); - rcon = 1.; - info = 0; - - if (typ == MatrixType::Permuted_Upper) - { - (*current_liboctave_error_handler) - ("permuted triangular matrix not implemented"); - } - else - { - const FloatComplex *tmp_data = fortran_vec (); - - if (calc_cond) - { - char norm = '1'; - char uplo = 'U'; - char dia = 'N'; - - Array<FloatComplex> z (2 * nc); - FloatComplex *pz = z.fortran_vec (); - Array<float> rz (nc); - float *prz = rz.fortran_vec (); - - F77_XFCN (ctrcon, CTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1), - F77_CONST_CHAR_ARG2 (&uplo, 1), - F77_CONST_CHAR_ARG2 (&dia, 1), - nr, tmp_data, nr, rcon, - pz, prz, info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); - - if (info != 0) - info = -2; - - volatile float rcond_plus_one = rcon + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcon)) - { - info = -2; - - if (sing_handler) - sing_handler (rcon); - else - (*current_liboctave_error_handler) - ("matrix singular to machine precision, rcond = %g", - rcon); - } - } - - if (info == 0) - { - retval = b; - FloatComplex *result = retval.fortran_vec (); - - char uplo = 'U'; - char trans = get_blas_char (transt); - char dia = 'N'; - - F77_XFCN (ctrtrs, CTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1), - F77_CONST_CHAR_ARG2 (&trans, 1), - F77_CONST_CHAR_ARG2 (&dia, 1), - nr, b_nc, tmp_data, nr, - result, nr, info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); - } - } - } + typ == MatrixType::Upper) + { + octave_idx_type b_nc = b.cols (); + rcon = 1.; + info = 0; + + if (typ == MatrixType::Permuted_Upper) + { + (*current_liboctave_error_handler) + ("permuted triangular matrix not implemented"); + } + else + { + const FloatComplex *tmp_data = fortran_vec (); + + if (calc_cond) + { + char norm = '1'; + char uplo = 'U'; + char dia = 'N'; + + Array<FloatComplex> z (2 * nc); + FloatComplex *pz = z.fortran_vec (); + Array<float> rz (nc); + float *prz = rz.fortran_vec (); + + F77_XFCN (ctrcon, CTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1), + F77_CONST_CHAR_ARG2 (&uplo, 1), + F77_CONST_CHAR_ARG2 (&dia, 1), + nr, tmp_data, nr, rcon, + pz, prz, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (info != 0) + info = -2; + + volatile float rcond_plus_one = rcon + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcon)) + { + info = -2; + + if (sing_handler) + sing_handler (rcon); + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision, rcond = %g", + rcon); + } + } + + if (info == 0) + { + retval = b; + FloatComplex *result = retval.fortran_vec (); + + char uplo = 'U'; + char trans = get_blas_char (transt); + char dia = 'N'; + + F77_XFCN (ctrtrs, CTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1), + F77_CONST_CHAR_ARG2 (&trans, 1), + F77_CONST_CHAR_ARG2 (&dia, 1), + nr, b_nc, tmp_data, nr, + result, nr, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + } + } + } else - (*current_liboctave_error_handler) ("incorrect matrix type"); + (*current_liboctave_error_handler) ("incorrect matrix type"); } return retval; @@ -1959,9 +1959,9 @@ FloatComplexMatrix FloatComplexMatrix::ltsolve (MatrixType &mattype, const FloatComplexMatrix& b, - octave_idx_type& info, float& rcon, - solve_singularity_handler sing_handler, - bool calc_cond, blas_trans_type transt) const + octave_idx_type& info, float& rcon, + solve_singularity_handler sing_handler, + bool calc_cond, blas_trans_type transt) const { FloatComplexMatrix retval; @@ -1978,81 +1978,81 @@ volatile int typ = mattype.type (); if (typ == MatrixType::Permuted_Lower || - typ == MatrixType::Lower) - { - octave_idx_type b_nc = b.cols (); - rcon = 1.; - info = 0; - - if (typ == MatrixType::Permuted_Lower) - { - (*current_liboctave_error_handler) - ("permuted triangular matrix not implemented"); - } - else - { - const FloatComplex *tmp_data = fortran_vec (); - - if (calc_cond) - { - char norm = '1'; - char uplo = 'L'; - char dia = 'N'; - - Array<FloatComplex> z (2 * nc); - FloatComplex *pz = z.fortran_vec (); - Array<float> rz (nc); - float *prz = rz.fortran_vec (); - - F77_XFCN (ctrcon, CTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1), - F77_CONST_CHAR_ARG2 (&uplo, 1), - F77_CONST_CHAR_ARG2 (&dia, 1), - nr, tmp_data, nr, rcon, - pz, prz, info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); - - if (info != 0) - info = -2; - - volatile float rcond_plus_one = rcon + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcon)) - { - info = -2; - - if (sing_handler) - sing_handler (rcon); - else - (*current_liboctave_error_handler) - ("matrix singular to machine precision, rcond = %g", - rcon); - } - } - - if (info == 0) - { - retval = b; - FloatComplex *result = retval.fortran_vec (); - - char uplo = 'L'; - char trans = get_blas_char (transt); - char dia = 'N'; - - F77_XFCN (ctrtrs, CTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1), - F77_CONST_CHAR_ARG2 (&trans, 1), - F77_CONST_CHAR_ARG2 (&dia, 1), - nr, b_nc, tmp_data, nr, - result, nr, info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); - } - } - } + typ == MatrixType::Lower) + { + octave_idx_type b_nc = b.cols (); + rcon = 1.; + info = 0; + + if (typ == MatrixType::Permuted_Lower) + { + (*current_liboctave_error_handler) + ("permuted triangular matrix not implemented"); + } + else + { + const FloatComplex *tmp_data = fortran_vec (); + + if (calc_cond) + { + char norm = '1'; + char uplo = 'L'; + char dia = 'N'; + + Array<FloatComplex> z (2 * nc); + FloatComplex *pz = z.fortran_vec (); + Array<float> rz (nc); + float *prz = rz.fortran_vec (); + + F77_XFCN (ctrcon, CTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1), + F77_CONST_CHAR_ARG2 (&uplo, 1), + F77_CONST_CHAR_ARG2 (&dia, 1), + nr, tmp_data, nr, rcon, + pz, prz, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (info != 0) + info = -2; + + volatile float rcond_plus_one = rcon + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcon)) + { + info = -2; + + if (sing_handler) + sing_handler (rcon); + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision, rcond = %g", + rcon); + } + } + + if (info == 0) + { + retval = b; + FloatComplex *result = retval.fortran_vec (); + + char uplo = 'L'; + char trans = get_blas_char (transt); + char dia = 'N'; + + F77_XFCN (ctrtrs, CTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1), + F77_CONST_CHAR_ARG2 (&trans, 1), + F77_CONST_CHAR_ARG2 (&dia, 1), + nr, b_nc, tmp_data, nr, + result, nr, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + } + } + } else - (*current_liboctave_error_handler) ("incorrect matrix type"); + (*current_liboctave_error_handler) ("incorrect matrix type"); } return retval; @@ -2060,9 +2060,9 @@ FloatComplexMatrix FloatComplexMatrix::fsolve (MatrixType &mattype, const FloatComplexMatrix& b, - octave_idx_type& info, float& rcon, - solve_singularity_handler sing_handler, - bool calc_cond) const + octave_idx_type& info, float& rcon, + solve_singularity_handler sing_handler, + bool calc_cond) const { FloatComplexMatrix retval; @@ -2083,160 +2083,160 @@ float anorm = -1.; if (typ == MatrixType::Hermitian) - { - info = 0; - char job = 'L'; - FloatComplexMatrix atmp = *this; - FloatComplex *tmp_data = atmp.fortran_vec (); - anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max(); - - F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr, - tmp_data, nr, info - F77_CHAR_ARG_LEN (1))); - - // Throw-away extra info LAPACK gives so as to not change output. - rcon = 0.0; - if (info != 0) - { - info = -2; - - mattype.mark_as_unsymmetric (); - typ = MatrixType::Full; - } - else - { - if (calc_cond) - { - Array<FloatComplex> z (2 * nc); - FloatComplex *pz = z.fortran_vec (); - Array<float> rz (nc); - float *prz = rz.fortran_vec (); - - F77_XFCN (cpocon, CPOCON, (F77_CONST_CHAR_ARG2 (&job, 1), - nr, tmp_data, nr, anorm, - rcon, pz, prz, info - F77_CHAR_ARG_LEN (1))); - - if (info != 0) - info = -2; - - volatile float rcond_plus_one = rcon + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcon)) - { - info = -2; - - if (sing_handler) - sing_handler (rcon); - else - (*current_liboctave_error_handler) - ("matrix singular to machine precision, rcond = %g", - rcon); - } - } - - if (info == 0) - { - retval = b; - FloatComplex *result = retval.fortran_vec (); - - octave_idx_type b_nc = b.cols (); - - F77_XFCN (cpotrs, CPOTRS, (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; - } - } - } + { + info = 0; + char job = 'L'; + FloatComplexMatrix atmp = *this; + FloatComplex *tmp_data = atmp.fortran_vec (); + anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max(); + + F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr, + tmp_data, nr, info + F77_CHAR_ARG_LEN (1))); + + // Throw-away extra info LAPACK gives so as to not change output. + rcon = 0.0; + if (info != 0) + { + info = -2; + + mattype.mark_as_unsymmetric (); + typ = MatrixType::Full; + } + else + { + if (calc_cond) + { + Array<FloatComplex> z (2 * nc); + FloatComplex *pz = z.fortran_vec (); + Array<float> rz (nc); + float *prz = rz.fortran_vec (); + + F77_XFCN (cpocon, CPOCON, (F77_CONST_CHAR_ARG2 (&job, 1), + nr, tmp_data, nr, anorm, + rcon, pz, prz, info + F77_CHAR_ARG_LEN (1))); + + if (info != 0) + info = -2; + + volatile float rcond_plus_one = rcon + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcon)) + { + info = -2; + + if (sing_handler) + sing_handler (rcon); + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision, rcond = %g", + rcon); + } + } + + if (info == 0) + { + retval = b; + FloatComplex *result = retval.fortran_vec (); + + octave_idx_type b_nc = b.cols (); + + F77_XFCN (cpotrs, CPOTRS, (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; + } + } + } if (typ == MatrixType::Full) - { - info = 0; - - Array<octave_idx_type> ipvt (nr); - octave_idx_type *pipvt = ipvt.fortran_vec (); - - FloatComplexMatrix atmp = *this; - FloatComplex *tmp_data = atmp.fortran_vec (); - - Array<FloatComplex> z (2 * nc); - FloatComplex *pz = z.fortran_vec (); - Array<float> rz (2 * nc); - float *prz = rz.fortran_vec (); - - // Calculate the norm of the matrix, for later use. - if (anorm < 0.) - anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max(); - - F77_XFCN (cgetrf, CGETRF, (nr, nr, tmp_data, nr, pipvt, info)); - - // Throw-away extra info LAPACK gives so as to not change output. - rcon = 0.0; - if (info != 0) - { - info = -2; - - if (sing_handler) - sing_handler (rcon); - 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 (cgecon, CGECON, (F77_CONST_CHAR_ARG2 (&job, 1), - nc, tmp_data, nr, anorm, - rcon, pz, prz, info - F77_CHAR_ARG_LEN (1))); - - if (info != 0) - info = -2; - - volatile float rcond_plus_one = rcon + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcon)) - { - info = -2; - - if (sing_handler) - sing_handler (rcon); - else - (*current_liboctave_error_handler) - ("matrix singular to machine precision, rcond = %g", - rcon); - } - } - - if (info == 0) - { - retval = b; - FloatComplex *result = retval.fortran_vec (); - - octave_idx_type b_nc = b.cols (); - - char job = 'N'; - F77_XFCN (cgetrs, CGETRS, (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 (); - } - } + { + info = 0; + + Array<octave_idx_type> ipvt (nr); + octave_idx_type *pipvt = ipvt.fortran_vec (); + + FloatComplexMatrix atmp = *this; + FloatComplex *tmp_data = atmp.fortran_vec (); + + Array<FloatComplex> z (2 * nc); + FloatComplex *pz = z.fortran_vec (); + Array<float> rz (2 * nc); + float *prz = rz.fortran_vec (); + + // Calculate the norm of the matrix, for later use. + if (anorm < 0.) + anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max(); + + F77_XFCN (cgetrf, CGETRF, (nr, nr, tmp_data, nr, pipvt, info)); + + // Throw-away extra info LAPACK gives so as to not change output. + rcon = 0.0; + if (info != 0) + { + info = -2; + + if (sing_handler) + sing_handler (rcon); + 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 (cgecon, CGECON, (F77_CONST_CHAR_ARG2 (&job, 1), + nc, tmp_data, nr, anorm, + rcon, pz, prz, info + F77_CHAR_ARG_LEN (1))); + + if (info != 0) + info = -2; + + volatile float rcond_plus_one = rcon + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcon)) + { + info = -2; + + if (sing_handler) + sing_handler (rcon); + else + (*current_liboctave_error_handler) + ("matrix singular to machine precision, rcond = %g", + rcon); + } + } + + if (info == 0) + { + retval = b; + FloatComplex *result = retval.fortran_vec (); + + octave_idx_type b_nc = b.cols (); + + char job = 'N'; + F77_XFCN (cgetrs, CGETRS, (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 (); + } + } } return retval; @@ -2252,7 +2252,7 @@ FloatComplexMatrix FloatComplexMatrix::solve (MatrixType &typ, const FloatMatrix& b, - octave_idx_type& info) const + octave_idx_type& info) const { float rcon; return solve (typ, b, info, rcon, 0); @@ -2260,15 +2260,15 @@ FloatComplexMatrix FloatComplexMatrix::solve (MatrixType &typ, const FloatMatrix& b, octave_idx_type& info, - float& rcon) const + float& rcon) const { return solve (typ, b, info, rcon, 0); } FloatComplexMatrix FloatComplexMatrix::solve (MatrixType &typ, const FloatMatrix& b, octave_idx_type& info, - float& rcon, solve_singularity_handler sing_handler, - bool singular_fallback, blas_trans_type transt) const + float& rcon, solve_singularity_handler sing_handler, + bool singular_fallback, blas_trans_type transt) const { FloatComplexMatrix tmp (b); return solve (typ, tmp, info, rcon, sing_handler, singular_fallback, transt); @@ -2284,7 +2284,7 @@ FloatComplexMatrix FloatComplexMatrix::solve (MatrixType &typ, const FloatComplexMatrix& b, - octave_idx_type& info) const + octave_idx_type& info) const { float rcon; return solve (typ, b, info, rcon, 0); @@ -2292,16 +2292,16 @@ FloatComplexMatrix FloatComplexMatrix::solve (MatrixType &typ, const FloatComplexMatrix& b, - octave_idx_type& info, float& rcon) const + octave_idx_type& info, float& rcon) const { return solve (typ, b, info, rcon, 0); } FloatComplexMatrix FloatComplexMatrix::solve (MatrixType &mattype, const FloatComplexMatrix& b, - octave_idx_type& info, float& rcon, - solve_singularity_handler sing_handler, - bool singular_fallback, blas_trans_type transt) const + octave_idx_type& info, float& rcon, + solve_singularity_handler sing_handler, + bool singular_fallback, blas_trans_type transt) const { FloatComplexMatrix retval; int typ = mattype.type (); @@ -2346,7 +2346,7 @@ FloatComplexColumnVector FloatComplexMatrix::solve (MatrixType &typ, const FloatColumnVector& b, - octave_idx_type& info) const + octave_idx_type& info) const { float rcon; return solve (typ, FloatComplexColumnVector (b), info, rcon, 0); @@ -2354,15 +2354,15 @@ FloatComplexColumnVector FloatComplexMatrix::solve (MatrixType &typ, const FloatColumnVector& b, - octave_idx_type& info, float& rcon) const + octave_idx_type& info, float& rcon) const { return solve (typ, FloatComplexColumnVector (b), info, rcon, 0); } FloatComplexColumnVector FloatComplexMatrix::solve (MatrixType &typ, const FloatColumnVector& b, - octave_idx_type& info, float& rcon, - solve_singularity_handler sing_handler, blas_trans_type transt) const + octave_idx_type& info, float& rcon, + solve_singularity_handler sing_handler, blas_trans_type transt) const { return solve (typ, FloatComplexColumnVector (b), info, rcon, sing_handler, transt); } @@ -2377,7 +2377,7 @@ FloatComplexColumnVector FloatComplexMatrix::solve (MatrixType &typ, const FloatComplexColumnVector& b, - octave_idx_type& info) const + octave_idx_type& info) const { float rcon; return solve (typ, b, info, rcon, 0); @@ -2385,15 +2385,15 @@ FloatComplexColumnVector FloatComplexMatrix::solve (MatrixType &typ, const FloatComplexColumnVector& b, - octave_idx_type& info, float& rcon) const + octave_idx_type& info, float& rcon) const { return solve (typ, b, info, rcon, 0); } FloatComplexColumnVector FloatComplexMatrix::solve (MatrixType &typ, const FloatComplexColumnVector& b, - octave_idx_type& info, float& rcon, - solve_singularity_handler sing_handler, blas_trans_type transt) const + octave_idx_type& info, float& rcon, + solve_singularity_handler sing_handler, blas_trans_type transt) const { FloatComplexMatrix tmp (b); @@ -2423,7 +2423,7 @@ FloatComplexMatrix FloatComplexMatrix::solve (const FloatMatrix& b, octave_idx_type& info, float& rcon, - solve_singularity_handler sing_handler, blas_trans_type transt) const + solve_singularity_handler sing_handler, blas_trans_type transt) const { FloatComplexMatrix tmp (b); return solve (tmp, info, rcon, sing_handler, transt); @@ -2452,7 +2452,7 @@ FloatComplexMatrix FloatComplexMatrix::solve (const FloatComplexMatrix& b, octave_idx_type& info, float& rcon, - solve_singularity_handler sing_handler, blas_trans_type transt) const + solve_singularity_handler sing_handler, blas_trans_type transt) const { MatrixType mattype (*this); return solve (mattype, b, info, rcon, sing_handler, true, transt); @@ -2475,15 +2475,15 @@ FloatComplexColumnVector FloatComplexMatrix::solve (const FloatColumnVector& b, octave_idx_type& info, - float& rcon) const + float& rcon) const { return solve (FloatComplexColumnVector (b), info, rcon, 0); } FloatComplexColumnVector FloatComplexMatrix::solve (const FloatColumnVector& b, octave_idx_type& info, - float& rcon, - solve_singularity_handler sing_handler, blas_trans_type transt) const + float& rcon, + solve_singularity_handler sing_handler, blas_trans_type transt) const { return solve (FloatComplexColumnVector (b), info, rcon, sing_handler, transt); } @@ -2505,15 +2505,15 @@ FloatComplexColumnVector FloatComplexMatrix::solve (const FloatComplexColumnVector& b, octave_idx_type& info, - float& rcon) const + float& rcon) const { return solve (b, info, rcon, 0); } FloatComplexColumnVector FloatComplexMatrix::solve (const FloatComplexColumnVector& b, octave_idx_type& info, - float& rcon, - solve_singularity_handler sing_handler, blas_trans_type transt) const + float& rcon, + solve_singularity_handler sing_handler, blas_trans_type transt) const { MatrixType mattype (*this); return solve (mattype, b, info, rcon, sing_handler, transt); @@ -2538,7 +2538,7 @@ FloatComplexMatrix FloatComplexMatrix::lssolve (const FloatMatrix& b, octave_idx_type& info, - octave_idx_type& rank) const + octave_idx_type& rank) const { float rcon; return lssolve (FloatComplexMatrix (b), info, rank, rcon); @@ -2546,7 +2546,7 @@ FloatComplexMatrix FloatComplexMatrix::lssolve (const FloatMatrix& b, octave_idx_type& info, - octave_idx_type& rank, float& rcon) const + octave_idx_type& rank, float& rcon) const { return lssolve (FloatComplexMatrix (b), info, rank, rcon); } @@ -2570,7 +2570,7 @@ FloatComplexMatrix FloatComplexMatrix::lssolve (const FloatComplexMatrix& b, octave_idx_type& info, - octave_idx_type& rank) const + octave_idx_type& rank) const { float rcon; return lssolve (b, info, rank, rcon); @@ -2578,7 +2578,7 @@ FloatComplexMatrix FloatComplexMatrix::lssolve (const FloatComplexMatrix& b, octave_idx_type& info, - octave_idx_type& rank, float& rcon) const + octave_idx_type& rank, float& rcon) const { FloatComplexMatrix retval; @@ -2599,15 +2599,15 @@ rcon = -1.0; if (m != n) - { - retval = FloatComplexMatrix (maxmn, nrhs); - - for (octave_idx_type j = 0; j < nrhs; j++) - for (octave_idx_type i = 0; i < m; i++) - retval.elem (i, j) = b.elem (i, j); - } + { + retval = FloatComplexMatrix (maxmn, nrhs); + + for (octave_idx_type j = 0; j < nrhs; j++) + for (octave_idx_type i = 0; i < m; i++) + retval.elem (i, j) = b.elem (i, j); + } else - retval = b; + retval = b; FloatComplexMatrix atmp = *this; FloatComplex *tmp_data = atmp.fortran_vec (); @@ -2623,17 +2623,17 @@ octave_idx_type smlsiz; F77_FUNC (xilaenv, XILAENV) (9, F77_CONST_CHAR_ARG2 ("CGELSD", 6), - F77_CONST_CHAR_ARG2 (" ", 1), - 0, 0, 0, 0, smlsiz - F77_CHAR_ARG_LEN (6) - F77_CHAR_ARG_LEN (1)); + F77_CONST_CHAR_ARG2 (" ", 1), + 0, 0, 0, 0, smlsiz + F77_CHAR_ARG_LEN (6) + F77_CHAR_ARG_LEN (1)); octave_idx_type mnthr; F77_FUNC (xilaenv, XILAENV) (6, F77_CONST_CHAR_ARG2 ("CGELSD", 6), - F77_CONST_CHAR_ARG2 (" ", 1), - m, n, nrhs, -1, mnthr - F77_CHAR_ARG_LEN (6) - F77_CHAR_ARG_LEN (1)); + F77_CONST_CHAR_ARG2 (" ", 1), + m, n, nrhs, -1, mnthr + F77_CHAR_ARG_LEN (6) + F77_CHAR_ARG_LEN (1)); // We compute the size of rwork and iwork because ZGELSD in // older versions of LAPACK does not return them on a query @@ -2647,72 +2647,72 @@ #endif octave_idx_type nlvl = static_cast<octave_idx_type> (tmp) + 1; if (nlvl < 0) - nlvl = 0; + nlvl = 0; octave_idx_type lrwork = minmn*(10 + 2*smlsiz + 8*nlvl) - + 3*smlsiz*nrhs + (smlsiz+1)*(smlsiz+1); + + 3*smlsiz*nrhs + (smlsiz+1)*(smlsiz+1); if (lrwork < 1) - lrwork = 1; + lrwork = 1; Array<float> rwork (lrwork); float *prwork = rwork.fortran_vec (); octave_idx_type liwork = 3 * minmn * nlvl + 11 * minmn; if (liwork < 1) - liwork = 1; + liwork = 1; Array<octave_idx_type> iwork (liwork); octave_idx_type* piwork = iwork.fortran_vec (); F77_XFCN (cgelsd, CGELSD, (m, n, nrhs, tmp_data, m, pretval, maxmn, - ps, rcon, rank, work.fortran_vec (), - lwork, prwork, piwork, info)); + ps, rcon, rank, work.fortran_vec (), + lwork, prwork, piwork, info)); // The workspace query is broken in at least LAPACK 3.0.0 // through 3.1.1 when n >= mnthr. The obtuse formula below // should provide sufficient workspace for ZGELSD to operate // efficiently. if (n >= mnthr) - { - octave_idx_type addend = m; - - if (2*m-4 > addend) - addend = 2*m-4; - - if (nrhs > addend) - addend = nrhs; - - if (n-3*m > addend) - addend = n-3*m; - - const octave_idx_type lworkaround = 4*m + m*m + addend; - - if (std::real (work(0)) < lworkaround) - work(0) = lworkaround; - } + { + octave_idx_type addend = m; + + if (2*m-4 > addend) + addend = 2*m-4; + + if (nrhs > addend) + addend = nrhs; + + if (n-3*m > addend) + addend = n-3*m; + + const octave_idx_type lworkaround = 4*m + m*m + addend; + + if (std::real (work(0)) < lworkaround) + work(0) = lworkaround; + } else if (m >= n) - { - octave_idx_type lworkaround = 2*m + m*nrhs; - - if (std::real (work(0)) < lworkaround) - work(0) = lworkaround; - } + { + octave_idx_type lworkaround = 2*m + m*nrhs; + + if (std::real (work(0)) < lworkaround) + work(0) = lworkaround; + } lwork = static_cast<octave_idx_type> (std::real (work(0))); work.resize (lwork); F77_XFCN (cgelsd, CGELSD, (m, n, nrhs, tmp_data, m, pretval, - maxmn, ps, rcon, rank, - work.fortran_vec (), lwork, - prwork, piwork, info)); + maxmn, ps, rcon, 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, rcon); + (*current_liboctave_warning_handler) + ("zgelsd: rank deficient %dx%d matrix, rank = %d, tol = %e", + m, n, rank, rcon); if (s.elem (0) == 0.0) - rcon = 0.0; + rcon = 0.0; else - rcon = s.elem (minmn - 1) / s.elem (0); + rcon = s.elem (minmn - 1) / s.elem (0); retval.resize (n, nrhs); } @@ -2739,7 +2739,7 @@ FloatComplexColumnVector FloatComplexMatrix::lssolve (const FloatColumnVector& b, octave_idx_type& info, - octave_idx_type& rank) const + octave_idx_type& rank) const { float rcon; return lssolve (FloatComplexColumnVector (b), info, rank, rcon); @@ -2747,7 +2747,7 @@ FloatComplexColumnVector FloatComplexMatrix::lssolve (const FloatColumnVector& b, octave_idx_type& info, - octave_idx_type& rank, float& rcon) const + octave_idx_type& rank, float& rcon) const { return lssolve (FloatComplexColumnVector (b), info, rank, rcon); } @@ -2771,7 +2771,7 @@ FloatComplexColumnVector FloatComplexMatrix::lssolve (const FloatComplexColumnVector& b, octave_idx_type& info, - octave_idx_type& rank) const + octave_idx_type& rank) const { float rcon; return lssolve (b, info, rank, rcon); @@ -2780,7 +2780,7 @@ FloatComplexColumnVector FloatComplexMatrix::lssolve (const FloatComplexColumnVector& b, octave_idx_type& info, - octave_idx_type& rank, float& rcon) const + octave_idx_type& rank, float& rcon) const { FloatComplexColumnVector retval; @@ -2801,14 +2801,14 @@ rcon = -1.0; if (m != n) - { - retval = FloatComplexColumnVector (maxmn); - - for (octave_idx_type i = 0; i < m; i++) - retval.elem (i) = b.elem (i); - } + { + retval = FloatComplexColumnVector (maxmn); + + for (octave_idx_type i = 0; i < m; i++) + retval.elem (i) = b.elem (i); + } else - retval = b; + retval = b; FloatComplexMatrix atmp = *this; FloatComplex *tmp_data = atmp.fortran_vec (); @@ -2824,10 +2824,10 @@ octave_idx_type smlsiz; F77_FUNC (xilaenv, XILAENV) (9, F77_CONST_CHAR_ARG2 ("CGELSD", 6), - F77_CONST_CHAR_ARG2 (" ", 1), - 0, 0, 0, 0, smlsiz - F77_CHAR_ARG_LEN (6) - F77_CHAR_ARG_LEN (1)); + F77_CONST_CHAR_ARG2 (" ", 1), + 0, 0, 0, 0, smlsiz + F77_CHAR_ARG_LEN (6) + F77_CHAR_ARG_LEN (1)); // We compute the size of rwork and iwork because ZGELSD in // older versions of LAPACK does not return them on a query @@ -2841,24 +2841,24 @@ #endif octave_idx_type nlvl = static_cast<octave_idx_type> (tmp) + 1; if (nlvl < 0) - nlvl = 0; + nlvl = 0; octave_idx_type lrwork = minmn*(10 + 2*smlsiz + 8*nlvl) - + 3*smlsiz*nrhs + (smlsiz+1)*(smlsiz+1); + + 3*smlsiz*nrhs + (smlsiz+1)*(smlsiz+1); if (lrwork < 1) - lrwork = 1; + lrwork = 1; Array<float> rwork (lrwork); float *prwork = rwork.fortran_vec (); octave_idx_type liwork = 3 * minmn * nlvl + 11 * minmn; if (liwork < 1) - liwork = 1; + liwork = 1; Array<octave_idx_type> iwork (liwork); octave_idx_type* piwork = iwork.fortran_vec (); F77_XFCN (cgelsd, CGELSD, (m, n, nrhs, tmp_data, m, pretval, maxmn, - ps, rcon, rank, work.fortran_vec (), - lwork, prwork, piwork, info)); + ps, rcon, rank, work.fortran_vec (), + lwork, prwork, piwork, info)); lwork = static_cast<octave_idx_type> (std::real (work(0))); work.resize (lwork); @@ -2866,24 +2866,24 @@ iwork.resize (iwork(0)); F77_XFCN (cgelsd, CGELSD, (m, n, nrhs, tmp_data, m, pretval, - maxmn, ps, rcon, rank, - work.fortran_vec (), lwork, - prwork, piwork, info)); + maxmn, ps, rcon, rank, + work.fortran_vec (), lwork, + prwork, piwork, info)); if (rank < minmn) - { - if (rank < minmn) - (*current_liboctave_warning_handler) - ("zgelsd: rank deficient %dx%d matrix, rank = %d, tol = %e", - m, n, rank, rcon); - - if (s.elem (0) == 0.0) - rcon = 0.0; - else - rcon = 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, rcon); + + if (s.elem (0) == 0.0) + rcon = 0.0; + else + rcon = s.elem (minmn - 1) / s.elem (0); + + retval.resize (n, nrhs); + } } return retval; @@ -2920,11 +2920,11 @@ FloatComplex *c = retval.fortran_vec (); F77_XFCN (cgemm, CGEMM, (F77_CONST_CHAR_ARG2 ("N", 1), - F77_CONST_CHAR_ARG2 ("N", 1), - len, a_len, 1, 1.0, v.data (), len, - a.data (), 1, 0.0, c, len - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); + F77_CONST_CHAR_ARG2 ("N", 1), + len, a_len, 1, 1.0, v.data (), len, + a.data (), 1, 0.0, c, len + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); } return retval; @@ -3085,9 +3085,9 @@ for (octave_idx_type j = 0; j < nc; j++) for (octave_idx_type i = 0; i < nr; i++) { - FloatComplex val = elem (i, j); - if (xisnan (val)) - return true; + FloatComplex val = elem (i, j); + if (xisnan (val)) + return true; } return false; @@ -3102,9 +3102,9 @@ for (octave_idx_type j = 0; j < nc; j++) for (octave_idx_type i = 0; i < nr; i++) { - FloatComplex val = elem (i, j); - if (xisinf (val) || xisnan (val)) - return true; + FloatComplex val = elem (i, j); + if (xisinf (val) || xisnan (val)) + return true; } return false; @@ -3139,10 +3139,10 @@ min_val = r_val; if (i_val > max_val) - max_val = i_val; + max_val = i_val; if (i_val < max_val) - min_val = i_val; + min_val = i_val; } else return false; @@ -3150,25 +3150,25 @@ for (octave_idx_type j = 0; j < nc; j++) for (octave_idx_type i = 0; i < nr; i++) { - FloatComplex val = elem (i, j); - - float r_val = std::real (val); - float i_val = std::imag (val); - - if (r_val > max_val) - max_val = r_val; - - if (i_val > max_val) - max_val = i_val; - - if (r_val < min_val) - min_val = r_val; - - if (i_val < min_val) - min_val = i_val; - - if (D_NINT (r_val) != r_val || D_NINT (i_val) != i_val) - return false; + FloatComplex val = elem (i, j); + + float r_val = std::real (val); + float i_val = std::imag (val); + + if (r_val > max_val) + max_val = r_val; + + if (i_val > max_val) + max_val = i_val; + + if (r_val < min_val) + min_val = r_val; + + if (i_val < min_val) + min_val = i_val; + + if (D_NINT (r_val) != r_val || D_NINT (i_val) != i_val) + return false; } return true; @@ -3183,16 +3183,16 @@ for (octave_idx_type j = 0; j < nc; j++) for (octave_idx_type i = 0; i < nr; i++) { - FloatComplex val = elem (i, j); - - float r_val = std::real (val); - float i_val = std::imag (val); - - if ((! (xisnan (r_val) || xisinf (r_val)) - && fabs (r_val) > FLT_MAX) - || (! (xisnan (i_val) || xisinf (i_val)) - && fabs (i_val) > FLT_MAX)) - return true; + FloatComplex val = elem (i, j); + + float r_val = std::real (val); + float i_val = std::imag (val); + + if ((! (xisnan (r_val) || xisinf (r_val)) + && fabs (r_val) > FLT_MAX) + || (! (xisnan (i_val) || xisinf (i_val)) + && fabs (i_val) > FLT_MAX)) + return true; } return false; @@ -3265,13 +3265,13 @@ for (octave_idx_type j = 0; j < nc; j++) { if (std::imag (elem (i, j)) != 0.0) - { - retval = false; - break; - } + { + retval = false; + break; + } } - return retval; + return retval; } bool @@ -3284,13 +3284,13 @@ for (octave_idx_type i = 0; i < nr; i++) { if (std::imag (elem (i, j)) != 0.0) - { - retval = false; - break; - } + { + retval = false; + break; + } } - return retval; + return retval; } FloatComplexColumnVector @@ -3315,52 +3315,52 @@ for (octave_idx_type i = 0; i < nr; i++) { - bool real_only = row_is_real_only (i); - - octave_idx_type idx_j; - - FloatComplex tmp_min; - - float abs_min = octave_Float_NaN; - - for (idx_j = 0; idx_j < nc; idx_j++) - { - tmp_min = elem (i, idx_j); - - if (! xisnan (tmp_min)) - { - abs_min = real_only ? std::real (tmp_min) : std::abs (tmp_min); - break; - } - } - - for (octave_idx_type j = idx_j+1; j < nc; j++) - { - FloatComplex tmp = elem (i, j); - - if (xisnan (tmp)) - continue; - - float abs_tmp = real_only ? std::real (tmp) : std::abs (tmp); - - if (abs_tmp < abs_min) - { - idx_j = j; - tmp_min = tmp; - abs_min = abs_tmp; - } - } - - if (xisnan (tmp_min)) - { - result.elem (i) = FloatComplex_NaN_result; - idx_arg.elem (i) = 0; - } - else - { - result.elem (i) = tmp_min; - idx_arg.elem (i) = idx_j; - } + bool real_only = row_is_real_only (i); + + octave_idx_type idx_j; + + FloatComplex tmp_min; + + float abs_min = octave_Float_NaN; + + for (idx_j = 0; idx_j < nc; idx_j++) + { + tmp_min = elem (i, idx_j); + + if (! xisnan (tmp_min)) + { + abs_min = real_only ? std::real (tmp_min) : std::abs (tmp_min); + break; + } + } + + for (octave_idx_type j = idx_j+1; j < nc; j++) + { + FloatComplex tmp = elem (i, j); + + if (xisnan (tmp)) + continue; + + float abs_tmp = real_only ? std::real (tmp) : std::abs (tmp); + + if (abs_tmp < abs_min) + { + idx_j = j; + tmp_min = tmp; + abs_min = abs_tmp; + } + } + + if (xisnan (tmp_min)) + { + result.elem (i) = FloatComplex_NaN_result; + idx_arg.elem (i) = 0; + } + else + { + result.elem (i) = tmp_min; + idx_arg.elem (i) = idx_j; + } } } @@ -3389,52 +3389,52 @@ for (octave_idx_type i = 0; i < nr; i++) { - bool real_only = row_is_real_only (i); - - octave_idx_type idx_j; - - FloatComplex tmp_max; - - float abs_max = octave_Float_NaN; - - for (idx_j = 0; idx_j < nc; idx_j++) - { - tmp_max = elem (i, idx_j); - - if (! xisnan (tmp_max)) - { - abs_max = real_only ? std::real (tmp_max) : std::abs (tmp_max); - break; - } - } - - for (octave_idx_type j = idx_j+1; j < nc; j++) - { - FloatComplex tmp = elem (i, j); - - if (xisnan (tmp)) - continue; - - float abs_tmp = real_only ? std::real (tmp) : std::abs (tmp); - - if (abs_tmp > abs_max) - { - idx_j = j; - tmp_max = tmp; - abs_max = abs_tmp; - } - } - - if (xisnan (tmp_max)) - { - result.elem (i) = FloatComplex_NaN_result; - idx_arg.elem (i) = 0; - } - else - { - result.elem (i) = tmp_max; - idx_arg.elem (i) = idx_j; - } + bool real_only = row_is_real_only (i); + + octave_idx_type idx_j; + + FloatComplex tmp_max; + + float abs_max = octave_Float_NaN; + + for (idx_j = 0; idx_j < nc; idx_j++) + { + tmp_max = elem (i, idx_j); + + if (! xisnan (tmp_max)) + { + abs_max = real_only ? std::real (tmp_max) : std::abs (tmp_max); + break; + } + } + + for (octave_idx_type j = idx_j+1; j < nc; j++) + { + FloatComplex tmp = elem (i, j); + + if (xisnan (tmp)) + continue; + + float abs_tmp = real_only ? std::real (tmp) : std::abs (tmp); + + if (abs_tmp > abs_max) + { + idx_j = j; + tmp_max = tmp; + abs_max = abs_tmp; + } + } + + if (xisnan (tmp_max)) + { + result.elem (i) = FloatComplex_NaN_result; + idx_arg.elem (i) = 0; + } + else + { + result.elem (i) = tmp_max; + idx_arg.elem (i) = idx_j; + } } } @@ -3463,52 +3463,52 @@ for (octave_idx_type j = 0; j < nc; j++) { - bool real_only = column_is_real_only (j); - - octave_idx_type idx_i; - - FloatComplex tmp_min; - - float abs_min = octave_Float_NaN; - - for (idx_i = 0; idx_i < nr; idx_i++) - { - tmp_min = elem (idx_i, j); - - if (! xisnan (tmp_min)) - { - abs_min = real_only ? std::real (tmp_min) : std::abs (tmp_min); - break; - } - } - - for (octave_idx_type i = idx_i+1; i < nr; i++) - { - FloatComplex tmp = elem (i, j); - - if (xisnan (tmp)) - continue; - - float abs_tmp = real_only ? std::real (tmp) : std::abs (tmp); - - if (abs_tmp < abs_min) - { - idx_i = i; - tmp_min = tmp; - abs_min = abs_tmp; - } - } - - if (xisnan (tmp_min)) - { - result.elem (j) = FloatComplex_NaN_result; - idx_arg.elem (j) = 0; - } - else - { - result.elem (j) = tmp_min; - idx_arg.elem (j) = idx_i; - } + bool real_only = column_is_real_only (j); + + octave_idx_type idx_i; + + FloatComplex tmp_min; + + float abs_min = octave_Float_NaN; + + for (idx_i = 0; idx_i < nr; idx_i++) + { + tmp_min = elem (idx_i, j); + + if (! xisnan (tmp_min)) + { + abs_min = real_only ? std::real (tmp_min) : std::abs (tmp_min); + break; + } + } + + for (octave_idx_type i = idx_i+1; i < nr; i++) + { + FloatComplex tmp = elem (i, j); + + if (xisnan (tmp)) + continue; + + float abs_tmp = real_only ? std::real (tmp) : std::abs (tmp); + + if (abs_tmp < abs_min) + { + idx_i = i; + tmp_min = tmp; + abs_min = abs_tmp; + } + } + + if (xisnan (tmp_min)) + { + result.elem (j) = FloatComplex_NaN_result; + idx_arg.elem (j) = 0; + } + else + { + result.elem (j) = tmp_min; + idx_arg.elem (j) = idx_i; + } } } @@ -3537,52 +3537,52 @@ for (octave_idx_type j = 0; j < nc; j++) { - bool real_only = column_is_real_only (j); - - octave_idx_type idx_i; - - FloatComplex tmp_max; - - float abs_max = octave_Float_NaN; - - for (idx_i = 0; idx_i < nr; idx_i++) - { - tmp_max = elem (idx_i, j); - - if (! xisnan (tmp_max)) - { - abs_max = real_only ? std::real (tmp_max) : std::abs (tmp_max); - break; - } - } - - for (octave_idx_type i = idx_i+1; i < nr; i++) - { - FloatComplex tmp = elem (i, j); - - if (xisnan (tmp)) - continue; - - float abs_tmp = real_only ? std::real (tmp) : std::abs (tmp); - - if (abs_tmp > abs_max) - { - idx_i = i; - tmp_max = tmp; - abs_max = abs_tmp; - } - } - - if (xisnan (tmp_max)) - { - result.elem (j) = FloatComplex_NaN_result; - idx_arg.elem (j) = 0; - } - else - { - result.elem (j) = tmp_max; - idx_arg.elem (j) = idx_i; - } + bool real_only = column_is_real_only (j); + + octave_idx_type idx_i; + + FloatComplex tmp_max; + + float abs_max = octave_Float_NaN; + + for (idx_i = 0; idx_i < nr; idx_i++) + { + tmp_max = elem (idx_i, j); + + if (! xisnan (tmp_max)) + { + abs_max = real_only ? std::real (tmp_max) : std::abs (tmp_max); + break; + } + } + + for (octave_idx_type i = idx_i+1; i < nr; i++) + { + FloatComplex tmp = elem (i, j); + + if (xisnan (tmp)) + continue; + + float abs_tmp = real_only ? std::real (tmp) : std::abs (tmp); + + if (abs_tmp > abs_max) + { + idx_i = i; + tmp_max = tmp; + abs_max = abs_tmp; + } + } + + if (xisnan (tmp_max)) + { + result.elem (j) = FloatComplex_NaN_result; + idx_arg.elem (j) = 0; + } + else + { + result.elem (j) = tmp_max; + idx_arg.elem (j) = idx_i; + } } } @@ -3597,10 +3597,10 @@ for (octave_idx_type i = 0; i < a.rows (); i++) { for (octave_idx_type j = 0; j < a.cols (); j++) - { - os << " "; - octave_write_complex (os, a.elem (i, j)); - } + { + os << " "; + octave_write_complex (os, a.elem (i, j)); + } os << "\n"; } return os; @@ -3616,14 +3616,14 @@ { FloatComplex tmp; for (octave_idx_type i = 0; i < nr; i++) - for (octave_idx_type j = 0; j < nc; j++) - { - tmp = octave_read_value<FloatComplex> (is); - if (is) - a.elem (i, j) = tmp; - else - goto done; - } + for (octave_idx_type j = 0; j < nc; j++) + { + tmp = octave_read_value<FloatComplex> (is); + if (is) + a.elem (i, j) = tmp; + else + goto done; + } } done: @@ -3651,7 +3651,7 @@ FloatComplexMatrix Sylvester (const FloatComplexMatrix& a, const FloatComplexMatrix& b, - const FloatComplexMatrix& c) + const FloatComplexMatrix& c) { FloatComplexMatrix retval; @@ -3687,11 +3687,11 @@ FloatComplex *px = cx.fortran_vec (); F77_XFCN (ctrsyl, CTRSYL, (F77_CONST_CHAR_ARG2 ("N", 1), - F77_CONST_CHAR_ARG2 ("N", 1), - 1, a_nr, b_nr, pa, a_nr, pb, - b_nr, px, a_nr, scale, info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); + F77_CONST_CHAR_ARG2 ("N", 1), + 1, a_nr, b_nr, pa, a_nr, pb, + b_nr, px, a_nr, scale, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); // FIXME -- check info? @@ -3768,13 +3768,13 @@ else { if (a_nr == 0 || a_nc == 0 || b_nc == 0) - retval = FloatComplexMatrix (a_nr, b_nc, 0.0); + retval = FloatComplexMatrix (a_nr, b_nc, 0.0); else if (a.data () == b.data () && a_nr == b_nc && tra != trb) { - octave_idx_type lda = a.rows (); + octave_idx_type lda = a.rows (); retval = FloatComplexMatrix (a_nr, b_nc); - FloatComplex *c = retval.fortran_vec (); + FloatComplex *c = retval.fortran_vec (); const char *ctra = get_blas_trans_arg (tra, cja); if (cja || cjb) @@ -3805,15 +3805,15 @@ } else - { - octave_idx_type lda = a.rows (), tda = a.cols (); - octave_idx_type ldb = b.rows (), tdb = b.cols (); - - retval = FloatComplexMatrix (a_nr, b_nc); - FloatComplex *c = retval.fortran_vec (); - - if (b_nc == 1 && a_nr == 1) - { + { + octave_idx_type lda = a.rows (), tda = a.cols (); + octave_idx_type ldb = b.rows (), tdb = b.cols (); + + retval = FloatComplexMatrix (a_nr, b_nc); + FloatComplex *c = retval.fortran_vec (); + + if (b_nc == 1 && a_nr == 1) + { if (cja == cjb) { F77_FUNC (xcdotu, XCDOTU) (a_nc, a.data (), 1, b.data (), 1, *c); @@ -3840,18 +3840,18 @@ a.data (), 1, 0.0, c, 1 F77_CHAR_ARG_LEN (1))); } - else - { + else + { const char *ctra = get_blas_trans_arg (tra, cja); const char *ctrb = get_blas_trans_arg (trb, cjb); - F77_XFCN (cgemm, CGEMM, (F77_CONST_CHAR_ARG2 (ctra, 1), - F77_CONST_CHAR_ARG2 (ctrb, 1), - a_nr, b_nc, a_nc, 1.0, a.data (), - lda, b.data (), ldb, 0.0, c, a_nr - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); - } - } + F77_XFCN (cgemm, CGEMM, (F77_CONST_CHAR_ARG2 (ctra, 1), + F77_CONST_CHAR_ARG2 (ctrb, 1), + a_nr, b_nc, a_nc, 1.0, a.data (), + lda, b.data (), ldb, 0.0, c, a_nr + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + } + } } return retval; @@ -3883,8 +3883,8 @@ for (octave_idx_type j = 0; j < nc; j++) for (octave_idx_type i = 0; i < nr; i++) { - octave_quit (); - result (i, j) = xmin (c, m (i, j)); + octave_quit (); + result (i, j) = xmin (c, m (i, j)); } return result; @@ -3903,8 +3903,8 @@ for (octave_idx_type j = 0; j < nc; j++) for (octave_idx_type i = 0; i < nr; i++) { - octave_quit (); - result (i, j) = xmin (m (i, j), c); + octave_quit (); + result (i, j) = xmin (m (i, j), c); } return result; @@ -3919,7 +3919,7 @@ if (nr != b.rows () || nc != b.columns ()) { (*current_liboctave_error_handler) - ("two-arg min expecting args of same size"); + ("two-arg min expecting args of same size"); return FloatComplexMatrix (); } @@ -3931,28 +3931,28 @@ { int columns_are_real_only = 1; for (octave_idx_type i = 0; i < nr; i++) - { - octave_quit (); - if (std::imag (a (i, j)) != 0.0 || std::imag (b (i, j)) != 0.0) - { - columns_are_real_only = 0; - break; - } - } + { + octave_quit (); + if (std::imag (a (i, j)) != 0.0 || std::imag (b (i, j)) != 0.0) + { + columns_are_real_only = 0; + break; + } + } if (columns_are_real_only) - { - for (octave_idx_type i = 0; i < nr; i++) - result (i, j) = xmin (std::real (a (i, j)), std::real (b (i, j))); - } + { + for (octave_idx_type i = 0; i < nr; i++) + result (i, j) = xmin (std::real (a (i, j)), std::real (b (i, j))); + } else - { - for (octave_idx_type i = 0; i < nr; i++) - { - octave_quit (); - result (i, j) = xmin (a (i, j), b (i, j)); - } - } + { + for (octave_idx_type i = 0; i < nr; i++) + { + octave_quit (); + result (i, j) = xmin (a (i, j), b (i, j)); + } + } } return result; @@ -3971,8 +3971,8 @@ for (octave_idx_type j = 0; j < nc; j++) for (octave_idx_type i = 0; i < nr; i++) { - octave_quit (); - result (i, j) = xmax (c, m (i, j)); + octave_quit (); + result (i, j) = xmax (c, m (i, j)); } return result; @@ -3991,8 +3991,8 @@ for (octave_idx_type j = 0; j < nc; j++) for (octave_idx_type i = 0; i < nr; i++) { - octave_quit (); - result (i, j) = xmax (m (i, j), c); + octave_quit (); + result (i, j) = xmax (m (i, j), c); } return result; @@ -4007,7 +4007,7 @@ if (nr != b.rows () || nc != b.columns ()) { (*current_liboctave_error_handler) - ("two-arg max expecting args of same size"); + ("two-arg max expecting args of same size"); return FloatComplexMatrix (); } @@ -4019,31 +4019,31 @@ { int columns_are_real_only = 1; for (octave_idx_type i = 0; i < nr; i++) - { - octave_quit (); - if (std::imag (a (i, j)) != 0.0 || std::imag (b (i, j)) != 0.0) - { - columns_are_real_only = 0; - break; - } - } + { + octave_quit (); + if (std::imag (a (i, j)) != 0.0 || std::imag (b (i, j)) != 0.0) + { + columns_are_real_only = 0; + break; + } + } if (columns_are_real_only) - { - for (octave_idx_type i = 0; i < nr; i++) - { - octave_quit (); - result (i, j) = xmax (std::real (a (i, j)), std::real (b (i, j))); - } - } + { + for (octave_idx_type i = 0; i < nr; i++) + { + octave_quit (); + result (i, j) = xmax (std::real (a (i, j)), std::real (b (i, j))); + } + } else - { - for (octave_idx_type i = 0; i < nr; i++) - { - octave_quit (); - result (i, j) = xmax (a (i, j), b (i, j)); - } - } + { + for (octave_idx_type i = 0; i < nr; i++) + { + octave_quit (); + result (i, j) = xmax (a (i, j), b (i, j)); + } + } } return result;