jwe/octave: liboctave/array/fCMatrix.cc comparison

comparison liboctave/array/fCMatrix.cc @ 21136:7cac4e7458f2

maint: clean up code around calls to current_liboctave_error_handler. Remove statements after call to handler that are no longer reachable. Place input validation first and immediately call handler if necessary. Change if/error_handler/else to if/error_handler and re-indent code. * Array-util.cc, Array.cc, CColVector.cc, CDiagMatrix.cc, CMatrix.cc, CNDArray.cc, CRowVector.cc, CSparse.cc, DiagArray2.cc, MArray.cc, PermMatrix.cc, Sparse.cc, Sparse.h, chMatrix.cc, chNDArray.cc, dColVector.cc, dDiagMatrix.cc, dMatrix.cc, dNDArray.cc, dRowVector.cc, dSparse.cc, fCColVector.cc, fCDiagMatrix.cc, fCMatrix.cc, fCNDArray.cc, fCRowVector.cc, fColVector.cc, fDiagMatrix.cc, fMatrix.cc, fNDArray.cc, fRowVector.cc, idx-vector.cc, CmplxAEPBAL.cc, CmplxCHOL.cc, CmplxGEPBAL.cc, CmplxHESS.cc, CmplxLU.cc, CmplxQR.cc, CmplxSCHUR.cc, CmplxSVD.cc, DASPK.cc, EIG.cc, LSODE.cc, Quad.cc, SparseCmplxCHOL.cc, SparseCmplxLU.cc, SparseCmplxQR.cc, SparseQR.cc, SparsedbleCHOL.cc, SparsedbleLU.cc, base-lu.cc, bsxfun-defs.cc, dbleAEPBAL.cc, dbleCHOL.cc, dbleGEPBAL.cc, dbleHESS.cc, dbleLU.cc, dbleQR.cc, dbleSCHUR.cc, dbleSVD.cc, eigs-base.cc, fCmplxAEPBAL.cc, fCmplxCHOL.cc, fCmplxLU.cc, fCmplxQR.cc, fCmplxSCHUR.cc, fEIG.cc, floatAEPBAL.cc, floatCHOL.cc, floatGEPBAL.cc, floatHESS.cc, floatLU.cc, floatQR.cc, floatSCHUR.cc, floatSVD.cc, lo-specfun.cc, oct-fftw.cc, oct-rand.cc, oct-spparms.cc, sparse-base-chol.cc, sparse-dmsolve.cc, file-ops.cc, lo-sysdep.cc, mach-info.cc, oct-env.cc, oct-syscalls.cc, cmd-edit.cc, cmd-hist.cc, data-conv.cc, lo-ieee.cc, lo-regexp.cc, oct-base64.cc, oct-shlib.cc, pathsearch.cc, singleton-cleanup.cc, sparse-util.cc, unwind-prot.cc: Remove statements after call to handler that are no longer reachable. Place input validation first and immediately call handler if necessary. Change if/error_handler/else to if/error_handler and re-indent code.

author	Rik <rik@octave.org>
date	Sat, 23 Jan 2016 13:52:03 -0800
parents	228b65504557
children	40051830f89b

comparison

equal deleted inserted replaced

-:95da3bc8a281
+:7cac4e7458f2
 {
 octave_idx_type a_nr = a.rows ();
 octave_idx_type a_nc = a.cols ();
 if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ())
-{
+(*current_liboctave_error_handler) ("range error for insert");
-(*current_liboctave_error_handler) ("range error for insert");
-return *this;
-}
 if (a_nr >0 && a_nc > 0)
 {
 make_unique ();
 octave_idx_type r, octave_idx_type c)
 {
 octave_idx_type a_len = a.numel ();
 if (r < 0 || r >= rows () || c < 0 || c + a_len > cols ())
-{
+(*current_liboctave_error_handler) ("range error for insert");
-(*current_liboctave_error_handler) ("range error for insert");
-return *this;
-}
 if (a_len > 0)
 {
 make_unique ();
 octave_idx_type r, octave_idx_type c)
 {
 octave_idx_type a_len = a.numel ();
 if (r < 0 || r + a_len > rows () || c < 0 || c >= cols ())
-{
+(*current_liboctave_error_handler) ("range error for insert");
-(*current_liboctave_error_handler) ("range error for insert");
-return *this;
-}
 if (a_len > 0)
 {
 make_unique ();
 {
 octave_idx_type a_nr = a.rows ();
 octave_idx_type a_nc = a.cols ();
 if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ())
-{
+(*current_liboctave_error_handler) ("range error for insert");
-(*current_liboctave_error_handler) ("range error for insert");
-return *this;
-}
 fill (0.0, r, c, r + a_nr - 1, c + a_nc - 1);
 octave_idx_type a_len = a.length ();
 FloatComplexMatrix::insert (const FloatComplexRowVector& a,
 octave_idx_type r, octave_idx_type c)
 {
 octave_idx_type a_len = a.numel ();
 if (r < 0 || r >= rows () || c < 0 || c + a_len > cols ())
-{
+(*current_liboctave_error_handler) ("range error for insert");
-(*current_liboctave_error_handler) ("range error for insert");
-return *this;
-}
 for (octave_idx_type i = 0; i < a_len; i++)
 elem (r, c+i) = a.elem (i);
 return *this;
 octave_idx_type r, octave_idx_type c)
 {
 octave_idx_type a_len = a.numel ();
 if (r < 0 || r + a_len > rows () || c < 0 || c >= cols ())
-{
+(*current_liboctave_error_handler) ("range error for insert");
-(*current_liboctave_error_handler) ("range error for insert");
-return *this;
-}
 if (a_len > 0)
 {
 make_unique ();
 {
 octave_idx_type a_nr = a.rows ();
 octave_idx_type a_nc = a.cols ();
 if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ())
-{
+(*current_liboctave_error_handler) ("range error for insert");
-(*current_liboctave_error_handler) ("range error for insert");
-return *this;
-}
 fill (0.0, r, c, r + a_nr - 1, c + a_nc - 1);
 octave_idx_type a_len = a.length ();
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0
 || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc)
-{
+(*current_liboctave_error_handler) ("range error for fill");
-(*current_liboctave_error_handler) ("range error for fill");
-return *this;
-}
 if (r1 > r2) { std::swap (r1, r2); }
 if (c1 > c2) { std::swap (c1, c2); }
 if (r2 >= r1 && c2 >= c1)
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0
 || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc)
-{
+(*current_liboctave_error_handler) ("range error for fill");
-(*current_liboctave_error_handler) ("range error for fill");
-return *this;
-}
 if (r1 > r2) { std::swap (r1, r2); }
 if (c1 > c2) { std::swap (c1, c2); }
 if (r2 >= r1 && c2 >=c1)
 FloatComplexMatrix::append (const FloatMatrix& a) const
 {
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 if (nr != a.rows ())
-{
+(*current_liboctave_error_handler) ("row dimension mismatch for append");
-(*current_liboctave_error_handler) ("row dimension mismatch for append");
-return *this;
-}
 octave_idx_type nc_insert = nc;
 FloatComplexMatrix retval (nr, nc + a.cols ());
 retval.insert (*this, 0, 0);
 retval.insert (a, 0, nc_insert);
 FloatComplexMatrix::append (const FloatRowVector& a) const
 {
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 if (nr != 1)
-{
+(*current_liboctave_error_handler) ("row dimension mismatch for append");
-(*current_liboctave_error_handler) ("row dimension mismatch for append");
-return *this;
-}
 octave_idx_type nc_insert = nc;
 FloatComplexMatrix retval (nr, nc + a.numel ());
 retval.insert (*this, 0, 0);
 retval.insert (a, 0, nc_insert);
 FloatComplexMatrix::append (const FloatColumnVector& a) const
 {
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 if (nr != a.numel ())
-{
+(*current_liboctave_error_handler) ("row dimension mismatch for append");
-(*current_liboctave_error_handler) ("row dimension mismatch for append");
-return *this;
-}
 octave_idx_type nc_insert = nc;
 FloatComplexMatrix retval (nr, nc + 1);
 retval.insert (*this, 0, 0);
 retval.insert (a, 0, nc_insert);
 FloatComplexMatrix::append (const FloatDiagMatrix& a) const
 {
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 if (nr != a.rows ())
-{
+(*current_liboctave_error_handler) ("row dimension mismatch for append");
-(*current_liboctave_error_handler) ("row dimension mismatch for append");
-return *this;
-}
 octave_idx_type nc_insert = nc;
 FloatComplexMatrix retval (nr, nc + a.cols ());
 retval.insert (*this, 0, 0);
 retval.insert (a, 0, nc_insert);
 FloatComplexMatrix::append (const FloatComplexMatrix& a) const
 {
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 if (nr != a.rows ())
-{
+(*current_liboctave_error_handler) ("row dimension mismatch for append");
-(*current_liboctave_error_handler) ("row dimension mismatch for append");
-return *this;
-}
 octave_idx_type nc_insert = nc;
 FloatComplexMatrix retval (nr, nc + a.cols ());
 retval.insert (*this, 0, 0);
 retval.insert (a, 0, nc_insert);
 FloatComplexMatrix::append (const FloatComplexRowVector& a) const
 {
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 if (nr != 1)
-{
+(*current_liboctave_error_handler) ("row dimension mismatch for append");
-(*current_liboctave_error_handler) ("row dimension mismatch for append");
-return *this;
-}
 octave_idx_type nc_insert = nc;
 FloatComplexMatrix retval (nr, nc + a.numel ());
 retval.insert (*this, 0, 0);
 retval.insert (a, 0, nc_insert);
 FloatComplexMatrix::append (const FloatComplexColumnVector& a) const
 {
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 if (nr != a.numel ())
-{
+(*current_liboctave_error_handler) ("row dimension mismatch for append");
-(*current_liboctave_error_handler) ("row dimension mismatch for append");
-return *this;
-}
 octave_idx_type nc_insert = nc;
 FloatComplexMatrix retval (nr, nc + 1);
 retval.insert (*this, 0, 0);
 retval.insert (a, 0, nc_insert);
 FloatComplexMatrix::append (const FloatComplexDiagMatrix& a) const
 {
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 if (nr != a.rows ())
-{
+(*current_liboctave_error_handler) ("row dimension mismatch for append");
-(*current_liboctave_error_handler) ("row dimension mismatch for append");
-return *this;
-}
 octave_idx_type nc_insert = nc;
 FloatComplexMatrix retval (nr, nc + a.cols ());
 retval.insert (*this, 0, 0);
 retval.insert (a, 0, nc_insert);
 FloatComplexMatrix::stack (const FloatMatrix& a) const
 {
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 if (nc != a.cols ())
-{
+(*current_liboctave_error_handler) ("column dimension mismatch for stack");
-(*current_liboctave_error_handler)
-("column dimension mismatch for stack");
-return *this;
-}
 octave_idx_type nr_insert = nr;
 FloatComplexMatrix retval (nr + a.rows (), nc);
 retval.insert (*this, 0, 0);
 retval.insert (a, nr_insert, 0);
 FloatComplexMatrix::stack (const FloatRowVector& a) const
 {
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 if (nc != a.numel ())
-{
+(*current_liboctave_error_handler) ("column dimension mismatch for stack");
-(*current_liboctave_error_handler)
-("column dimension mismatch for stack");
-return *this;
-}
 octave_idx_type nr_insert = nr;
 FloatComplexMatrix retval (nr + 1, nc);
 retval.insert (*this, 0, 0);
 retval.insert (a, nr_insert, 0);
 FloatComplexMatrix::stack (const FloatColumnVector& a) const
 {
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 if (nc != 1)
-{
+(*current_liboctave_error_handler) ("column dimension mismatch for stack");
-(*current_liboctave_error_handler)
-("column dimension mismatch for stack");
-return *this;
-}
 octave_idx_type nr_insert = nr;
 FloatComplexMatrix retval (nr + a.numel (), nc);
 retval.insert (*this, 0, 0);
 retval.insert (a, nr_insert, 0);
 FloatComplexMatrix::stack (const FloatDiagMatrix& a) const
 {
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 if (nc != a.cols ())
-{
+(*current_liboctave_error_handler) ("column dimension mismatch for stack");
-(*current_liboctave_error_handler)
-("column dimension mismatch for stack");
-return *this;
-}
 octave_idx_type nr_insert = nr;
 FloatComplexMatrix retval (nr + a.rows (), nc);
 retval.insert (*this, 0, 0);
 retval.insert (a, nr_insert, 0);
 FloatComplexMatrix::stack (const FloatComplexMatrix& a) const
 {
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 if (nc != a.cols ())
-{
+(*current_liboctave_error_handler) ("column dimension mismatch for stack");
-(*current_liboctave_error_handler)
-("column dimension mismatch for stack");
-return *this;
-}
 octave_idx_type nr_insert = nr;
 FloatComplexMatrix retval (nr + a.rows (), nc);
 retval.insert (*this, 0, 0);
 retval.insert (a, nr_insert, 0);
 FloatComplexMatrix::stack (const FloatComplexRowVector& a) const
 {
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 if (nc != a.numel ())
-{
+(*current_liboctave_error_handler) ("column dimension mismatch for stack");
-(*current_liboctave_error_handler)
-("column dimension mismatch for stack");
-return *this;
-}
 octave_idx_type nr_insert = nr;
 FloatComplexMatrix retval (nr + 1, nc);
 retval.insert (*this, 0, 0);
 retval.insert (a, nr_insert, 0);
 FloatComplexMatrix::stack (const FloatComplexColumnVector& a) const
 {
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 if (nc != 1)
-{
+(*current_liboctave_error_handler) ("column dimension mismatch for stack");
-(*current_liboctave_error_handler)
-("column dimension mismatch for stack");
-return *this;
-}
 octave_idx_type nr_insert = nr;
 FloatComplexMatrix retval (nr + a.numel (), nc);
 retval.insert (*this, 0, 0);
 retval.insert (a, nr_insert, 0);
 FloatComplexMatrix::stack (const FloatComplexDiagMatrix& a) const
 {
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 if (nc != a.cols ())
-{
+(*current_liboctave_error_handler) ("column dimension mismatch for stack");
-(*current_liboctave_error_handler)
-("column dimension mismatch for stack");
-return *this;
-}
 octave_idx_type nr_insert = nr;
 FloatComplexMatrix retval (nr + a.rows (), nc);
 retval.insert (*this, 0, 0);
 retval.insert (a, nr_insert, 0);
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 if (nr != nc || nr == 0 || nc == 0)
 (*current_liboctave_error_handler) ("inverse requires square matrix");
-else
-{
+int typ = mattype.type ();
-int typ = mattype.type ();
+char uplo = (typ == MatrixType::Lower ? 'L' : 'U');
-char uplo = (typ == MatrixType::Lower ? 'L' : 'U');
+char udiag = 'N';
-char udiag = 'N';
+retval = *this;
-retval = *this;
+FloatComplex *tmp_data = retval.fortran_vec ();
-FloatComplex *tmp_data = retval.fortran_vec ();
+F77_XFCN (ctrtri, CTRTRI, (F77_CONST_CHAR_ARG2 (&uplo, 1),
-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)));
+// Throw-away extra info LAPACK gives so as to not change output.
+rcon = 0.0;
+if (info != 0)
+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, info
+nr, tmp_data, nr, rcon,
+cwork, rwork, ztrcon_info
+F77_CHAR_ARG_LEN (1)
 F77_CHAR_ARG_LEN (1)
 F77_CHAR_ARG_LEN (1)));
-// Throw-away extra info LAPACK gives so as to not change output.
+if (ztrcon_info != 0)
-rcon = 0.0;
-if (info != 0)
 info = -1;
-else if (calc_cond)
+}
-{
-octave_idx_type ztrcon_info = 0;
+if (info == -1 && ! force)
-char job = '1';
+retval = *this; // Restore matrix contents.
-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.
-}
 return retval;
 }
 FloatComplexMatrix
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 if (nr != nc)
 (*current_liboctave_error_handler) ("inverse requires square matrix");
+Array<octave_idx_type> ipvt (dim_vector (nr, 1));
+octave_idx_type *pipvt = ipvt.fortran_vec ();
+retval = *this;
+FloatComplex *tmp_data = retval.fortran_vec ();
+Array<FloatComplex> z (dim_vector (1, 1));
+octave_idx_type lwork = -1;
+// Query the optimum work array size.
+F77_XFCN (cgetri, CGETRI, (nc, tmp_data, nr, pipvt,
+z.fortran_vec (), lwork, info));
+lwork = static_cast<octave_idx_type> (std::real (z(0)));
+lwork = (lwork <  2 *nc ? 2*nc : lwork);
+z.resize (dim_vector (lwork, 1));
+FloatComplex *pz = z.fortran_vec ();
+info = 0;
+// 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 ();
+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;
+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 (dim_vector (2 * nc, 1));
+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.
 else
 {
-Array<octave_idx_type> ipvt (dim_vector (nr, 1));
+octave_idx_type zgetri_info = 0;
-octave_idx_type *pipvt = ipvt.fortran_vec ();
-retval = *this;
-FloatComplex *tmp_data = retval.fortran_vec ();
-Array<FloatComplex> z (dim_vector (1, 1));
-octave_idx_type lwork = -1;
-// Query the optimum work array size.
 F77_XFCN (cgetri, CGETRI, (nc, tmp_data, nr, pipvt,
-z.fortran_vec (), lwork, info));
+pz, lwork, zgetri_info));
-lwork = static_cast<octave_idx_type> (std::real (z(0)));
+if (zgetri_info != 0)
-lwork = (lwork <  2 *nc ? 2*nc : lwork);
-z.resize (dim_vector (lwork, 1));
-FloatComplex *pz = z.fortran_vec ();
-info = 0;
-// 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 ();
-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;
-else if (calc_cond)
+}
-{
-// Now calculate the condition number for non-singular matrix.
+if (info != 0)
-octave_idx_type zgecon_info = 0;
+mattype.mark_as_rectangular ();
-char job = '1';
-Array<float> rz (dim_vector (2 * nc, 1));
-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.
-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;
-}
-if (info != 0)
-mattype.mark_as_rectangular ();
-}
 return retval;
 }
 FloatComplexMatrix
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 if (nr != nc)
 (*current_liboctave_error_handler) ("matrix must be square");
-else
-{
+volatile int typ = mattype.type ();
-volatile int typ = mattype.type ();
+// Even though the matrix is marked as singular (Rectangular), we may
-// Even though the matrix is marked as singular (Rectangular), we may
+// still get a useful number from the LU factorization, because it always
-// still get a useful number from the LU factorization, because it always
+// completes.
-// completes.
+if (typ == MatrixType::Unknown)
-if (typ == MatrixType::Unknown)
+typ = mattype.type (*this);
-typ = mattype.type (*this);
+else if (typ == MatrixType::Rectangular)
-else if (typ == MatrixType::Rectangular)
+typ = MatrixType::Full;
-typ = MatrixType::Full;
+if (typ == MatrixType::Lower || typ == MatrixType::Upper)
-if (typ == MatrixType::Lower || typ == MatrixType::Upper)
+{
+for (octave_idx_type i = 0; i < nc; i++)
+retval *= elem (i,i);
+}
+else if (typ == MatrixType::Hermitian)
+{
+FloatComplexMatrix atmp = *this;
+FloatComplex *tmp_data = atmp.fortran_vec ();
+float anorm = 0;
+if (calc_cond) anorm = xnorm (*this, 1);
+char job = 'L';
+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 (dim_vector (2 * nc, 1));
+FloatComplex *pz = z.fortran_vec ();
+Array<float> rz (dim_vector (nc, 1));
+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;
 for (octave_idx_type i = 0; i < nc; i++)
-retval *= elem (i,i);
+retval *= atmp (i,i);
+retval = retval.square ();
 }
-else if (typ == MatrixType::Hermitian)
+}
+else if (typ != MatrixType::Full)
+(*current_liboctave_error_handler) ("det: invalid dense matrix type");
+if (typ == MatrixType::Full)
+{
+Array<octave_idx_type> ipvt (dim_vector (nr, 1));
+octave_idx_type *pipvt = ipvt.fortran_vec ();
+FloatComplexMatrix atmp = *this;
+FloatComplex *tmp_data = atmp.fortran_vec ();
+info = 0;
+// Calculate the norm of the matrix, for later use.
+float anorm = 0;
+if (calc_cond) anorm = xnorm (*this, 1);
+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)
 {
-FloatComplexMatrix atmp = *this;
+info = -1;
-FloatComplex *tmp_data = atmp.fortran_vec ();
+retval = FloatComplexDET ();
+}
-float anorm = 0;
+else
-if (calc_cond) anorm = xnorm (*this, 1);
+{
+if (calc_cond)
-char job = 'L';
-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;
+// Now calc the condition number for non-singular matrix.
-mattype.mark_as_unsymmetric ();
+char job = '1';
-typ = MatrixType::Full;
-}
-else
-{
 Array<FloatComplex> z (dim_vector (2 * nc, 1));
 FloatComplex *pz = z.fortran_vec ();
-Array<float> rz (dim_vector (nc, 1));
+Array<float> rz (dim_vector (2 * nc, 1));
 float *prz = rz.fortran_vec ();
-F77_XFCN (cpocon, CPOCON, (F77_CONST_CHAR_ARG2 (&job, 1),
+F77_XFCN (cgecon, CGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
-nr, tmp_data, nr, anorm,
+nc, tmp_data, nr, anorm,
 rcon, pz, prz, info
 F77_CHAR_ARG_LEN (1)));
-if (info != 0)
-rcon = 0.0;
-for (octave_idx_type i = 0; i < nc; i++)
-retval *= atmp (i,i);
-retval = retval.square ();
 }
-}
-else if (typ != MatrixType::Full)
-(*current_liboctave_error_handler) ("det: invalid dense matrix type");
-if (typ == MatrixType::Full)
-{
-Array<octave_idx_type> ipvt (dim_vector (nr, 1));
-octave_idx_type *pipvt = ipvt.fortran_vec ();
-FloatComplexMatrix atmp = *this;
-FloatComplex *tmp_data = atmp.fortran_vec ();
-info = 0;
-// Calculate the norm of the matrix, for later use.
-float anorm = 0;
-if (calc_cond) anorm = xnorm (*this, 1);
-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 = -1;
 retval = FloatComplexDET ();
 }
 else
 {
-if (calc_cond)
+for (octave_idx_type i = 0; i < nc; i++)
 {
-// Now calc the condition number for non-singular matrix.
+FloatComplex c = atmp(i,i);
-char job = '1';
+retval *= (ipvt(i) != (i+1)) ? -c : c;
-Array<FloatComplex> z (dim_vector (2 * nc, 1));
-FloatComplex *pz = z.fortran_vec ();
-Array<float> rz (dim_vector (2 * nc, 1));
-float *prz = rz.fortran_vec ();
-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 = -1;
-retval = FloatComplexDET ();
-}
-else
-{
-for (octave_idx_type i = 0; i < nc; i++)
-{
-FloatComplex c = atmp(i,i);
-retval *= (ipvt(i) != (i+1)) ? -c : c;
-}
 }
 }
 }
 }
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 if (nr != nc)
 (*current_liboctave_error_handler) ("matrix must be square");
-else if (nr == 0 || nc == 0)
+if (nr == 0 || nc == 0)
 rcon = octave_Inf;
 else
 {
 volatile int typ = mattype.type ();
 octave_idx_type nc = cols ();
 if (nr != b.rows ())
 (*current_liboctave_error_handler)
 ("matrix dimension mismatch solution of linear equations");
-else if (nr == 0 || nc == 0 || b.cols () == 0)
+if (nr == 0 || nc == 0 || b.cols () == 0)
 retval = FloatComplexMatrix (nc, b.cols (), FloatComplex (0.0, 0.0));
 else
 {
 volatile int typ = mattype.type ();
 octave_idx_type b_nc = b.cols ();
 rcon = 1.;
 info = 0;
 if (typ == MatrixType::Permuted_Upper)
-{
+(*current_liboctave_error_handler)
-(*current_liboctave_error_handler)
+("permuted triangular matrix not implemented");
-("permuted triangular matrix not implemented");
-}
 else
 {
 const FloatComplex *tmp_data = fortran_vec ();
 retval = b;
 octave_idx_type nc = cols ();
 if (nr != b.rows ())
 (*current_liboctave_error_handler)
 ("matrix dimension mismatch solution of linear equations");
-else if (nr == 0 || nc == 0 || b.cols () == 0)
+if (nr == 0 || nc == 0 || b.cols () == 0)
 retval = FloatComplexMatrix (nc, b.cols (), FloatComplex (0.0, 0.0));
 else
 {
 volatile int typ = mattype.type ();
 octave_idx_type b_nc = b.cols ();
 rcon = 1.;
 info = 0;
 if (typ == MatrixType::Permuted_Lower)
-{
+(*current_liboctave_error_handler)
-(*current_liboctave_error_handler)
+("permuted triangular matrix not implemented");
-("permuted triangular matrix not implemented");
-}
 else
 {
 const FloatComplex *tmp_data = fortran_vec ();
 retval = b;
 if (nr != nc || nr != b.rows ())
 (*current_liboctave_error_handler)
 ("matrix dimension mismatch solution of linear equations");
-else if (nr == 0 || b.cols () == 0)
+if (nr == 0 || b.cols () == 0)
 retval = FloatComplexMatrix (nc, b.cols (), FloatComplex (0.0, 0.0));
 else
 {
 volatile int typ = mattype.type ();
 retval = hermitian ().solve (mattype, b, info, rcon, sing_handler,
 singular_fallback);
 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
 retval = fsolve (mattype, b, info, rcon, sing_handler, true);
 else if (typ != MatrixType::Rectangular)
-{
+(*current_liboctave_error_handler) ("unknown matrix type");
-(*current_liboctave_error_handler) ("unknown matrix type");
-return FloatComplexMatrix ();
-}
 // Rectangular or one of the above solvers flags a singular matrix
 if (singular_fallback && mattype.type () == MatrixType::Rectangular)
 {
 octave_idx_type rank;
 octave_idx_type n = cols ();
 if (m != b.rows ())
 (*current_liboctave_error_handler)
 ("matrix dimension mismatch solution of linear equations");
-else if (m== 0 || n == 0 || b.cols () == 0)
+if (m== 0 || n == 0 || b.cols () == 0)
 retval = FloatComplexMatrix (n, b.cols (), FloatComplex (0.0, 0.0));
 else
 {
 volatile octave_idx_type minmn = (m < n ? m : n);
 octave_idx_type maxmn = m > n ? m : n;
 octave_idx_type n = cols ();
 if (m != b.numel ())
 (*current_liboctave_error_handler)
 ("matrix dimension mismatch solution of linear equations");
-else if (m == 0 || n == 0 || b.cols () == 0)
+if (m == 0 || n == 0 || b.cols () == 0)
 retval = FloatComplexColumnVector (n, FloatComplex (0.0, 0.0));
 else
 {
 volatile octave_idx_type minmn = (m < n ? m : n);
 octave_idx_type maxmn = m > n ? m : n;
 octave_idx_type nc = cols ();
 if (nr == 1 || nc == 1)
 retval = FloatComplexDiagMatrix (*this, m, n);
 else
-(*current_liboctave_error_handler)
+(*current_liboctave_error_handler) ("diag: expecting vector argument");
-("diag: expecting vector argument");
 return retval;
 }
 bool
 {
 octave_idx_type nr = a.rows ();
 octave_idx_type nc = a.columns ();
 if (nr != b.rows () || nc != b.columns ())
-{
+(*current_liboctave_error_handler)
-(*current_liboctave_error_handler)
+("two-arg min requires same size arguments");
-("two-arg min requires same size arguments");
-return FloatComplexMatrix ();
-}
 EMPTY_RETURN_CHECK (FloatComplexMatrix);
 FloatComplexMatrix result (nr, nc);
 {
 octave_idx_type nr = a.rows ();
 octave_idx_type nc = a.columns ();
 if (nr != b.rows () || nc != b.columns ())
-{
+(*current_liboctave_error_handler)
-(*current_liboctave_error_handler)
+("two-arg max requires same size arguments");
-("two-arg max requires same size arguments");
-return FloatComplexMatrix ();
-}
 EMPTY_RETURN_CHECK (FloatComplexMatrix);
 FloatComplexMatrix result (nr, nc);

Mercurial > jwe > octave

comparison liboctave/array/fCMatrix.cc @ 21136:7cac4e7458f2