Mercurial > octave-nkf
diff liboctave/fMatrix.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/fMatrix.cc Thu Feb 11 12:16:43 2010 -0500 +++ b/liboctave/fMatrix.cc Thu Feb 11 12:23:32 2010 -0500 @@ -63,151 +63,151 @@ { 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 (sgebal, SGEBAL) (F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, float*, const octave_idx_type&, octave_idx_type&, - octave_idx_type&, float*, octave_idx_type& - F77_CHAR_ARG_LEN_DECL); + const octave_idx_type&, float*, 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 (sgemm, SGEMM) (F77_CONST_CHAR_ARG_DECL, - F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, - const float&, const float*, const octave_idx_type&, - const float*, const octave_idx_type&, const float&, - float*, 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 float&, const float*, const octave_idx_type&, + const float*, const octave_idx_type&, const float&, + float*, const octave_idx_type& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); F77_RET_T F77_FUNC (sgemv, SGEMV) (F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, const octave_idx_type&, const float&, - const float*, const octave_idx_type&, const float*, - const octave_idx_type&, const float&, float*, - const octave_idx_type& - F77_CHAR_ARG_LEN_DECL); + const octave_idx_type&, const octave_idx_type&, const float&, + const float*, const octave_idx_type&, const float*, + const octave_idx_type&, const float&, float*, + const octave_idx_type& + F77_CHAR_ARG_LEN_DECL); F77_RET_T F77_FUNC (xsdot, XSDOT) (const octave_idx_type&, const float*, const octave_idx_type&, - const float*, const octave_idx_type&, float&); + const float*, const octave_idx_type&, float&); F77_RET_T F77_FUNC (ssyrk, SSYRK) (F77_CONST_CHAR_ARG_DECL, - F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, const octave_idx_type&, - const float&, const float*, const octave_idx_type&, - const float&, float*, 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 float*, const octave_idx_type&, + const float&, float*, const octave_idx_type& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); F77_RET_T F77_FUNC (sgetrf, SGETRF) (const octave_idx_type&, const octave_idx_type&, float*, const octave_idx_type&, - octave_idx_type*, octave_idx_type&); + octave_idx_type*, octave_idx_type&); F77_RET_T F77_FUNC (sgetrs, SGETRS) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, const octave_idx_type&, - const float*, const octave_idx_type&, - const octave_idx_type*, float*, const octave_idx_type&, octave_idx_type& - F77_CHAR_ARG_LEN_DECL); + const float*, const octave_idx_type&, + const octave_idx_type*, float*, const octave_idx_type&, octave_idx_type& + F77_CHAR_ARG_LEN_DECL); F77_RET_T F77_FUNC (sgetri, SGETRI) (const octave_idx_type&, float*, const octave_idx_type&, const octave_idx_type*, - float*, const octave_idx_type&, octave_idx_type&); + float*, const octave_idx_type&, octave_idx_type&); F77_RET_T F77_FUNC (sgecon, SGECON) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, float*, - const octave_idx_type&, const float&, float&, - float*, octave_idx_type*, octave_idx_type& - F77_CHAR_ARG_LEN_DECL); + const octave_idx_type&, const float&, float&, + float*, octave_idx_type*, octave_idx_type& + F77_CHAR_ARG_LEN_DECL); F77_RET_T F77_FUNC (sgelsy, SGELSY) (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*, float&, octave_idx_type&, - float*, const octave_idx_type&, octave_idx_type&); + float*, const octave_idx_type&, float*, + const octave_idx_type&, octave_idx_type*, float&, octave_idx_type&, + float*, const octave_idx_type&, octave_idx_type&); F77_RET_T F77_FUNC (sgelsd, SGELSD) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, - float*, const octave_idx_type&, float*, - const octave_idx_type&, float*, float&, octave_idx_type&, - float*, const octave_idx_type&, octave_idx_type*, - octave_idx_type&); + float*, const octave_idx_type&, float*, + const octave_idx_type&, float*, float&, octave_idx_type&, + float*, const octave_idx_type&, octave_idx_type*, + octave_idx_type&); F77_RET_T F77_FUNC (spotrf, SPOTRF) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, - float *, const octave_idx_type&, - octave_idx_type& F77_CHAR_ARG_LEN_DECL); + float *, const octave_idx_type&, + octave_idx_type& F77_CHAR_ARG_LEN_DECL); F77_RET_T F77_FUNC (spocon, SPOCON) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, - float*, const octave_idx_type&, const float&, - float&, float*, octave_idx_type*, - octave_idx_type& F77_CHAR_ARG_LEN_DECL); + float*, const octave_idx_type&, const float&, + float&, float*, octave_idx_type*, + octave_idx_type& F77_CHAR_ARG_LEN_DECL); F77_RET_T F77_FUNC (spotrs, SPOTRS) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, - const octave_idx_type&, const float*, - const octave_idx_type&, float*, - const octave_idx_type&, octave_idx_type& - F77_CHAR_ARG_LEN_DECL); + const octave_idx_type&, const float*, + const octave_idx_type&, float*, + const octave_idx_type&, octave_idx_type& + F77_CHAR_ARG_LEN_DECL); F77_RET_T F77_FUNC (strtri, STRTRI) (F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, const float*, - const octave_idx_type&, octave_idx_type& - F77_CHAR_ARG_LEN_DECL - F77_CHAR_ARG_LEN_DECL); + const octave_idx_type&, const float*, + const octave_idx_type&, octave_idx_type& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); F77_RET_T F77_FUNC (strcon, STRCON) (F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, - F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, - const float*, const octave_idx_type&, float&, - float*, 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 float*, const octave_idx_type&, float&, + float*, 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 (strtrs, STRTRS) (F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL, - F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, - const octave_idx_type&, const float*, - const octave_idx_type&, float*, - 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 float*, + const octave_idx_type&, float*, + 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 (slartg, SLARTG) (const float&, const float&, float&, - float&, float&); + float&, float&); F77_RET_T F77_FUNC (strsyl, STRSYL) (F77_CONST_CHAR_ARG_DECL, - F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, - const float*, const octave_idx_type&, const float*, - const octave_idx_type&, const float*, 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 float*, const octave_idx_type&, const float*, + const octave_idx_type&, const float*, const octave_idx_type&, + float&, octave_idx_type& + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); F77_RET_T F77_FUNC (xslange, XSLANGE) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, - const octave_idx_type&, const float*, - const octave_idx_type&, float*, float& - F77_CHAR_ARG_LEN_DECL); + const octave_idx_type&, const float*, + const octave_idx_type&, float*, float& + F77_CHAR_ARG_LEN_DECL); } // Matrix class. @@ -282,9 +282,9 @@ if (is_square () && rows () > 0) { for (octave_idx_type i = 0; i < rows (); i++) - for (octave_idx_type j = i+1; j < cols (); j++) - if (elem (i, j) != elem (j, i)) - return false; + for (octave_idx_type j = i+1; j < cols (); j++) + if (elem (i, j) != elem (j, i)) + return false; return true; } @@ -315,7 +315,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; @@ -337,7 +337,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; @@ -364,7 +364,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; @@ -381,8 +381,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; @@ -409,8 +409,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; @@ -496,7 +496,7 @@ if (nc != a.cols ()) { (*current_liboctave_error_handler) - ("column dimension mismatch for stack"); + ("column dimension mismatch for stack"); return FloatMatrix (); } @@ -515,7 +515,7 @@ if (nc != a.length ()) { (*current_liboctave_error_handler) - ("column dimension mismatch for stack"); + ("column dimension mismatch for stack"); return FloatMatrix (); } @@ -534,7 +534,7 @@ if (nc != 1) { (*current_liboctave_error_handler) - ("column dimension mismatch for stack"); + ("column dimension mismatch for stack"); return FloatMatrix (); } @@ -553,7 +553,7 @@ if (nc != a.cols ()) { (*current_liboctave_error_handler) - ("column dimension mismatch for stack"); + ("column dimension mismatch for stack"); return FloatMatrix (); } @@ -641,7 +641,7 @@ FloatMatrix FloatMatrix::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); @@ -664,7 +664,7 @@ FloatMatrix FloatMatrix::tinverse (MatrixType &mattype, octave_idx_type& info, float& rcon, - int force, int calc_cond) const + int force, int calc_cond) const { FloatMatrix retval; @@ -682,38 +682,38 @@ float *tmp_data = retval.fortran_vec (); F77_XFCN (strtri, STRTRI, (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 dtrcon_info = 0; - char job = '1'; - - OCTAVE_LOCAL_BUFFER (float, work, 3 * nr); - OCTAVE_LOCAL_BUFFER (octave_idx_type, iwork, nr); - - F77_XFCN (strcon, STRCON, (F77_CONST_CHAR_ARG2 (&job, 1), - F77_CONST_CHAR_ARG2 (&uplo, 1), - F77_CONST_CHAR_ARG2 (&udiag, 1), - nr, tmp_data, nr, rcon, - 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; - } + { + octave_idx_type dtrcon_info = 0; + char job = '1'; + + OCTAVE_LOCAL_BUFFER (float, work, 3 * nr); + OCTAVE_LOCAL_BUFFER (octave_idx_type, iwork, nr); + + F77_XFCN (strcon, STRCON, (F77_CONST_CHAR_ARG2 (&job, 1), + F77_CONST_CHAR_ARG2 (&uplo, 1), + F77_CONST_CHAR_ARG2 (&udiag, 1), + nr, tmp_data, nr, rcon, + 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; + } if (info == -1 && ! force) - retval = *this; // Restore matrix contents. + retval = *this; // Restore matrix contents. } return retval; @@ -722,7 +722,7 @@ FloatMatrix FloatMatrix::finverse (MatrixType &mattype, octave_idx_type& info, float& rcon, - int force, int calc_cond) const + int force, int calc_cond) const { FloatMatrix retval; @@ -744,7 +744,7 @@ // Query the optimum work array size. F77_XFCN (sgetri, SGETRI, (nc, tmp_data, nr, pipvt, - z.fortran_vec (), lwork, info)); + z.fortran_vec (), lwork, info)); lwork = static_cast<octave_idx_type> (z(0)); lwork = (lwork < 2 *nc ? 2*nc : lwork); @@ -756,46 +756,46 @@ // Calculate the norm of the matrix, for later use. float anorm = 0; 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 (sgetrf, SGETRF, (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) - { - 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 (sgecon, SGECON, (F77_CONST_CHAR_ARG2 (&job, 1), - nc, tmp_data, nr, anorm, - rcon, pz, piz, dgecon_info - F77_CHAR_ARG_LEN (1))); - - if (dgecon_info != 0) - info = -1; - } + { + 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 (sgecon, SGECON, (F77_CONST_CHAR_ARG2 (&job, 1), + nc, tmp_data, nr, anorm, + rcon, pz, piz, dgecon_info + F77_CHAR_ARG_LEN (1))); + + if (dgecon_info != 0) + info = -1; + } if (info == -1 && ! force) - retval = *this; // Restore matrix contents. + retval = *this; // Restore matrix contents. else - { - octave_idx_type dgetri_info = 0; - - F77_XFCN (sgetri, SGETRI, (nc, tmp_data, nr, pipvt, - pz, lwork, dgetri_info)); - - if (dgetri_info != 0) - info = -1; - } + { + octave_idx_type dgetri_info = 0; + + F77_XFCN (sgetri, SGETRI, (nc, tmp_data, nr, pipvt, + pz, lwork, dgetri_info)); + + if (dgetri_info != 0) + info = -1; + } if (info != 0) - mattype.mark_as_rectangular(); + mattype.mark_as_rectangular(); } return retval; @@ -803,7 +803,7 @@ FloatMatrix FloatMatrix::inverse (MatrixType &mattype, octave_idx_type& info, float& rcon, - int force, int calc_cond) const + int force, int calc_cond) const { int typ = mattype.type (false); FloatMatrix ret; @@ -816,25 +816,25 @@ else { if (mattype.is_hermitian ()) - { - FloatCHOL 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 (); - } + { + FloatCHOL 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 = FloatMatrix (rows (), columns (), octave_Float_Inf); + ret = FloatMatrix (rows (), columns (), octave_Float_Inf); } return ret; @@ -858,9 +858,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) @@ -1122,12 +1122,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; @@ -1191,12 +1191,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; @@ -1383,143 +1383,143 @@ 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 float *tmp_data = fortran_vec (); - octave_idx_type info = 0; - char norm = '1'; - char uplo = 'U'; - char dia = 'N'; - - Array<float> z (3 * nc); - float *pz = z.fortran_vec (); - Array<octave_idx_type> iz (nc); - octave_idx_type *piz = iz.fortran_vec (); - - F77_XFCN (strcon, STRCON, (F77_CONST_CHAR_ARG2 (&norm, 1), - F77_CONST_CHAR_ARG2 (&uplo, 1), - F77_CONST_CHAR_ARG2 (&dia, 1), - nr, tmp_data, nr, rcon, - pz, piz, info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); - - if (info != 0) - rcon = 0.0; - } + { + const float *tmp_data = fortran_vec (); + octave_idx_type info = 0; + char norm = '1'; + char uplo = 'U'; + char dia = 'N'; + + Array<float> z (3 * nc); + float *pz = z.fortran_vec (); + Array<octave_idx_type> iz (nc); + octave_idx_type *piz = iz.fortran_vec (); + + F77_XFCN (strcon, STRCON, (F77_CONST_CHAR_ARG2 (&norm, 1), + F77_CONST_CHAR_ARG2 (&uplo, 1), + F77_CONST_CHAR_ARG2 (&dia, 1), + nr, tmp_data, nr, rcon, + pz, piz, 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_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 float *tmp_data = fortran_vec (); - octave_idx_type info = 0; - char norm = '1'; - char uplo = 'L'; - char dia = 'N'; - - Array<float> z (3 * nc); - float *pz = z.fortran_vec (); - Array<octave_idx_type> iz (nc); - octave_idx_type *piz = iz.fortran_vec (); - - F77_XFCN (strcon, STRCON, (F77_CONST_CHAR_ARG2 (&norm, 1), - F77_CONST_CHAR_ARG2 (&uplo, 1), - F77_CONST_CHAR_ARG2 (&dia, 1), - nr, tmp_data, nr, rcon, - pz, piz, info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); - - if (info != 0) - rcon = 0.0; - } + { + const float *tmp_data = fortran_vec (); + octave_idx_type info = 0; + char norm = '1'; + char uplo = 'L'; + char dia = 'N'; + + Array<float> z (3 * nc); + float *pz = z.fortran_vec (); + Array<octave_idx_type> iz (nc); + octave_idx_type *piz = iz.fortran_vec (); + + F77_XFCN (strcon, STRCON, (F77_CONST_CHAR_ARG2 (&norm, 1), + F77_CONST_CHAR_ARG2 (&uplo, 1), + F77_CONST_CHAR_ARG2 (&dia, 1), + nr, tmp_data, nr, rcon, + pz, piz, 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; - FloatMatrix atmp = *this; - float *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 (spotrf, SPOTRF, (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<float> z (3 * nc); - float *pz = z.fortran_vec (); - Array<octave_idx_type> iz (nc); - octave_idx_type *piz = iz.fortran_vec (); - - F77_XFCN (spocon, SPOCON, (F77_CONST_CHAR_ARG2 (&job, 1), - nr, tmp_data, nr, anorm, - rcon, pz, piz, 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<float> z (4 * nc); - float *pz = z.fortran_vec (); - Array<octave_idx_type> iz (nc); - octave_idx_type *piz = iz.fortran_vec (); - - F77_XFCN (sgetrf, SGETRF, (nr, nr, tmp_data, nr, pipvt, info)); - - if (info != 0) - { - rcon = 0.0; - mattype.mark_as_rectangular (); - } - else - { - char job = '1'; - F77_XFCN (sgecon, SGECON, (F77_CONST_CHAR_ARG2 (&job, 1), - nc, tmp_data, nr, anorm, - rcon, pz, piz, info - F77_CHAR_ARG_LEN (1))); - - if (info != 0) - rcon = 0.0; - } - } - } + { + float anorm = -1.0; + FloatMatrix atmp = *this; + float *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 (spotrf, SPOTRF, (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<float> z (3 * nc); + float *pz = z.fortran_vec (); + Array<octave_idx_type> iz (nc); + octave_idx_type *piz = iz.fortran_vec (); + + F77_XFCN (spocon, SPOCON, (F77_CONST_CHAR_ARG2 (&job, 1), + nr, tmp_data, nr, anorm, + rcon, pz, piz, 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<float> z (4 * nc); + float *pz = z.fortran_vec (); + Array<octave_idx_type> iz (nc); + octave_idx_type *piz = iz.fortran_vec (); + + F77_XFCN (sgetrf, SGETRF, (nr, nr, tmp_data, nr, pipvt, info)); + + if (info != 0) + { + rcon = 0.0; + mattype.mark_as_rectangular (); + } + else + { + char job = '1'; + F77_XFCN (sgecon, SGECON, (F77_CONST_CHAR_ARG2 (&job, 1), + nc, tmp_data, nr, anorm, + rcon, pz, piz, info + F77_CHAR_ARG_LEN (1))); + + if (info != 0) + rcon = 0.0; + } + } + } else - rcon = 0.0; + rcon = 0.0; } return rcon; @@ -1527,8 +1527,8 @@ FloatMatrix FloatMatrix::utsolve (MatrixType &mattype, const FloatMatrix& b, octave_idx_type& info, - float& rcon, solve_singularity_handler sing_handler, - bool calc_cond, blas_trans_type transt) const + float& rcon, solve_singularity_handler sing_handler, + bool calc_cond, blas_trans_type transt) const { FloatMatrix retval; @@ -1545,81 +1545,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 float *tmp_data = fortran_vec (); - - if (calc_cond) - { - char norm = '1'; - char uplo = 'U'; - char dia = 'N'; - - Array<float> z (3 * nc); - float *pz = z.fortran_vec (); - Array<octave_idx_type> iz (nc); - octave_idx_type *piz = iz.fortran_vec (); - - F77_XFCN (strcon, STRCON, (F77_CONST_CHAR_ARG2 (&norm, 1), - F77_CONST_CHAR_ARG2 (&uplo, 1), - F77_CONST_CHAR_ARG2 (&dia, 1), - nr, tmp_data, nr, rcon, - pz, piz, 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; - float *result = retval.fortran_vec (); - - char uplo = 'U'; - char trans = get_blas_char (transt); - char dia = 'N'; - - F77_XFCN (strtrs, STRTRS, (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 float *tmp_data = fortran_vec (); + + if (calc_cond) + { + char norm = '1'; + char uplo = 'U'; + char dia = 'N'; + + Array<float> z (3 * nc); + float *pz = z.fortran_vec (); + Array<octave_idx_type> iz (nc); + octave_idx_type *piz = iz.fortran_vec (); + + F77_XFCN (strcon, STRCON, (F77_CONST_CHAR_ARG2 (&norm, 1), + F77_CONST_CHAR_ARG2 (&uplo, 1), + F77_CONST_CHAR_ARG2 (&dia, 1), + nr, tmp_data, nr, rcon, + pz, piz, 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; + float *result = retval.fortran_vec (); + + char uplo = 'U'; + char trans = get_blas_char (transt); + char dia = 'N'; + + F77_XFCN (strtrs, STRTRS, (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; @@ -1627,8 +1627,8 @@ FloatMatrix FloatMatrix::ltsolve (MatrixType &mattype, const FloatMatrix& b, octave_idx_type& info, - float& rcon, solve_singularity_handler sing_handler, - bool calc_cond, blas_trans_type transt) const + float& rcon, solve_singularity_handler sing_handler, + bool calc_cond, blas_trans_type transt) const { FloatMatrix retval; @@ -1645,81 +1645,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 float *tmp_data = fortran_vec (); - - if (calc_cond) - { - char norm = '1'; - char uplo = 'L'; - char dia = 'N'; - - Array<float> z (3 * nc); - float *pz = z.fortran_vec (); - Array<octave_idx_type> iz (nc); - octave_idx_type *piz = iz.fortran_vec (); - - F77_XFCN (strcon, STRCON, (F77_CONST_CHAR_ARG2 (&norm, 1), - F77_CONST_CHAR_ARG2 (&uplo, 1), - F77_CONST_CHAR_ARG2 (&dia, 1), - nr, tmp_data, nr, rcon, - pz, piz, 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; - float *result = retval.fortran_vec (); - - char uplo = 'L'; - char trans = get_blas_char (transt); - char dia = 'N'; - - F77_XFCN (strtrs, STRTRS, (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 float *tmp_data = fortran_vec (); + + if (calc_cond) + { + char norm = '1'; + char uplo = 'L'; + char dia = 'N'; + + Array<float> z (3 * nc); + float *pz = z.fortran_vec (); + Array<octave_idx_type> iz (nc); + octave_idx_type *piz = iz.fortran_vec (); + + F77_XFCN (strcon, STRCON, (F77_CONST_CHAR_ARG2 (&norm, 1), + F77_CONST_CHAR_ARG2 (&uplo, 1), + F77_CONST_CHAR_ARG2 (&dia, 1), + nr, tmp_data, nr, rcon, + pz, piz, 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; + float *result = retval.fortran_vec (); + + char uplo = 'L'; + char trans = get_blas_char (transt); + char dia = 'N'; + + F77_XFCN (strtrs, STRTRS, (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; @@ -1727,8 +1727,8 @@ FloatMatrix FloatMatrix::fsolve (MatrixType &mattype, const FloatMatrix& b, octave_idx_type& info, - float& rcon, solve_singularity_handler sing_handler, - bool calc_cond) const + float& rcon, solve_singularity_handler sing_handler, + bool calc_cond) const { FloatMatrix retval; @@ -1748,160 +1748,160 @@ float anorm = -1.; if (typ == MatrixType::Hermitian) - { - info = 0; - char job = 'L'; - FloatMatrix atmp = *this; - float *tmp_data = atmp.fortran_vec (); - anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max(); - - F77_XFCN (spotrf, SPOTRF, (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<float> z (3 * nc); - float *pz = z.fortran_vec (); - Array<octave_idx_type> iz (nc); - octave_idx_type *piz = iz.fortran_vec (); - - F77_XFCN (spocon, SPOCON, (F77_CONST_CHAR_ARG2 (&job, 1), - nr, tmp_data, nr, anorm, - rcon, pz, piz, 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; - float *result = retval.fortran_vec (); - - octave_idx_type b_nc = b.cols (); - - F77_XFCN (spotrs, SPOTRS, (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'; + FloatMatrix atmp = *this; + float *tmp_data = atmp.fortran_vec (); + anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max(); + + F77_XFCN (spotrf, SPOTRF, (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<float> z (3 * nc); + float *pz = z.fortran_vec (); + Array<octave_idx_type> iz (nc); + octave_idx_type *piz = iz.fortran_vec (); + + F77_XFCN (spocon, SPOCON, (F77_CONST_CHAR_ARG2 (&job, 1), + nr, tmp_data, nr, anorm, + rcon, pz, piz, 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; + float *result = retval.fortran_vec (); + + octave_idx_type b_nc = b.cols (); + + F77_XFCN (spotrs, SPOTRS, (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 (); - - FloatMatrix atmp = *this; - float *tmp_data = atmp.fortran_vec (); - if(anorm < 0.) - anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max(); - - Array<float> z (4 * nc); - float *pz = z.fortran_vec (); - Array<octave_idx_type> iz (nc); - octave_idx_type *piz = iz.fortran_vec (); - - F77_XFCN (sgetrf, SGETRF, (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 (sgecon, SGECON, (F77_CONST_CHAR_ARG2 (&job, 1), - nc, tmp_data, nr, anorm, - rcon, pz, piz, 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; - float *result = retval.fortran_vec (); - - octave_idx_type b_nc = b.cols (); - - char job = 'N'; - F77_XFCN (sgetrs, SGETRS, (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 (); + + FloatMatrix atmp = *this; + float *tmp_data = atmp.fortran_vec (); + if(anorm < 0.) + anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max(); + + Array<float> z (4 * nc); + float *pz = z.fortran_vec (); + Array<octave_idx_type> iz (nc); + octave_idx_type *piz = iz.fortran_vec (); + + F77_XFCN (sgetrf, SGETRF, (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 (sgecon, SGECON, (F77_CONST_CHAR_ARG2 (&job, 1), + nc, tmp_data, nr, anorm, + rcon, pz, piz, 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; + float *result = retval.fortran_vec (); + + octave_idx_type b_nc = b.cols (); + + char job = 'N'; + F77_XFCN (sgetrs, SGETRS, (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) - (*current_liboctave_error_handler) ("incorrect matrix type"); + (*current_liboctave_error_handler) ("incorrect matrix type"); } return retval; @@ -1924,15 +1924,15 @@ FloatMatrix FloatMatrix::solve (MatrixType &typ, const FloatMatrix& b, octave_idx_type& info, - float& rcon) const + float& rcon) const { return solve (typ, b, info, rcon, 0); } FloatMatrix FloatMatrix::solve (MatrixType &mattype, 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 { FloatMatrix retval; int typ = mattype.type (); @@ -1983,7 +1983,7 @@ FloatComplexMatrix FloatMatrix::solve (MatrixType &typ, const FloatComplexMatrix& b, octave_idx_type& info, - float& rcon) const + float& rcon) const { return solve (typ, b, info, rcon, 0); } @@ -2017,8 +2017,8 @@ FloatComplexMatrix FloatMatrix::solve (MatrixType &typ, const FloatComplexMatrix& 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 { FloatMatrix tmp = stack_complex_matrix (b); tmp = solve (typ, tmp, info, rcon, sing_handler, singular_fallback, transt); @@ -2034,7 +2034,7 @@ FloatColumnVector FloatMatrix::solve (MatrixType &typ, const FloatColumnVector& b, - octave_idx_type& info) const + octave_idx_type& info) const { float rcon; return solve (typ, b, info, rcon); @@ -2042,14 +2042,14 @@ FloatColumnVector FloatMatrix::solve (MatrixType &typ, const FloatColumnVector& b, octave_idx_type& info, - float& rcon) const + float& rcon) const { return solve (typ, b, info, rcon, 0); } FloatColumnVector FloatMatrix::solve (MatrixType &typ, 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 { FloatMatrix tmp (b); return solve (typ, tmp, info, rcon, sing_handler, transt).column(static_cast<octave_idx_type> (0)); @@ -2064,7 +2064,7 @@ FloatComplexColumnVector FloatMatrix::solve (MatrixType &typ, const FloatComplexColumnVector& b, - octave_idx_type& info) const + octave_idx_type& info) const { FloatComplexMatrix tmp (*this); return tmp.solve (typ, b, info); @@ -2072,7 +2072,7 @@ FloatComplexColumnVector FloatMatrix::solve (MatrixType &typ, const FloatComplexColumnVector& b, - octave_idx_type& info, float& rcon) const + octave_idx_type& info, float& rcon) const { FloatComplexMatrix tmp (*this); return tmp.solve (typ, b, info, rcon); @@ -2080,8 +2080,8 @@ FloatComplexColumnVector FloatMatrix::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 (*this); return tmp.solve(typ, b, info, rcon, sing_handler, transt); @@ -2110,7 +2110,7 @@ FloatMatrix FloatMatrix::solve (const FloatMatrix& 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, true, transt); @@ -2139,7 +2139,7 @@ FloatComplexMatrix FloatMatrix::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 { FloatComplexMatrix tmp (*this); return tmp.solve (b, info, rcon, sing_handler, transt); @@ -2167,7 +2167,7 @@ FloatColumnVector FloatMatrix::solve (const FloatColumnVector& 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, transt); @@ -2196,7 +2196,7 @@ FloatComplexColumnVector FloatMatrix::solve (const FloatComplexColumnVector& 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 (*this); return tmp.solve (b, info, rcon, sing_handler, transt); @@ -2221,7 +2221,7 @@ FloatMatrix FloatMatrix::lssolve (const FloatMatrix& b, octave_idx_type& info, - octave_idx_type& rank) const + octave_idx_type& rank) const { float rcon; return lssolve (b, info, rank, rcon); @@ -2229,7 +2229,7 @@ FloatMatrix FloatMatrix::lssolve (const FloatMatrix& b, octave_idx_type& info, - octave_idx_type& rank, float &rcon) const + octave_idx_type& rank, float &rcon) const { FloatMatrix retval; @@ -2249,15 +2249,15 @@ octave_idx_type maxmn = m > n ? m : n; rcon = -1.0; if (m != n) - { - retval = FloatMatrix (maxmn, nrhs, 0.0); - - 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 = FloatMatrix (maxmn, nrhs, 0.0); + + 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; FloatMatrix atmp = *this; float *tmp_data = atmp.fortran_vec (); @@ -2273,17 +2273,17 @@ octave_idx_type smlsiz; F77_FUNC (xilaenv, XILAENV) (9, F77_CONST_CHAR_ARG2 ("SGELSD", 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 ("SGELSD", 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 iwork because DGELSD in older versions // of LAPACK does not return it on a query call. @@ -2296,70 +2296,70 @@ #endif octave_idx_type nlvl = static_cast<octave_idx_type> (tmp) + 1; if (nlvl < 0) - nlvl = 0; + nlvl = 0; 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 (sgelsd, SGELSD, (m, n, nrhs, tmp_data, m, pretval, maxmn, - ps, rcon, rank, work.fortran_vec (), - lwork, piwork, info)); + ps, rcon, rank, work.fortran_vec (), + lwork, 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 DGELSD to operate // efficiently. if (n >= mnthr) - { - const octave_idx_type wlalsd - = 9*m + 2*m*smlsiz + 8*m*nlvl + m*nrhs + (smlsiz+1)*(smlsiz+1); - - 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; - - if (wlalsd > addend) - addend = wlalsd; - - const octave_idx_type lworkaround = 4*m + m*m + addend; - - if (work(0) < lworkaround) - work(0) = lworkaround; - } + { + const octave_idx_type wlalsd + = 9*m + 2*m*smlsiz + 8*m*nlvl + m*nrhs + (smlsiz+1)*(smlsiz+1); + + 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; + + if (wlalsd > addend) + addend = wlalsd; + + const octave_idx_type lworkaround = 4*m + m*m + addend; + + if (work(0) < lworkaround) + work(0) = lworkaround; + } else if (m >= n) - { - octave_idx_type lworkaround - = 12*n + 2*n*smlsiz + 8*n*nlvl + n*nrhs + (smlsiz+1)*(smlsiz+1); - - if (work(0) < lworkaround) - work(0) = lworkaround; - } + { + octave_idx_type lworkaround + = 12*n + 2*n*smlsiz + 8*n*nlvl + n*nrhs + (smlsiz+1)*(smlsiz+1); + + if (work(0) < lworkaround) + work(0) = lworkaround; + } lwork = static_cast<octave_idx_type> (work(0)); work.resize (lwork); F77_XFCN (sgelsd, SGELSD, (m, n, nrhs, tmp_data, m, pretval, - maxmn, ps, rcon, rank, - work.fortran_vec (), lwork, - piwork, info)); + maxmn, ps, rcon, 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); + (*current_liboctave_warning_handler) + ("dgelsd: rank deficient %dx%d matrix, rank = %d", m, n, rank); 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); } @@ -2388,7 +2388,7 @@ FloatComplexMatrix FloatMatrix::lssolve (const FloatComplexMatrix& b, octave_idx_type& info, - octave_idx_type& rank) const + octave_idx_type& rank) const { FloatComplexMatrix tmp (*this); float rcon; @@ -2397,7 +2397,7 @@ FloatComplexMatrix FloatMatrix::lssolve (const FloatComplexMatrix& b, octave_idx_type& info, - octave_idx_type& rank, float& rcon) const + octave_idx_type& rank, float& rcon) const { FloatComplexMatrix tmp (*this); return tmp.lssolve (b, info, rank, rcon); @@ -2422,7 +2422,7 @@ FloatColumnVector FloatMatrix::lssolve (const FloatColumnVector& b, octave_idx_type& info, - octave_idx_type& rank) const + octave_idx_type& rank) const { float rcon; return lssolve (b, info, rank, rcon); @@ -2430,7 +2430,7 @@ FloatColumnVector FloatMatrix::lssolve (const FloatColumnVector& b, octave_idx_type& info, - octave_idx_type& rank, float &rcon) const + octave_idx_type& rank, float &rcon) const { FloatColumnVector retval; @@ -2451,14 +2451,14 @@ rcon = -1.0; if (m != n) - { - retval = FloatColumnVector (maxmn, 0.0); - - for (octave_idx_type i = 0; i < m; i++) - retval.elem (i) = b.elem (i); - } + { + retval = FloatColumnVector (maxmn, 0.0); + + for (octave_idx_type i = 0; i < m; i++) + retval.elem (i) = b.elem (i); + } else - retval = b; + retval = b; FloatMatrix atmp = *this; float *tmp_data = atmp.fortran_vec (); @@ -2474,10 +2474,10 @@ octave_idx_type smlsiz; F77_FUNC (xilaenv, XILAENV) (9, F77_CONST_CHAR_ARG2 ("SGELSD", 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 iwork because DGELSD in older versions // of LAPACK does not return it on a query call. @@ -2490,36 +2490,36 @@ #endif octave_idx_type nlvl = static_cast<octave_idx_type> (tmp) + 1; if (nlvl < 0) - nlvl = 0; + nlvl = 0; 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 (sgelsd, SGELSD, (m, n, nrhs, tmp_data, m, pretval, maxmn, - ps, rcon, rank, work.fortran_vec (), - lwork, piwork, info)); + ps, rcon, rank, work.fortran_vec (), + lwork, piwork, info)); lwork = static_cast<octave_idx_type> (work(0)); work.resize (lwork); F77_XFCN (sgelsd, SGELSD, (m, n, nrhs, tmp_data, m, pretval, - maxmn, ps, rcon, rank, - work.fortran_vec (), lwork, - piwork, info)); + maxmn, ps, rcon, rank, + work.fortran_vec (), lwork, + piwork, info)); 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) - rcon = 0.0; - else - rcon = s.elem (minmn - 1) / s.elem (0); - } + { + if (rank < minmn) + (*current_liboctave_warning_handler) + ("dgelsd: rank deficient %dx%d matrix, rank = %d", m, n, rank); + if (s.elem (0) == 0.0) + rcon = 0.0; + else + rcon = s.elem (minmn - 1) / s.elem (0); + } retval.resize (n, nrhs); } @@ -2548,7 +2548,7 @@ FloatComplexColumnVector FloatMatrix::lssolve (const FloatComplexColumnVector& b, octave_idx_type& info, - octave_idx_type& rank) const + octave_idx_type& rank) const { FloatComplexMatrix tmp (*this); float rcon; @@ -2557,7 +2557,7 @@ FloatComplexColumnVector FloatMatrix::lssolve (const FloatComplexColumnVector& b, octave_idx_type& info, - octave_idx_type& rank, float &rcon) const + octave_idx_type& rank, float &rcon) const { FloatComplexMatrix tmp (*this); return tmp.lssolve (b, info, rank, rcon); @@ -2628,13 +2628,13 @@ retval = FloatMatrix (len, a_len); float *c = retval.fortran_vec (); - + F77_XFCN (sgemm, SGEMM, (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; @@ -2650,14 +2650,14 @@ if (neg_zero) { for (octave_idx_type i = 0; i < nel; i++) - if (lo_ieee_signbit (elem (i))) - return true; + if (lo_ieee_signbit (elem (i))) + return true; } else { for (octave_idx_type i = 0; i < nel; i++) - if (elem (i) < 0) - return true; + if (elem (i) < 0) + return true; } return false; @@ -2672,7 +2672,7 @@ { float val = elem (i); if (xisnan (val)) - return true; + return true; } return false; @@ -2687,7 +2687,7 @@ { float val = elem (i); if (xisinf (val) || xisnan (val)) - return true; + return true; } return false; @@ -2702,7 +2702,7 @@ { float val = elem (i); if (val != 0 && val != 1) - return true; + return true; } return false; @@ -2717,9 +2717,9 @@ { float val = elem (i); if (xisnan (val) || D_NINT (val) == val) - continue; + continue; else - return false; + return false; } return true; @@ -2746,13 +2746,13 @@ float val = elem (i); if (val > max_val) - max_val = val; + max_val = val; if (val < min_val) - min_val = val; + min_val = val; if (D_NINT (val) != val) - return false; + return false; } return true; @@ -2768,8 +2768,8 @@ float val = elem (i); if (! (xisnan (val) || xisinf (val)) - && fabs (val) > FLT_MAX) - return true; + && fabs (val) > FLT_MAX) + return true; } return false; @@ -2855,33 +2855,33 @@ for (octave_idx_type i = 0; i < nr; i++) { - octave_idx_type idx_j; - - float tmp_min = octave_Float_NaN; - - for (idx_j = 0; idx_j < nc; idx_j++) - { - tmp_min = elem (i, idx_j); - - if (! xisnan (tmp_min)) - break; - } - - for (octave_idx_type j = idx_j+1; j < nc; j++) - { - float tmp = elem (i, j); - - if (xisnan (tmp)) - continue; - else if (tmp < tmp_min) - { - idx_j = j; - tmp_min = tmp; - } - } - - result.elem (i) = tmp_min; - idx_arg.elem (i) = xisnan (tmp_min) ? 0 : idx_j; + octave_idx_type idx_j; + + float tmp_min = octave_Float_NaN; + + for (idx_j = 0; idx_j < nc; idx_j++) + { + tmp_min = elem (i, idx_j); + + if (! xisnan (tmp_min)) + break; + } + + for (octave_idx_type j = idx_j+1; j < nc; j++) + { + float tmp = elem (i, j); + + if (xisnan (tmp)) + continue; + else if (tmp < tmp_min) + { + idx_j = j; + tmp_min = tmp; + } + } + + result.elem (i) = tmp_min; + idx_arg.elem (i) = xisnan (tmp_min) ? 0 : idx_j; } } @@ -2910,33 +2910,33 @@ for (octave_idx_type i = 0; i < nr; i++) { - octave_idx_type idx_j; - - float tmp_max = octave_Float_NaN; - - for (idx_j = 0; idx_j < nc; idx_j++) - { - tmp_max = elem (i, idx_j); - - if (! xisnan (tmp_max)) - break; - } - - for (octave_idx_type j = idx_j+1; j < nc; j++) - { - float tmp = elem (i, j); - - if (xisnan (tmp)) - continue; - else if (tmp > tmp_max) - { - idx_j = j; - tmp_max = tmp; - } - } - - result.elem (i) = tmp_max; - idx_arg.elem (i) = xisnan (tmp_max) ? 0 : idx_j; + octave_idx_type idx_j; + + float tmp_max = octave_Float_NaN; + + for (idx_j = 0; idx_j < nc; idx_j++) + { + tmp_max = elem (i, idx_j); + + if (! xisnan (tmp_max)) + break; + } + + for (octave_idx_type j = idx_j+1; j < nc; j++) + { + float tmp = elem (i, j); + + if (xisnan (tmp)) + continue; + else if (tmp > tmp_max) + { + idx_j = j; + tmp_max = tmp; + } + } + + result.elem (i) = tmp_max; + idx_arg.elem (i) = xisnan (tmp_max) ? 0 : idx_j; } } @@ -2965,33 +2965,33 @@ for (octave_idx_type j = 0; j < nc; j++) { - octave_idx_type idx_i; - - float tmp_min = octave_Float_NaN; - - for (idx_i = 0; idx_i < nr; idx_i++) - { - tmp_min = elem (idx_i, j); - - if (! xisnan (tmp_min)) - break; - } - - for (octave_idx_type i = idx_i+1; i < nr; i++) - { - float tmp = elem (i, j); - - if (xisnan (tmp)) - continue; - else if (tmp < tmp_min) - { - idx_i = i; - tmp_min = tmp; - } - } - - result.elem (j) = tmp_min; - idx_arg.elem (j) = xisnan (tmp_min) ? 0 : idx_i; + octave_idx_type idx_i; + + float tmp_min = octave_Float_NaN; + + for (idx_i = 0; idx_i < nr; idx_i++) + { + tmp_min = elem (idx_i, j); + + if (! xisnan (tmp_min)) + break; + } + + for (octave_idx_type i = idx_i+1; i < nr; i++) + { + float tmp = elem (i, j); + + if (xisnan (tmp)) + continue; + else if (tmp < tmp_min) + { + idx_i = i; + tmp_min = tmp; + } + } + + result.elem (j) = tmp_min; + idx_arg.elem (j) = xisnan (tmp_min) ? 0 : idx_i; } } @@ -3020,33 +3020,33 @@ for (octave_idx_type j = 0; j < nc; j++) { - octave_idx_type idx_i; - - float tmp_max = octave_Float_NaN; - - for (idx_i = 0; idx_i < nr; idx_i++) - { - tmp_max = elem (idx_i, j); - - if (! xisnan (tmp_max)) - break; - } - - for (octave_idx_type i = idx_i+1; i < nr; i++) - { - float tmp = elem (i, j); - - if (xisnan (tmp)) - continue; - else if (tmp > tmp_max) - { - idx_i = i; - tmp_max = tmp; - } - } - - result.elem (j) = tmp_max; - idx_arg.elem (j) = xisnan (tmp_max) ? 0 : idx_i; + octave_idx_type idx_i; + + float tmp_max = octave_Float_NaN; + + for (idx_i = 0; idx_i < nr; idx_i++) + { + tmp_max = elem (idx_i, j); + + if (! xisnan (tmp_max)) + break; + } + + for (octave_idx_type i = idx_i+1; i < nr; i++) + { + float tmp = elem (i, j); + + if (xisnan (tmp)) + continue; + else if (tmp > tmp_max) + { + idx_i = i; + tmp_max = tmp; + } + } + + result.elem (j) = tmp_max; + idx_arg.elem (j) = xisnan (tmp_max) ? 0 : idx_i; } } @@ -3059,10 +3059,10 @@ for (octave_idx_type i = 0; i < a.rows (); i++) { for (octave_idx_type j = 0; j < a.cols (); j++) - { - os << " "; - octave_write_float (os, a.elem (i, j)); - } + { + os << " "; + octave_write_float (os, a.elem (i, j)); + } os << "\n"; } return os; @@ -3078,14 +3078,14 @@ { float tmp; for (octave_idx_type i = 0; i < nr; i++) - for (octave_idx_type j = 0; j < nc; j++) - { - tmp = octave_read_value<float> (is); - if (is) - a.elem (i, j) = tmp; - else - goto done; - } + for (octave_idx_type j = 0; j < nc; j++) + { + tmp = octave_read_value<float> (is); + if (is) + a.elem (i, j) = tmp; + else + goto done; + } } done: @@ -3147,11 +3147,11 @@ float *px = cx.fortran_vec (); F77_XFCN (strsyl, STRSYL, (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? @@ -3209,13 +3209,13 @@ else { if (a_nr == 0 || a_nc == 0 || b_nc == 0) - retval = FloatMatrix (a_nr, b_nc, 0.0); + retval = FloatMatrix (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 = FloatMatrix (a_nr, b_nc); - float *c = retval.fortran_vec (); + float *c = retval.fortran_vec (); const char *ctra = get_blas_trans_arg (tra); F77_XFCN (ssyrk, SSYRK, (F77_CONST_CHAR_ARG2 ("U", 1), @@ -3230,25 +3230,25 @@ } else - { - octave_idx_type lda = a.rows (), tda = a.cols (); - octave_idx_type ldb = b.rows (), tdb = b.cols (); - - retval = FloatMatrix (a_nr, b_nc); - float *c = retval.fortran_vec (); - - if (b_nc == 1) - { - if (a_nr == 1) - F77_FUNC (xsdot, XSDOT) (a_nc, a.data (), 1, b.data (), 1, *c); - else - { + { + octave_idx_type lda = a.rows (), tda = a.cols (); + octave_idx_type ldb = b.rows (), tdb = b.cols (); + + retval = FloatMatrix (a_nr, b_nc); + float *c = retval.fortran_vec (); + + if (b_nc == 1) + { + if (a_nr == 1) + F77_FUNC (xsdot, XSDOT) (a_nc, a.data (), 1, b.data (), 1, *c); + else + { const char *ctra = get_blas_trans_arg (tra); - F77_XFCN (sgemv, SGEMV, (F77_CONST_CHAR_ARG2 (ctra, 1), - lda, tda, 1.0, a.data (), lda, - b.data (), 1, 0.0, c, 1 - F77_CHAR_ARG_LEN (1))); - } + F77_XFCN (sgemv, SGEMV, (F77_CONST_CHAR_ARG2 (ctra, 1), + lda, tda, 1.0, a.data (), lda, + b.data (), 1, 0.0, c, 1 + F77_CHAR_ARG_LEN (1))); + } } else if (a_nr == 1) { @@ -3258,18 +3258,18 @@ a.data (), 1, 0.0, c, 1 F77_CHAR_ARG_LEN (1))); } - else - { + else + { const char *ctra = get_blas_trans_arg (tra); const char *ctrb = get_blas_trans_arg (trb); - F77_XFCN (sgemm, SGEMM, (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 (sgemm, SGEMM, (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; @@ -3301,8 +3301,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 (d, m (i, j)); + octave_quit (); + result (i, j) = xmin (d, m (i, j)); } return result; @@ -3321,8 +3321,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), d); + octave_quit (); + result (i, j) = xmin (m (i, j), d); } return result; @@ -3337,7 +3337,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 FloatMatrix (); } @@ -3348,8 +3348,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 (a (i, j), b (i, j)); + octave_quit (); + result (i, j) = xmin (a (i, j), b (i, j)); } return result; @@ -3368,8 +3368,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 (d, m (i, j)); + octave_quit (); + result (i, j) = xmax (d, m (i, j)); } return result; @@ -3388,8 +3388,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), d); + octave_quit (); + result (i, j) = xmax (m (i, j), d); } return result; @@ -3404,7 +3404,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 FloatMatrix (); } @@ -3415,8 +3415,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 (a (i, j), b (i, j)); + octave_quit (); + result (i, j) = xmax (a (i, j), b (i, j)); } return result;