Mercurial > octave-nkf
diff liboctave/eigs-base.cc @ 10314:07ebe522dac2
untabify liboctave C++ sources
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Thu, 11 Feb 2010 12:23:32 -0500 |
| parents | 4c0cdbe0acca |
| children | 12884915a8e4 |
line wrap: on
line diff
--- a/liboctave/eigs-base.cc Thu Feb 11 12:16:43 2010 -0500 +++ b/liboctave/eigs-base.cc Thu Feb 11 12:23:32 2010 -0500 @@ -54,82 +54,82 @@ { F77_RET_T F77_FUNC (dsaupd, DSAUPD) (octave_idx_type&, F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, const double&, - double*, const octave_idx_type&, double*, - const octave_idx_type&, octave_idx_type*, - octave_idx_type*, double*, double*, - const octave_idx_type&, octave_idx_type& - F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL); + const octave_idx_type&, F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, const double&, + double*, const octave_idx_type&, double*, + const octave_idx_type&, octave_idx_type*, + octave_idx_type*, double*, double*, + const octave_idx_type&, octave_idx_type& + F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL); F77_RET_T F77_FUNC (dseupd, DSEUPD) (const int&, F77_CONST_CHAR_ARG_DECL, - octave_idx_type*, double*, double*, - const octave_idx_type&, const double&, - F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, - F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, - const double&, double*, const octave_idx_type&, - double*, const octave_idx_type&, octave_idx_type*, - octave_idx_type*, double*, double*, - const octave_idx_type&, octave_idx_type& - F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL - F77_CHAR_ARG_LEN_DECL); + octave_idx_type*, double*, double*, + const octave_idx_type&, const double&, + F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, + F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, + const double&, double*, const octave_idx_type&, + double*, const octave_idx_type&, octave_idx_type*, + octave_idx_type*, double*, double*, + 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 (dnaupd, DNAUPD) (octave_idx_type&, F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, F77_CONST_CHAR_ARG_DECL, - octave_idx_type&, const double&, - double*, const octave_idx_type&, double*, - const octave_idx_type&, octave_idx_type*, - octave_idx_type*, double*, double*, - const octave_idx_type&, octave_idx_type& - F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL); + const octave_idx_type&, F77_CONST_CHAR_ARG_DECL, + octave_idx_type&, const double&, + double*, const octave_idx_type&, double*, + const octave_idx_type&, octave_idx_type*, + octave_idx_type*, double*, double*, + const octave_idx_type&, octave_idx_type& + F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL); F77_RET_T F77_FUNC (dneupd, DNEUPD) (const int&, F77_CONST_CHAR_ARG_DECL, - octave_idx_type*, double*, double*, - double*, const octave_idx_type&, const double&, - const double&, double*, F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, F77_CONST_CHAR_ARG_DECL, - octave_idx_type&, const double&, double*, - const octave_idx_type&, double*, - const octave_idx_type&, octave_idx_type*, - octave_idx_type*, double*, double*, - const octave_idx_type&, octave_idx_type& - F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL - F77_CHAR_ARG_LEN_DECL); + octave_idx_type*, double*, double*, + double*, const octave_idx_type&, const double&, + const double&, double*, F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, F77_CONST_CHAR_ARG_DECL, + octave_idx_type&, const double&, double*, + const octave_idx_type&, double*, + const octave_idx_type&, octave_idx_type*, + octave_idx_type*, double*, double*, + 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 (znaupd, ZNAUPD) (octave_idx_type&, F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, const double&, - Complex*, const octave_idx_type&, Complex*, - const octave_idx_type&, octave_idx_type*, - octave_idx_type*, Complex*, Complex*, - const octave_idx_type&, double *, octave_idx_type& - F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL); + const octave_idx_type&, F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, const double&, + Complex*, const octave_idx_type&, Complex*, + const octave_idx_type&, octave_idx_type*, + octave_idx_type*, Complex*, Complex*, + const octave_idx_type&, double *, octave_idx_type& + F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL); F77_RET_T F77_FUNC (zneupd, ZNEUPD) (const int&, F77_CONST_CHAR_ARG_DECL, - octave_idx_type*, Complex*, Complex*, - const octave_idx_type&, const Complex&, - Complex*, F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, const double&, - Complex*, const octave_idx_type&, Complex*, - const octave_idx_type&, octave_idx_type*, - octave_idx_type*, Complex*, Complex*, - const octave_idx_type&, double *, octave_idx_type& - F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL - F77_CHAR_ARG_LEN_DECL); + octave_idx_type*, Complex*, Complex*, + const octave_idx_type&, const Complex&, + Complex*, F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, F77_CONST_CHAR_ARG_DECL, + const octave_idx_type&, const double&, + Complex*, const octave_idx_type&, Complex*, + const octave_idx_type&, octave_idx_type*, + octave_idx_type*, Complex*, Complex*, + const octave_idx_type&, double *, octave_idx_type& + F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL); F77_RET_T F77_FUNC (dgemv, DGEMV) (F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, const octave_idx_type&, const double&, - const double*, const octave_idx_type&, const double*, - const octave_idx_type&, const double&, double*, - const octave_idx_type& - F77_CHAR_ARG_LEN_DECL); + const octave_idx_type&, const octave_idx_type&, const double&, + const double*, const octave_idx_type&, const double*, + const octave_idx_type&, const double&, double*, + const octave_idx_type& + F77_CHAR_ARG_LEN_DECL); F77_RET_T @@ -148,7 +148,7 @@ static octave_idx_type lusolve (const SparseComplexMatrix&, const SparseComplexMatrix&, - ComplexMatrix&); + ComplexMatrix&); static octave_idx_type lusolve (const Matrix&, const Matrix&, Matrix&); @@ -158,7 +158,7 @@ static ComplexMatrix ltsolve (const SparseComplexMatrix&, const ColumnVector&, - const ComplexMatrix&); + const ComplexMatrix&); static Matrix ltsolve (const SparseMatrix&, const ColumnVector&, const Matrix&,); @@ -218,11 +218,11 @@ { retval.resize (n, b_nc); for (octave_idx_type j = 0; j < b_nc; j++) - { - for (octave_idx_type i = 0; i < n; i++) - retval.elem(static_cast<octave_idx_type>(qv[i]), j) = - tmp.elem(i,j); - } + { + for (octave_idx_type i = 0; i < n; i++) + retval.elem(static_cast<octave_idx_type>(qv[i]), j) = + tmp.elem(i,j); + } } return retval; @@ -243,7 +243,7 @@ for (octave_idx_type j = 0; j < b_nc; j++) { for (octave_idx_type i = 0; i < n; i++) - retval.elem(i,j) = m.elem(static_cast<octave_idx_type>(qv[i]), j); + retval.elem(i,j) = m.elem(static_cast<octave_idx_type>(qv[i]), j); } return U.solve (utyp, retval, err, rcond, 0); } @@ -270,14 +270,14 @@ octave_idx_type nc = m.cols (); F77_XFCN (dgemv, DGEMV, (F77_CONST_CHAR_ARG2 ("N", 1), - nr, nc, 1.0, m.data (), nr, - x, 1, 0.0, y, 1 - F77_CHAR_ARG_LEN (1))); + nr, nc, 1.0, m.data (), nr, + x, 1, 0.0, y, 1 + F77_CHAR_ARG_LEN (1))); if (f77_exception_encountered) { (*current_liboctave_error_handler) - ("eigs: unrecoverable error in dgemv"); + ("eigs: unrecoverable error in dgemv"); return false; } else @@ -286,7 +286,7 @@ static bool vector_product (const SparseComplexMatrix& m, const Complex* x, - Complex* y) + Complex* y) { octave_idx_type nc = m.cols (); @@ -307,14 +307,14 @@ octave_idx_type nc = m.cols (); F77_XFCN (zgemv, ZGEMV, (F77_CONST_CHAR_ARG2 ("N", 1), - nr, nc, 1.0, m.data (), nr, - x, 1, 0.0, y, 1 - F77_CHAR_ARG_LEN (1))); + nr, nc, 1.0, m.data (), nr, + x, 1, 0.0, y, 1 + F77_CHAR_ARG_LEN (1))); if (f77_exception_encountered) { (*current_liboctave_error_handler) - ("eigs: unrecoverable error in zgemv"); + ("eigs: unrecoverable error in zgemv"); return false; } else @@ -336,7 +336,7 @@ b = bt.transpose(); permB = ColumnVector(n); for (octave_idx_type i = 0; i < n; i++) - permB(i) = i; + permB(i) = i; return true; } } @@ -373,14 +373,14 @@ b = bt.hermitian(); permB = ColumnVector(n); for (octave_idx_type i = 0; i < n; i++) - permB(i) = i; + permB(i) = i; return true; } } static bool make_cholb (SparseComplexMatrix& b, SparseComplexMatrix& bt, - ColumnVector& permB) + ColumnVector& permB) { octave_idx_type info; SparseComplexCHOL fact (b, info, false); @@ -398,9 +398,9 @@ static bool LuAminusSigmaB (const SparseMatrix &m, const SparseMatrix &b, - bool cholB, const ColumnVector& permB, double sigma, - SparseMatrix &L, SparseMatrix &U, octave_idx_type *P, - octave_idx_type *Q) + bool cholB, const ColumnVector& permB, double sigma, + SparseMatrix &L, SparseMatrix &U, octave_idx_type *P, + octave_idx_type *Q) { bool have_b = ! b.is_empty (); octave_idx_type n = m.rows(); @@ -411,44 +411,44 @@ if (have_b) { if (cholB) - { - if (permB.length()) - { - SparseMatrix tmp(n,n,n); - for (octave_idx_type i = 0; i < n; i++) - { - tmp.xcidx(i) = i; - tmp.xridx(i) = - static_cast<octave_idx_type>(permB(i)); - tmp.xdata(i) = 1; - } - tmp.xcidx(n) = n; - - AminusSigmaB = AminusSigmaB - sigma * tmp * - b.transpose() * b * tmp.transpose(); - } - else - AminusSigmaB = AminusSigmaB - sigma * - b.transpose() * b; - } + { + if (permB.length()) + { + SparseMatrix tmp(n,n,n); + for (octave_idx_type i = 0; i < n; i++) + { + tmp.xcidx(i) = i; + tmp.xridx(i) = + static_cast<octave_idx_type>(permB(i)); + tmp.xdata(i) = 1; + } + tmp.xcidx(n) = n; + + AminusSigmaB = AminusSigmaB - sigma * tmp * + b.transpose() * b * tmp.transpose(); + } + else + AminusSigmaB = AminusSigmaB - sigma * + b.transpose() * b; + } else - AminusSigmaB = AminusSigmaB - sigma * b; + AminusSigmaB = AminusSigmaB - sigma * b; } else { SparseMatrix sigmat (n, n, n); - // Create sigma * speye(n,n) - sigmat.xcidx (0) = 0; - for (octave_idx_type i = 0; i < n; i++) - { - sigmat.xdata(i) = sigma; - sigmat.xridx(i) = i; - sigmat.xcidx(i+1) = i + 1; - } - - AminusSigmaB = AminusSigmaB - sigmat; - } + // Create sigma * speye(n,n) + sigmat.xcidx (0) = 0; + for (octave_idx_type i = 0; i < n; i++) + { + sigmat.xdata(i) = sigma; + sigmat.xridx(i) = i; + sigmat.xcidx(i+1) = i + 1; + } + + AminusSigmaB = AminusSigmaB - sigmat; + } SparseLU fact (AminusSigmaB); @@ -470,14 +470,14 @@ { double d = 0.; if (U.xcidx(j+1) > U.xcidx(j) && - U.xridx (U.xcidx(j+1)-1) == j) - d = std::abs (U.xdata (U.xcidx(j+1)-1)); + U.xridx (U.xcidx(j+1)-1) == j) + d = std::abs (U.xdata (U.xcidx(j+1)-1)); if (xisnan (minU) || d < minU) - minU = d; + minU = d; if (xisnan (maxU) || d > maxU) - maxU = d; + maxU = d; } double rcond = (minU / maxU); @@ -486,9 +486,9 @@ if (rcond_plus_one == 1.0 || xisnan (rcond)) { (*current_liboctave_warning_handler) - ("eigs: 'A - sigma*B' is singular, indicating sigma is exactly"); + ("eigs: 'A - sigma*B' is singular, indicating sigma is exactly"); (*current_liboctave_warning_handler) - (" an eigenvalue. Convergence is not guaranteed"); + (" an eigenvalue. Convergence is not guaranteed"); } return true; @@ -496,9 +496,9 @@ static bool LuAminusSigmaB (const Matrix &m, const Matrix &b, - bool cholB, const ColumnVector& permB, double sigma, - Matrix &L, Matrix &U, octave_idx_type *P, - octave_idx_type *Q) + bool cholB, const ColumnVector& permB, double sigma, + Matrix &L, Matrix &U, octave_idx_type *P, + octave_idx_type *Q) { bool have_b = ! b.is_empty (); octave_idx_type n = m.cols(); @@ -509,32 +509,32 @@ if (have_b) { if (cholB) - { - Matrix tmp = sigma * b.transpose() * b; - const double *pB = permB.fortran_vec(); - double *p = AminusSigmaB.fortran_vec(); - - if (permB.length()) - { - for (octave_idx_type j = 0; - j < b.cols(); j++) - for (octave_idx_type i = 0; - i < b.rows(); i++) - *p++ -= tmp.xelem (static_cast<octave_idx_type>(pB[i]), - static_cast<octave_idx_type>(pB[j])); - } - else - AminusSigmaB = AminusSigmaB - tmp; - } + { + Matrix tmp = sigma * b.transpose() * b; + const double *pB = permB.fortran_vec(); + double *p = AminusSigmaB.fortran_vec(); + + if (permB.length()) + { + for (octave_idx_type j = 0; + j < b.cols(); j++) + for (octave_idx_type i = 0; + i < b.rows(); i++) + *p++ -= tmp.xelem (static_cast<octave_idx_type>(pB[i]), + static_cast<octave_idx_type>(pB[j])); + } + else + AminusSigmaB = AminusSigmaB - tmp; + } else - AminusSigmaB = AminusSigmaB - sigma * b; + AminusSigmaB = AminusSigmaB - sigma * b; } else { double *p = AminusSigmaB.fortran_vec(); for (octave_idx_type i = 0; i < n; i++) - p[i*(n+1)] -= sigma; + p[i*(n+1)] -= sigma; } LU fact (AminusSigmaB); @@ -551,10 +551,10 @@ { double d = std::abs (U.xelem(j,j)); if (xisnan (minU) || d < minU) - minU = d; + minU = d; if (xisnan (maxU) || d > maxU) - maxU = d; + maxU = d; } double rcond = (minU / maxU); @@ -563,9 +563,9 @@ if (rcond_plus_one == 1.0 || xisnan (rcond)) { (*current_liboctave_warning_handler) - ("eigs: 'A - sigma*B' is singular, indicating sigma is exactly"); + ("eigs: 'A - sigma*B' is singular, indicating sigma is exactly"); (*current_liboctave_warning_handler) - (" an eigenvalue. Convergence is not guaranteed"); + (" an eigenvalue. Convergence is not guaranteed"); } return true; @@ -573,9 +573,9 @@ static bool LuAminusSigmaB (const SparseComplexMatrix &m, const SparseComplexMatrix &b, - bool cholB, const ColumnVector& permB, Complex sigma, - SparseComplexMatrix &L, SparseComplexMatrix &U, - octave_idx_type *P, octave_idx_type *Q) + bool cholB, const ColumnVector& permB, Complex sigma, + SparseComplexMatrix &L, SparseComplexMatrix &U, + octave_idx_type *P, octave_idx_type *Q) { bool have_b = ! b.is_empty (); octave_idx_type n = m.rows(); @@ -586,27 +586,27 @@ if (have_b) { if (cholB) - { - if (permB.length()) - { - SparseMatrix tmp(n,n,n); - for (octave_idx_type i = 0; i < n; i++) - { - tmp.xcidx(i) = i; - tmp.xridx(i) = - static_cast<octave_idx_type>(permB(i)); - tmp.xdata(i) = 1; - } - tmp.xcidx(n) = n; - - AminusSigmaB = AminusSigmaB - tmp * b.hermitian() * b * - tmp.transpose() * sigma; - } - else - AminusSigmaB = AminusSigmaB - sigma * b.hermitian() * b; - } + { + if (permB.length()) + { + SparseMatrix tmp(n,n,n); + for (octave_idx_type i = 0; i < n; i++) + { + tmp.xcidx(i) = i; + tmp.xridx(i) = + static_cast<octave_idx_type>(permB(i)); + tmp.xdata(i) = 1; + } + tmp.xcidx(n) = n; + + AminusSigmaB = AminusSigmaB - tmp * b.hermitian() * b * + tmp.transpose() * sigma; + } + else + AminusSigmaB = AminusSigmaB - sigma * b.hermitian() * b; + } else - AminusSigmaB = AminusSigmaB - sigma * b; + AminusSigmaB = AminusSigmaB - sigma * b; } else { @@ -615,11 +615,11 @@ // Create sigma * speye(n,n) sigmat.xcidx (0) = 0; for (octave_idx_type i = 0; i < n; i++) - { - sigmat.xdata(i) = sigma; - sigmat.xridx(i) = i; - sigmat.xcidx(i+1) = i + 1; - } + { + sigmat.xdata(i) = sigma; + sigmat.xridx(i) = i; + sigmat.xcidx(i+1) = i + 1; + } AminusSigmaB = AminusSigmaB - sigmat; } @@ -644,14 +644,14 @@ { double d = 0.; if (U.xcidx(j+1) > U.xcidx(j) && - U.xridx (U.xcidx(j+1)-1) == j) - d = std::abs (U.xdata (U.xcidx(j+1)-1)); + U.xridx (U.xcidx(j+1)-1) == j) + d = std::abs (U.xdata (U.xcidx(j+1)-1)); if (xisnan (minU) || d < minU) - minU = d; + minU = d; if (xisnan (maxU) || d > maxU) - maxU = d; + maxU = d; } double rcond = (minU / maxU); @@ -660,9 +660,9 @@ if (rcond_plus_one == 1.0 || xisnan (rcond)) { (*current_liboctave_warning_handler) - ("eigs: 'A - sigma*B' is singular, indicating sigma is exactly"); + ("eigs: 'A - sigma*B' is singular, indicating sigma is exactly"); (*current_liboctave_warning_handler) - (" an eigenvalue. Convergence is not guaranteed"); + (" an eigenvalue. Convergence is not guaranteed"); } return true; @@ -670,9 +670,9 @@ static bool LuAminusSigmaB (const ComplexMatrix &m, const ComplexMatrix &b, - bool cholB, const ColumnVector& permB, Complex sigma, - ComplexMatrix &L, ComplexMatrix &U, octave_idx_type *P, - octave_idx_type *Q) + bool cholB, const ColumnVector& permB, Complex sigma, + ComplexMatrix &L, ComplexMatrix &U, octave_idx_type *P, + octave_idx_type *Q) { bool have_b = ! b.is_empty (); octave_idx_type n = m.cols(); @@ -683,32 +683,32 @@ if (have_b) { if (cholB) - { - ComplexMatrix tmp = sigma * b.hermitian() * b; - const double *pB = permB.fortran_vec(); - Complex *p = AminusSigmaB.fortran_vec(); - - if (permB.length()) - { - for (octave_idx_type j = 0; - j < b.cols(); j++) - for (octave_idx_type i = 0; - i < b.rows(); i++) - *p++ -= tmp.xelem (static_cast<octave_idx_type>(pB[i]), - static_cast<octave_idx_type>(pB[j])); - } - else - AminusSigmaB = AminusSigmaB - tmp; - } + { + ComplexMatrix tmp = sigma * b.hermitian() * b; + const double *pB = permB.fortran_vec(); + Complex *p = AminusSigmaB.fortran_vec(); + + if (permB.length()) + { + for (octave_idx_type j = 0; + j < b.cols(); j++) + for (octave_idx_type i = 0; + i < b.rows(); i++) + *p++ -= tmp.xelem (static_cast<octave_idx_type>(pB[i]), + static_cast<octave_idx_type>(pB[j])); + } + else + AminusSigmaB = AminusSigmaB - tmp; + } else - AminusSigmaB = AminusSigmaB - sigma * b; + AminusSigmaB = AminusSigmaB - sigma * b; } else { Complex *p = AminusSigmaB.fortran_vec(); for (octave_idx_type i = 0; i < n; i++) - p[i*(n+1)] -= sigma; + p[i*(n+1)] -= sigma; } ComplexLU fact (AminusSigmaB); @@ -725,10 +725,10 @@ { double d = std::abs (U.xelem(j,j)); if (xisnan (minU) || d < minU) - minU = d; + minU = d; if (xisnan (maxU) || d > maxU) - maxU = d; + maxU = d; } double rcond = (minU / maxU); @@ -737,9 +737,9 @@ if (rcond_plus_one == 1.0 || xisnan (rcond)) { (*current_liboctave_warning_handler) - ("eigs: 'A - sigma*B' is singular, indicating sigma is exactly"); + ("eigs: 'A - sigma*B' is singular, indicating sigma is exactly"); (*current_liboctave_warning_handler) - (" an eigenvalue. Convergence is not guaranteed"); + (" an eigenvalue. Convergence is not guaranteed"); } return true; @@ -748,12 +748,12 @@ template <class M> octave_idx_type EigsRealSymmetricMatrix (const M& m, const std::string typ, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, Matrix &eig_vec, - ColumnVector &eig_val, const M& _b, - ColumnVector &permB, ColumnVector &resid, - std::ostream& os, double tol, int rvec, - bool cholB, int disp, int maxit) + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, Matrix &eig_vec, + ColumnVector &eig_val, const M& _b, + ColumnVector &permB, ColumnVector &resid, + std::ostream& os, double tol, int rvec, + bool cholB, int disp, int maxit) { M b(_b); octave_idx_type n = m.cols (); @@ -772,7 +772,7 @@ if (have_b && (m.rows() != b.rows() || m.rows() != b.cols())) { (*current_liboctave_error_handler) - ("eigs: B must be square and the same size as A"); + ("eigs: B must be square and the same size as A"); return -1; } @@ -787,7 +787,7 @@ if (n < 3) { (*current_liboctave_error_handler) - ("eigs: n must be at least 3"); + ("eigs: n must be at least 3"); return -1; } @@ -796,24 +796,24 @@ p = k * 2; if (p < 20) - p = 20; + p = 20; if (p > n - 1) - p = n - 1 ; + p = n - 1 ; } if (k < 1 || k > n - 2) { (*current_liboctave_error_handler) - ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1-1).\n" - " Use 'eig(full(A))' instead"); + ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1-1).\n" + " Use 'eig(full(A))' instead"); return -1; } if (p <= k || p >= n) { (*current_liboctave_error_handler) - ("eigs: opts.p must be greater than k and less than n"); + ("eigs: opts.p must be greater than k and less than n"); return -1; } @@ -821,27 +821,27 @@ { // Check the we really have a permutation vector if (permB.length() != n) - { - (*current_liboctave_error_handler) - ("eigs: permB vector invalid"); - return -1; - } + { + (*current_liboctave_error_handler) + ("eigs: permB vector invalid"); + return -1; + } else - { - Array<bool> checked(n,false); - for (octave_idx_type i = 0; i < n; i++) - { - octave_idx_type bidx = - static_cast<octave_idx_type> (permB(i)); - if (checked(bidx) || bidx < 0 || - bidx >= n || D_NINT (bidx) != bidx) - { - (*current_liboctave_error_handler) - ("eigs: permB vector invalid"); - return -1; - } - } - } + { + Array<bool> checked(n,false); + for (octave_idx_type i = 0; i < n; i++) + { + octave_idx_type bidx = + static_cast<octave_idx_type> (permB(i)); + if (checked(bidx) || bidx < 0 || + bidx >= n || D_NINT (bidx) != bidx) + { + (*current_liboctave_error_handler) + ("eigs: permB vector invalid"); + return -1; + } + } + } } if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" && @@ -849,14 +849,14 @@ typ != "SI") { (*current_liboctave_error_handler) - ("eigs: unrecognized sigma value"); + ("eigs: unrecognized sigma value"); return -1; } if (typ == "LI" || typ == "SI" || typ == "LR" || typ == "SR") { (*current_liboctave_error_handler) - ("eigs: invalid sigma value for real symmetric problem"); + ("eigs: invalid sigma value for real symmetric problem"); return -1; } @@ -865,25 +865,25 @@ // See Note 3 dsaupd note3 = true; if (cholB) - { - bt = b; - b = b.transpose(); - if (permB.length() == 0) - { - permB = ColumnVector(n); - for (octave_idx_type i = 0; i < n; i++) - permB(i) = i; - } - } + { + bt = b; + b = b.transpose(); + if (permB.length() == 0) + { + permB = ColumnVector(n); + for (octave_idx_type i = 0; i < n; i++) + permB(i) = i; + } + } else - { - if (! make_cholb(b, bt, permB)) - { - (*current_liboctave_error_handler) - ("eigs: The matrix B is not positive definite"); - return -1; - } - } + { + if (! make_cholb(b, bt, permB)) + { + (*current_liboctave_error_handler) + ("eigs: The matrix B is not positive definite"); + return -1; + } + } } Array<octave_idx_type> ip (11); @@ -917,66 +917,66 @@ do { F77_FUNC (dsaupd, DSAUPD) - (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, - F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), - k, tol, presid, p, v, n, iparam, - ipntr, workd, workl, lwork, info - F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); + (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, + F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), + k, tol, presid, p, v, n, iparam, + ipntr, workd, workl, lwork, info + F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("eigs: unrecoverable exception encountered in dsaupd"); - return -1; - } + { + (*current_liboctave_error_handler) + ("eigs: unrecoverable exception encountered in dsaupd"); + return -1; + } if (disp > 0 && !xisnan(workl[iptr(5)-1])) - { - if (iter++) - { - os << "Iteration " << iter - 1 << - ": a few Ritz values of the " << p << "-by-" << - p << " matrix\n"; - for (int i = 0 ; i < k; i++) - os << " " << workl[iptr(5)+i-1] << "\n"; - } - - // This is a kludge, as ARPACK doesn't give its - // iteration pointer. But as workl[iptr(5)-1] is - // an output value updated at each iteration, setting - // a value in this array to NaN and testing for it - // is a way of obtaining the iteration counter. - if (ido != 99) - workl[iptr(5)-1] = octave_NaN; - } + { + if (iter++) + { + os << "Iteration " << iter - 1 << + ": a few Ritz values of the " << p << "-by-" << + p << " matrix\n"; + for (int i = 0 ; i < k; i++) + os << " " << workl[iptr(5)+i-1] << "\n"; + } + + // This is a kludge, as ARPACK doesn't give its + // iteration pointer. But as workl[iptr(5)-1] is + // an output value updated at each iteration, setting + // a value in this array to NaN and testing for it + // is a way of obtaining the iteration counter. + if (ido != 99) + workl[iptr(5)-1] = octave_NaN; + } if (ido == -1 || ido == 1 || ido == 2) - { - if (have_b) - { - Matrix mtmp (n,1); - for (octave_idx_type i = 0; i < n; i++) - mtmp(i,0) = workd[i + iptr(0) - 1]; - - mtmp = utsolve(bt, permB, m * ltsolve(b, permB, mtmp)); - - for (octave_idx_type i = 0; i < n; i++) - workd[i+iptr(1)-1] = mtmp(i,0); - } - else if (!vector_product (m, workd + iptr(0) - 1, - workd + iptr(1) - 1)) - break; - } + { + if (have_b) + { + Matrix mtmp (n,1); + for (octave_idx_type i = 0; i < n; i++) + mtmp(i,0) = workd[i + iptr(0) - 1]; + + mtmp = utsolve(bt, permB, m * ltsolve(b, permB, mtmp)); + + for (octave_idx_type i = 0; i < n; i++) + workd[i+iptr(1)-1] = mtmp(i,0); + } + else if (!vector_product (m, workd + iptr(0) - 1, + workd + iptr(1) - 1)) + break; + } else - { - if (info < 0) - { - (*current_liboctave_error_handler) - ("eigs: error %d in dsaupd", info); - return -1; - } - break; - } + { + if (info < 0) + { + (*current_liboctave_error_handler) + ("eigs: error %d in dsaupd", info); + return -1; + } + break; + } } while (1); @@ -1008,59 +1008,59 @@ if (f77_exception_encountered) { (*current_liboctave_error_handler) - ("eigs: unrecoverable exception encountered in dseupd"); + ("eigs: unrecoverable exception encountered in dseupd"); return -1; } else { if (info2 == 0) - { - octave_idx_type k2 = k / 2; - if (typ != "SM" && typ != "BE") - { - for (octave_idx_type i = 0; i < k2; i++) - { - double dtmp = d[i]; - d[i] = d[k - i - 1]; - d[k - i - 1] = dtmp; - } - } - - if (rvec) - { - if (typ != "SM" && typ != "BE") - { - OCTAVE_LOCAL_BUFFER (double, dtmp, n); - - for (octave_idx_type i = 0; i < k2; i++) - { - octave_idx_type off1 = i * n; - octave_idx_type off2 = (k - i - 1) * n; - - if (off1 == off2) - continue; - - for (octave_idx_type j = 0; j < n; j++) - dtmp[j] = z[off1 + j]; - - for (octave_idx_type j = 0; j < n; j++) - z[off1 + j] = z[off2 + j]; - - for (octave_idx_type j = 0; j < n; j++) - z[off2 + j] = dtmp[j]; - } - } - - if (note3) - eig_vec = ltsolve(b, permB, eig_vec); - } - } + { + octave_idx_type k2 = k / 2; + if (typ != "SM" && typ != "BE") + { + for (octave_idx_type i = 0; i < k2; i++) + { + double dtmp = d[i]; + d[i] = d[k - i - 1]; + d[k - i - 1] = dtmp; + } + } + + if (rvec) + { + if (typ != "SM" && typ != "BE") + { + OCTAVE_LOCAL_BUFFER (double, dtmp, n); + + for (octave_idx_type i = 0; i < k2; i++) + { + octave_idx_type off1 = i * n; + octave_idx_type off2 = (k - i - 1) * n; + + if (off1 == off2) + continue; + + for (octave_idx_type j = 0; j < n; j++) + dtmp[j] = z[off1 + j]; + + for (octave_idx_type j = 0; j < n; j++) + z[off1 + j] = z[off2 + j]; + + for (octave_idx_type j = 0; j < n; j++) + z[off2 + j] = dtmp[j]; + } + } + + if (note3) + eig_vec = ltsolve(b, permB, eig_vec); + } + } else - { - (*current_liboctave_error_handler) - ("eigs: error %d in dseupd", info2); - return -1; - } + { + (*current_liboctave_error_handler) + ("eigs: error %d in dseupd", info2); + return -1; + } } return ip(4); @@ -1069,12 +1069,12 @@ template <class M> octave_idx_type EigsRealSymmetricMatrixShift (const M& m, double sigma, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, Matrix &eig_vec, - ColumnVector &eig_val, const M& _b, - ColumnVector &permB, ColumnVector &resid, - std::ostream& os, double tol, int rvec, - bool cholB, int disp, int maxit) + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, Matrix &eig_vec, + ColumnVector &eig_val, const M& _b, + ColumnVector &permB, ColumnVector &resid, + std::ostream& os, double tol, int rvec, + bool cholB, int disp, int maxit) { M b(_b); octave_idx_type n = m.cols (); @@ -1090,15 +1090,15 @@ if (have_b && (m.rows() != b.rows() || m.rows() != b.cols())) { (*current_liboctave_error_handler) - ("eigs: B must be square and the same size as A"); + ("eigs: B must be square and the same size as A"); return -1; } // FIXME: The "SM" type for mode 1 seems unstable though faster!! //if (! std::abs (sigma)) // return EigsRealSymmetricMatrix (m, "SM", k, p, info, eig_vec, eig_val, - // _b, permB, resid, os, tol, rvec, cholB, - // disp, maxit); + // _b, permB, resid, os, tol, rvec, cholB, + // disp, maxit); if (resid.is_empty()) { @@ -1111,15 +1111,15 @@ if (n < 3) { (*current_liboctave_error_handler) - ("eigs: n must be at least 3"); + ("eigs: n must be at least 3"); return -1; } if (k <= 0 || k >= n - 1) { (*current_liboctave_error_handler) - ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1-1).\n" - " Use 'eig(full(A))' instead"); + ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1-1).\n" + " Use 'eig(full(A))' instead"); return -1; } @@ -1128,16 +1128,16 @@ p = k * 2; if (p < 20) - p = 20; + p = 20; if (p > n - 1) - p = n - 1 ; + p = n - 1 ; } if (p <= k || p >= n) { (*current_liboctave_error_handler) - ("eigs: opts.p must be greater than k and less than n"); + ("eigs: opts.p must be greater than k and less than n"); return -1; } @@ -1145,26 +1145,26 @@ { // Check the we really have a permutation vector if (permB.length() != n) - { - (*current_liboctave_error_handler) ("eigs: permB vector invalid"); - return -1; - } + { + (*current_liboctave_error_handler) ("eigs: permB vector invalid"); + return -1; + } else - { - Array<bool> checked(n,false); - for (octave_idx_type i = 0; i < n; i++) - { - octave_idx_type bidx = - static_cast<octave_idx_type> (permB(i)); - if (checked(bidx) || bidx < 0 || - bidx >= n || D_NINT (bidx) != bidx) - { - (*current_liboctave_error_handler) - ("eigs: permB vector invalid"); - return -1; - } - } - } + { + Array<bool> checked(n,false); + for (octave_idx_type i = 0; i < n; i++) + { + octave_idx_type bidx = + static_cast<octave_idx_type> (permB(i)); + if (checked(bidx) || bidx < 0 || + bidx >= n || D_NINT (bidx) != bidx) + { + (*current_liboctave_error_handler) + ("eigs: permB vector invalid"); + return -1; + } + } + } } char bmat = 'I'; @@ -1210,110 +1210,110 @@ do { F77_FUNC (dsaupd, DSAUPD) - (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, - F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), - k, tol, presid, p, v, n, iparam, - ipntr, workd, workl, lwork, info - F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); + (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, + F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), + k, tol, presid, p, v, n, iparam, + ipntr, workd, workl, lwork, info + F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("eigs: unrecoverable exception encountered in dsaupd"); - return -1; - } + { + (*current_liboctave_error_handler) + ("eigs: unrecoverable exception encountered in dsaupd"); + return -1; + } if (disp > 0 && !xisnan(workl[iptr(5)-1])) - { - if (iter++) - { - os << "Iteration " << iter - 1 << - ": a few Ritz values of the " << p << "-by-" << - p << " matrix\n"; - for (int i = 0 ; i < k; i++) - os << " " << workl[iptr(5)+i-1] << "\n"; - } - - // This is a kludge, as ARPACK doesn't give its - // iteration pointer. But as workl[iptr(5)-1] is - // an output value updated at each iteration, setting - // a value in this array to NaN and testing for it - // is a way of obtaining the iteration counter. - if (ido != 99) - workl[iptr(5)-1] = octave_NaN; - } + { + if (iter++) + { + os << "Iteration " << iter - 1 << + ": a few Ritz values of the " << p << "-by-" << + p << " matrix\n"; + for (int i = 0 ; i < k; i++) + os << " " << workl[iptr(5)+i-1] << "\n"; + } + + // This is a kludge, as ARPACK doesn't give its + // iteration pointer. But as workl[iptr(5)-1] is + // an output value updated at each iteration, setting + // a value in this array to NaN and testing for it + // is a way of obtaining the iteration counter. + if (ido != 99) + workl[iptr(5)-1] = octave_NaN; + } if (ido == -1 || ido == 1 || ido == 2) - { - if (have_b) - { - if (ido == -1) - { - OCTAVE_LOCAL_BUFFER (double, dtmp, n); - - vector_product (m, workd+iptr(0)-1, dtmp); - - Matrix tmp(n, 1); - - for (octave_idx_type i = 0; i < n; i++) - tmp(i,0) = dtmp[P[i]]; - - lusolve (L, U, tmp); - - double *ip2 = workd+iptr(1)-1; - for (octave_idx_type i = 0; i < n; i++) - ip2[Q[i]] = tmp(i,0); - } - else if (ido == 2) - vector_product (b, workd+iptr(0)-1, workd+iptr(1)-1); - else - { - double *ip2 = workd+iptr(2)-1; - Matrix tmp(n, 1); - - for (octave_idx_type i = 0; i < n; i++) - tmp(i,0) = ip2[P[i]]; - - lusolve (L, U, tmp); - - ip2 = workd+iptr(1)-1; - for (octave_idx_type i = 0; i < n; i++) - ip2[Q[i]] = tmp(i,0); - } - } - else - { - if (ido == 2) - { - for (octave_idx_type i = 0; i < n; i++) - workd[iptr(0) + i - 1] = workd[iptr(1) + i - 1]; - } - else - { - double *ip2 = workd+iptr(0)-1; - Matrix tmp(n, 1); - - for (octave_idx_type i = 0; i < n; i++) - tmp(i,0) = ip2[P[i]]; - - lusolve (L, U, tmp); - - ip2 = workd+iptr(1)-1; - for (octave_idx_type i = 0; i < n; i++) - ip2[Q[i]] = tmp(i,0); - } - } - } + { + if (have_b) + { + if (ido == -1) + { + OCTAVE_LOCAL_BUFFER (double, dtmp, n); + + vector_product (m, workd+iptr(0)-1, dtmp); + + Matrix tmp(n, 1); + + for (octave_idx_type i = 0; i < n; i++) + tmp(i,0) = dtmp[P[i]]; + + lusolve (L, U, tmp); + + double *ip2 = workd+iptr(1)-1; + for (octave_idx_type i = 0; i < n; i++) + ip2[Q[i]] = tmp(i,0); + } + else if (ido == 2) + vector_product (b, workd+iptr(0)-1, workd+iptr(1)-1); + else + { + double *ip2 = workd+iptr(2)-1; + Matrix tmp(n, 1); + + for (octave_idx_type i = 0; i < n; i++) + tmp(i,0) = ip2[P[i]]; + + lusolve (L, U, tmp); + + ip2 = workd+iptr(1)-1; + for (octave_idx_type i = 0; i < n; i++) + ip2[Q[i]] = tmp(i,0); + } + } + else + { + if (ido == 2) + { + for (octave_idx_type i = 0; i < n; i++) + workd[iptr(0) + i - 1] = workd[iptr(1) + i - 1]; + } + else + { + double *ip2 = workd+iptr(0)-1; + Matrix tmp(n, 1); + + for (octave_idx_type i = 0; i < n; i++) + tmp(i,0) = ip2[P[i]]; + + lusolve (L, U, tmp); + + ip2 = workd+iptr(1)-1; + for (octave_idx_type i = 0; i < n; i++) + ip2[Q[i]] = tmp(i,0); + } + } + } else - { - if (info < 0) - { - (*current_liboctave_error_handler) - ("eigs: error %d in dsaupd", info); - return -1; - } - break; - } + { + if (info < 0) + { + (*current_liboctave_error_handler) + ("eigs: error %d in dsaupd", info); + return -1; + } + break; + } } while (1); @@ -1328,7 +1328,7 @@ // avoid problems. Array<octave_idx_type> s (p); octave_idx_type *sel = s.fortran_vec (); - + eig_vec.resize (n, k); double *z = eig_vec.fortran_vec (); @@ -1345,50 +1345,50 @@ if (f77_exception_encountered) { (*current_liboctave_error_handler) - ("eigs: unrecoverable exception encountered in dseupd"); + ("eigs: unrecoverable exception encountered in dseupd"); return -1; } else { if (info2 == 0) - { - octave_idx_type k2 = k / 2; - for (octave_idx_type i = 0; i < k2; i++) - { - double dtmp = d[i]; - d[i] = d[k - i - 1]; - d[k - i - 1] = dtmp; - } - - if (rvec) - { - OCTAVE_LOCAL_BUFFER (double, dtmp, n); - - for (octave_idx_type i = 0; i < k2; i++) - { - octave_idx_type off1 = i * n; - octave_idx_type off2 = (k - i - 1) * n; - - if (off1 == off2) - continue; - - for (octave_idx_type j = 0; j < n; j++) - dtmp[j] = z[off1 + j]; - - for (octave_idx_type j = 0; j < n; j++) - z[off1 + j] = z[off2 + j]; - - for (octave_idx_type j = 0; j < n; j++) - z[off2 + j] = dtmp[j]; - } - } - } + { + octave_idx_type k2 = k / 2; + for (octave_idx_type i = 0; i < k2; i++) + { + double dtmp = d[i]; + d[i] = d[k - i - 1]; + d[k - i - 1] = dtmp; + } + + if (rvec) + { + OCTAVE_LOCAL_BUFFER (double, dtmp, n); + + for (octave_idx_type i = 0; i < k2; i++) + { + octave_idx_type off1 = i * n; + octave_idx_type off2 = (k - i - 1) * n; + + if (off1 == off2) + continue; + + for (octave_idx_type j = 0; j < n; j++) + dtmp[j] = z[off1 + j]; + + for (octave_idx_type j = 0; j < n; j++) + z[off1 + j] = z[off2 + j]; + + for (octave_idx_type j = 0; j < n; j++) + z[off2 + j] = dtmp[j]; + } + } + } else - { - (*current_liboctave_error_handler) - ("eigs: error %d in dseupd", info2); - return -1; - } + { + (*current_liboctave_error_handler) + ("eigs: error %d in dseupd", info2); + return -1; + } } return ip(4); @@ -1396,12 +1396,12 @@ octave_idx_type EigsRealSymmetricFunc (EigsFunc fun, octave_idx_type n, - const std::string &_typ, double sigma, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, Matrix &eig_vec, - ColumnVector &eig_val, ColumnVector &resid, - std::ostream& os, double tol, int rvec, - bool /* cholB */, int disp, int maxit) + const std::string &_typ, double sigma, + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, Matrix &eig_vec, + ColumnVector &eig_val, ColumnVector &resid, + std::ostream& os, double tol, int rvec, + bool /* cholB */, int disp, int maxit) { std::string typ (_typ); bool have_sigma = (sigma ? true : false); @@ -1420,7 +1420,7 @@ if (n < 3) { (*current_liboctave_error_handler) - ("eigs: n must be at least 3"); + ("eigs: n must be at least 3"); return -1; } @@ -1429,48 +1429,48 @@ p = k * 2; if (p < 20) - p = 20; + p = 20; if (p > n - 1) - p = n - 1 ; + p = n - 1 ; } if (k <= 0 || k >= n - 1) { (*current_liboctave_error_handler) - ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1).\n" - " Use 'eig(full(A))' instead"); + ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1).\n" + " Use 'eig(full(A))' instead"); return -1; } if (p <= k || p >= n) { (*current_liboctave_error_handler) - ("eigs: opts.p must be greater than k and less than n"); + ("eigs: opts.p must be greater than k and less than n"); return -1; } if (! have_sigma) { if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" && - typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" && - typ != "SI") - (*current_liboctave_error_handler) - ("eigs: unrecognized sigma value"); + typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" && + typ != "SI") + (*current_liboctave_error_handler) + ("eigs: unrecognized sigma value"); if (typ == "LI" || typ == "SI" || typ == "LR" || typ == "SR") - { - (*current_liboctave_error_handler) - ("eigs: invalid sigma value for real symmetric problem"); - return -1; - } + { + (*current_liboctave_error_handler) + ("eigs: invalid sigma value for real symmetric problem"); + return -1; + } if (typ == "SM") - { - typ = "LM"; - sigma = 0.; - mode = 3; - } + { + typ = "LM"; + sigma = 0.; + mode = 3; + } } else if (! std::abs (sigma)) typ = "SM"; @@ -1511,67 +1511,67 @@ do { F77_FUNC (dsaupd, DSAUPD) - (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, - F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), - k, tol, presid, p, v, n, iparam, - ipntr, workd, workl, lwork, info - F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); + (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, + F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), + k, tol, presid, p, v, n, iparam, + ipntr, workd, workl, lwork, info + F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("eigs: unrecoverable exception encountered in dsaupd"); - return -1; - } + { + (*current_liboctave_error_handler) + ("eigs: unrecoverable exception encountered in dsaupd"); + return -1; + } if (disp > 0 && !xisnan(workl[iptr(5)-1])) - { - if (iter++) - { - os << "Iteration " << iter - 1 << - ": a few Ritz values of the " << p << "-by-" << - p << " matrix\n"; - for (int i = 0 ; i < k; i++) - os << " " << workl[iptr(5)+i-1] << "\n"; - } - - // This is a kludge, as ARPACK doesn't give its - // iteration pointer. But as workl[iptr(5)-1] is - // an output value updated at each iteration, setting - // a value in this array to NaN and testing for it - // is a way of obtaining the iteration counter. - if (ido != 99) - workl[iptr(5)-1] = octave_NaN; - } + { + if (iter++) + { + os << "Iteration " << iter - 1 << + ": a few Ritz values of the " << p << "-by-" << + p << " matrix\n"; + for (int i = 0 ; i < k; i++) + os << " " << workl[iptr(5)+i-1] << "\n"; + } + + // This is a kludge, as ARPACK doesn't give its + // iteration pointer. But as workl[iptr(5)-1] is + // an output value updated at each iteration, setting + // a value in this array to NaN and testing for it + // is a way of obtaining the iteration counter. + if (ido != 99) + workl[iptr(5)-1] = octave_NaN; + } if (ido == -1 || ido == 1 || ido == 2) - { - double *ip2 = workd + iptr(0) - 1; - ColumnVector x(n); - - for (octave_idx_type i = 0; i < n; i++) - x(i) = *ip2++; - - ColumnVector y = fun (x, err); - - if (err) - return false; - - ip2 = workd + iptr(1) - 1; - for (octave_idx_type i = 0; i < n; i++) - *ip2++ = y(i); - } + { + double *ip2 = workd + iptr(0) - 1; + ColumnVector x(n); + + for (octave_idx_type i = 0; i < n; i++) + x(i) = *ip2++; + + ColumnVector y = fun (x, err); + + if (err) + return false; + + ip2 = workd + iptr(1) - 1; + for (octave_idx_type i = 0; i < n; i++) + *ip2++ = y(i); + } else - { - if (info < 0) - { - (*current_liboctave_error_handler) - ("eigs: error %d in dsaupd", info); - return -1; - } - break; - } + { + if (info < 0) + { + (*current_liboctave_error_handler) + ("eigs: error %d in dsaupd", info); + return -1; + } + break; + } } while (1); @@ -1586,7 +1586,7 @@ // avoid problems. Array<octave_idx_type> s (p); octave_idx_type *sel = s.fortran_vec (); - + eig_vec.resize (n, k); double *z = eig_vec.fortran_vec (); @@ -1603,56 +1603,56 @@ if (f77_exception_encountered) { (*current_liboctave_error_handler) - ("eigs: unrecoverable exception encountered in dseupd"); + ("eigs: unrecoverable exception encountered in dseupd"); return -1; } else { if (info2 == 0) - { - octave_idx_type k2 = k / 2; - if (typ != "SM" && typ != "BE") - { - for (octave_idx_type i = 0; i < k2; i++) - { - double dtmp = d[i]; - d[i] = d[k - i - 1]; - d[k - i - 1] = dtmp; - } - } - - if (rvec) - { - if (typ != "SM" && typ != "BE") - { - OCTAVE_LOCAL_BUFFER (double, dtmp, n); - - for (octave_idx_type i = 0; i < k2; i++) - { - octave_idx_type off1 = i * n; - octave_idx_type off2 = (k - i - 1) * n; - - if (off1 == off2) - continue; - - for (octave_idx_type j = 0; j < n; j++) - dtmp[j] = z[off1 + j]; - - for (octave_idx_type j = 0; j < n; j++) - z[off1 + j] = z[off2 + j]; - - for (octave_idx_type j = 0; j < n; j++) - z[off2 + j] = dtmp[j]; - } - } - } - } + { + octave_idx_type k2 = k / 2; + if (typ != "SM" && typ != "BE") + { + for (octave_idx_type i = 0; i < k2; i++) + { + double dtmp = d[i]; + d[i] = d[k - i - 1]; + d[k - i - 1] = dtmp; + } + } + + if (rvec) + { + if (typ != "SM" && typ != "BE") + { + OCTAVE_LOCAL_BUFFER (double, dtmp, n); + + for (octave_idx_type i = 0; i < k2; i++) + { + octave_idx_type off1 = i * n; + octave_idx_type off2 = (k - i - 1) * n; + + if (off1 == off2) + continue; + + for (octave_idx_type j = 0; j < n; j++) + dtmp[j] = z[off1 + j]; + + for (octave_idx_type j = 0; j < n; j++) + z[off1 + j] = z[off2 + j]; + + for (octave_idx_type j = 0; j < n; j++) + z[off2 + j] = dtmp[j]; + } + } + } + } else - { - (*current_liboctave_error_handler) - ("eigs: error %d in dseupd", info2); - return -1; - } + { + (*current_liboctave_error_handler) + ("eigs: error %d in dseupd", info2); + return -1; + } } return ip(4); @@ -1661,12 +1661,12 @@ template <class M> octave_idx_type EigsRealNonSymmetricMatrix (const M& m, const std::string typ, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, ComplexMatrix &eig_vec, - ComplexColumnVector &eig_val, const M& _b, - ColumnVector &permB, ColumnVector &resid, - std::ostream& os, double tol, int rvec, - bool cholB, int disp, int maxit) + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, ComplexMatrix &eig_vec, + ComplexColumnVector &eig_val, const M& _b, + ColumnVector &permB, ColumnVector &resid, + std::ostream& os, double tol, int rvec, + bool cholB, int disp, int maxit) { M b(_b); octave_idx_type n = m.cols (); @@ -1686,7 +1686,7 @@ if (have_b && (m.rows() != b.rows() || m.rows() != b.cols())) { (*current_liboctave_error_handler) - ("eigs: B must be square and the same size as A"); + ("eigs: B must be square and the same size as A"); return -1; } @@ -1701,7 +1701,7 @@ if (n < 3) { (*current_liboctave_error_handler) - ("eigs: n must be at least 3"); + ("eigs: n must be at least 3"); return -1; } @@ -1710,24 +1710,24 @@ p = k * 2 + 1; if (p < 20) - p = 20; + p = 20; if (p > n - 1) - p = n - 1 ; + p = n - 1 ; } if (k <= 0 || k >= n - 1) { (*current_liboctave_error_handler) - ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1).\n" - " Use 'eig(full(A))' instead"); + ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1).\n" + " Use 'eig(full(A))' instead"); return -1; } if (p <= k || p >= n) { (*current_liboctave_error_handler) - ("eigs: opts.p must be greater than k and less than n"); + ("eigs: opts.p must be greater than k and less than n"); return -1; } @@ -1735,27 +1735,27 @@ { // Check the we really have a permutation vector if (permB.length() != n) - { - (*current_liboctave_error_handler) - ("eigs: permB vector invalid"); - return -1; - } + { + (*current_liboctave_error_handler) + ("eigs: permB vector invalid"); + return -1; + } else - { - Array<bool> checked(n,false); - for (octave_idx_type i = 0; i < n; i++) - { - octave_idx_type bidx = - static_cast<octave_idx_type> (permB(i)); - if (checked(bidx) || bidx < 0 || - bidx >= n || D_NINT (bidx) != bidx) - { - (*current_liboctave_error_handler) - ("eigs: permB vector invalid"); - return -1; - } - } - } + { + Array<bool> checked(n,false); + for (octave_idx_type i = 0; i < n; i++) + { + octave_idx_type bidx = + static_cast<octave_idx_type> (permB(i)); + if (checked(bidx) || bidx < 0 || + bidx >= n || D_NINT (bidx) != bidx) + { + (*current_liboctave_error_handler) + ("eigs: permB vector invalid"); + return -1; + } + } + } } if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" && @@ -1763,14 +1763,14 @@ typ != "SI") { (*current_liboctave_error_handler) - ("eigs: unrecognized sigma value"); + ("eigs: unrecognized sigma value"); return -1; } if (typ == "LA" || typ == "SA" || typ == "BE") { (*current_liboctave_error_handler) - ("eigs: invalid sigma value for unsymmetric problem"); + ("eigs: invalid sigma value for unsymmetric problem"); return -1; } @@ -1779,25 +1779,25 @@ // See Note 3 dsaupd note3 = true; if (cholB) - { - bt = b; - b = b.transpose(); - if (permB.length() == 0) - { - permB = ColumnVector(n); - for (octave_idx_type i = 0; i < n; i++) - permB(i) = i; - } - } + { + bt = b; + b = b.transpose(); + if (permB.length() == 0) + { + permB = ColumnVector(n); + for (octave_idx_type i = 0; i < n; i++) + permB(i) = i; + } + } else - { - if (! make_cholb(b, bt, permB)) - { - (*current_liboctave_error_handler) - ("eigs: The matrix B is not positive definite"); - return -1; - } - } + { + if (! make_cholb(b, bt, permB)) + { + (*current_liboctave_error_handler) + ("eigs: The matrix B is not positive definite"); + return -1; + } + } } Array<octave_idx_type> ip (11); @@ -1831,66 +1831,66 @@ do { F77_FUNC (dnaupd, DNAUPD) - (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, - F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), - k, tol, presid, p, v, n, iparam, - ipntr, workd, workl, lwork, info - F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); + (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, + F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), + k, tol, presid, p, v, n, iparam, + ipntr, workd, workl, lwork, info + F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("eigs: unrecoverable exception encountered in dnaupd"); - return -1; - } + { + (*current_liboctave_error_handler) + ("eigs: unrecoverable exception encountered in dnaupd"); + return -1; + } if (disp > 0 && !xisnan(workl[iptr(5)-1])) - { - if (iter++) - { - os << "Iteration " << iter - 1 << - ": a few Ritz values of the " << p << "-by-" << - p << " matrix\n"; - for (int i = 0 ; i < k; i++) - os << " " << workl[iptr(5)+i-1] << "\n"; - } - - // This is a kludge, as ARPACK doesn't give its - // iteration pointer. But as workl[iptr(5)-1] is - // an output value updated at each iteration, setting - // a value in this array to NaN and testing for it - // is a way of obtaining the iteration counter. - if (ido != 99) - workl[iptr(5)-1] = octave_NaN; - } + { + if (iter++) + { + os << "Iteration " << iter - 1 << + ": a few Ritz values of the " << p << "-by-" << + p << " matrix\n"; + for (int i = 0 ; i < k; i++) + os << " " << workl[iptr(5)+i-1] << "\n"; + } + + // This is a kludge, as ARPACK doesn't give its + // iteration pointer. But as workl[iptr(5)-1] is + // an output value updated at each iteration, setting + // a value in this array to NaN and testing for it + // is a way of obtaining the iteration counter. + if (ido != 99) + workl[iptr(5)-1] = octave_NaN; + } if (ido == -1 || ido == 1 || ido == 2) - { - if (have_b) - { - Matrix mtmp (n,1); - for (octave_idx_type i = 0; i < n; i++) - mtmp(i,0) = workd[i + iptr(0) - 1]; - - mtmp = utsolve(bt, permB, m * ltsolve(b, permB, mtmp)); - - for (octave_idx_type i = 0; i < n; i++) - workd[i+iptr(1)-1] = mtmp(i,0); - } - else if (!vector_product (m, workd + iptr(0) - 1, - workd + iptr(1) - 1)) - break; - } + { + if (have_b) + { + Matrix mtmp (n,1); + for (octave_idx_type i = 0; i < n; i++) + mtmp(i,0) = workd[i + iptr(0) - 1]; + + mtmp = utsolve(bt, permB, m * ltsolve(b, permB, mtmp)); + + for (octave_idx_type i = 0; i < n; i++) + workd[i+iptr(1)-1] = mtmp(i,0); + } + else if (!vector_product (m, workd + iptr(0) - 1, + workd + iptr(1) - 1)) + break; + } else - { - if (info < 0) - { - (*current_liboctave_error_handler) - ("eigs: error %d in dnaupd", info); - return -1; - } - break; - } + { + if (info < 0) + { + (*current_liboctave_error_handler) + ("eigs: error %d in dnaupd", info); + return -1; + } + break; + } } while (1); @@ -1925,7 +1925,7 @@ if (f77_exception_encountered) { (*current_liboctave_error_handler) - ("eigs: unrecoverable exception encountered in dneupd"); + ("eigs: unrecoverable exception encountered in dneupd"); return -1; } else @@ -1934,87 +1934,87 @@ Complex *d = eig_val.fortran_vec (); if (info2 == 0) - { - octave_idx_type jj = 0; - for (octave_idx_type i = 0; i < k+1; i++) - { - if (dr[i] == 0.0 && di[i] == 0.0 && jj == 0) - jj++; - else - d [i-jj] = Complex (dr[i], di[i]); - } - if (jj == 0 && !rvec) - for (octave_idx_type i = 0; i < k; i++) - d[i] = d[i+1]; - - octave_idx_type k2 = k / 2; - for (octave_idx_type i = 0; i < k2; i++) - { - Complex dtmp = d[i]; - d[i] = d[k - i - 1]; - d[k - i - 1] = dtmp; - } - eig_val.resize(k); - - if (rvec) - { - OCTAVE_LOCAL_BUFFER (double, dtmp, n); - - for (octave_idx_type i = 0; i < k2; i++) - { - octave_idx_type off1 = i * n; - octave_idx_type off2 = (k - i - 1) * n; - - if (off1 == off2) - continue; - - for (octave_idx_type j = 0; j < n; j++) - dtmp[j] = z[off1 + j]; - - for (octave_idx_type j = 0; j < n; j++) - z[off1 + j] = z[off2 + j]; - - for (octave_idx_type j = 0; j < n; j++) - z[off2 + j] = dtmp[j]; - } - - eig_vec.resize (n, k); - octave_idx_type i = 0; - while (i < k) - { - octave_idx_type off1 = i * n; - octave_idx_type off2 = (i+1) * n; - if (std::imag(eig_val(i)) == 0) - { - for (octave_idx_type j = 0; j < n; j++) - eig_vec(j,i) = - Complex(z[j+off1],0.); - i++; - } - else - { - for (octave_idx_type j = 0; j < n; j++) - { - eig_vec(j,i) = - Complex(z[j+off1],z[j+off2]); - if (i < k - 1) - eig_vec(j,i+1) = - Complex(z[j+off1],-z[j+off2]); - } - i+=2; - } - } - - if (note3) - eig_vec = ltsolve(M (b), permB, eig_vec); - } - } + { + octave_idx_type jj = 0; + for (octave_idx_type i = 0; i < k+1; i++) + { + if (dr[i] == 0.0 && di[i] == 0.0 && jj == 0) + jj++; + else + d [i-jj] = Complex (dr[i], di[i]); + } + if (jj == 0 && !rvec) + for (octave_idx_type i = 0; i < k; i++) + d[i] = d[i+1]; + + octave_idx_type k2 = k / 2; + for (octave_idx_type i = 0; i < k2; i++) + { + Complex dtmp = d[i]; + d[i] = d[k - i - 1]; + d[k - i - 1] = dtmp; + } + eig_val.resize(k); + + if (rvec) + { + OCTAVE_LOCAL_BUFFER (double, dtmp, n); + + for (octave_idx_type i = 0; i < k2; i++) + { + octave_idx_type off1 = i * n; + octave_idx_type off2 = (k - i - 1) * n; + + if (off1 == off2) + continue; + + for (octave_idx_type j = 0; j < n; j++) + dtmp[j] = z[off1 + j]; + + for (octave_idx_type j = 0; j < n; j++) + z[off1 + j] = z[off2 + j]; + + for (octave_idx_type j = 0; j < n; j++) + z[off2 + j] = dtmp[j]; + } + + eig_vec.resize (n, k); + octave_idx_type i = 0; + while (i < k) + { + octave_idx_type off1 = i * n; + octave_idx_type off2 = (i+1) * n; + if (std::imag(eig_val(i)) == 0) + { + for (octave_idx_type j = 0; j < n; j++) + eig_vec(j,i) = + Complex(z[j+off1],0.); + i++; + } + else + { + for (octave_idx_type j = 0; j < n; j++) + { + eig_vec(j,i) = + Complex(z[j+off1],z[j+off2]); + if (i < k - 1) + eig_vec(j,i+1) = + Complex(z[j+off1],-z[j+off2]); + } + i+=2; + } + } + + if (note3) + eig_vec = ltsolve(M (b), permB, eig_vec); + } + } else - { - (*current_liboctave_error_handler) - ("eigs: error %d in dneupd", info2); - return -1; - } + { + (*current_liboctave_error_handler) + ("eigs: error %d in dneupd", info2); + return -1; + } } return ip(4); @@ -2023,13 +2023,13 @@ template <class M> octave_idx_type EigsRealNonSymmetricMatrixShift (const M& m, double sigmar, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, - ComplexMatrix &eig_vec, - ComplexColumnVector &eig_val, const M& _b, - ColumnVector &permB, ColumnVector &resid, - std::ostream& os, double tol, int rvec, - bool cholB, int disp, int maxit) + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, + ComplexMatrix &eig_vec, + ComplexColumnVector &eig_val, const M& _b, + ColumnVector &permB, ColumnVector &resid, + std::ostream& os, double tol, int rvec, + bool cholB, int disp, int maxit) { M b(_b); octave_idx_type n = m.cols (); @@ -2046,15 +2046,15 @@ if (have_b && (m.rows() != b.rows() || m.rows() != b.cols())) { (*current_liboctave_error_handler) - ("eigs: B must be square and the same size as A"); + ("eigs: B must be square and the same size as A"); return -1; } // FIXME: The "SM" type for mode 1 seems unstable though faster!! //if (! std::abs (sigmar)) // return EigsRealNonSymmetricMatrix (m, "SM", k, p, info, eig_vec, eig_val, - // _b, permB, resid, os, tol, rvec, cholB, - // disp, maxit); + // _b, permB, resid, os, tol, rvec, cholB, + // disp, maxit); if (resid.is_empty()) { @@ -2067,7 +2067,7 @@ if (n < 3) { (*current_liboctave_error_handler) - ("eigs: n must be at least 3"); + ("eigs: n must be at least 3"); return -1; } @@ -2076,24 +2076,24 @@ p = k * 2 + 1; if (p < 20) - p = 20; + p = 20; if (p > n - 1) - p = n - 1 ; + p = n - 1 ; } if (k <= 0 || k >= n - 1) { (*current_liboctave_error_handler) - ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1).\n" - " Use 'eig(full(A))' instead"); + ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1).\n" + " Use 'eig(full(A))' instead"); return -1; } if (p <= k || p >= n) { (*current_liboctave_error_handler) - ("eigs: opts.p must be greater than k and less than n"); + ("eigs: opts.p must be greater than k and less than n"); return -1; } @@ -2101,26 +2101,26 @@ { // Check that we really have a permutation vector if (permB.length() != n) - { - (*current_liboctave_error_handler) ("eigs: permB vector invalid"); - return -1; - } + { + (*current_liboctave_error_handler) ("eigs: permB vector invalid"); + return -1; + } else - { - Array<bool> checked(n,false); - for (octave_idx_type i = 0; i < n; i++) - { - octave_idx_type bidx = - static_cast<octave_idx_type> (permB(i)); - if (checked(bidx) || bidx < 0 || - bidx >= n || D_NINT (bidx) != bidx) - { - (*current_liboctave_error_handler) - ("eigs: permB vector invalid"); - return -1; - } - } - } + { + Array<bool> checked(n,false); + for (octave_idx_type i = 0; i < n; i++) + { + octave_idx_type bidx = + static_cast<octave_idx_type> (permB(i)); + if (checked(bidx) || bidx < 0 || + bidx >= n || D_NINT (bidx) != bidx) + { + (*current_liboctave_error_handler) + ("eigs: permB vector invalid"); + return -1; + } + } + } } char bmat = 'I'; @@ -2166,110 +2166,110 @@ do { F77_FUNC (dnaupd, DNAUPD) - (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, - F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), - k, tol, presid, p, v, n, iparam, - ipntr, workd, workl, lwork, info - F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); + (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, + F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), + k, tol, presid, p, v, n, iparam, + ipntr, workd, workl, lwork, info + F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("eigs: unrecoverable exception encountered in dsaupd"); - return -1; - } + { + (*current_liboctave_error_handler) + ("eigs: unrecoverable exception encountered in dsaupd"); + return -1; + } if (disp > 0 && !xisnan(workl[iptr(5)-1])) - { - if (iter++) - { - os << "Iteration " << iter - 1 << - ": a few Ritz values of the " << p << "-by-" << - p << " matrix\n"; - for (int i = 0 ; i < k; i++) - os << " " << workl[iptr(5)+i-1] << "\n"; - } - - // This is a kludge, as ARPACK doesn't give its - // iteration pointer. But as workl[iptr(5)-1] is - // an output value updated at each iteration, setting - // a value in this array to NaN and testing for it - // is a way of obtaining the iteration counter. - if (ido != 99) - workl[iptr(5)-1] = octave_NaN; - } + { + if (iter++) + { + os << "Iteration " << iter - 1 << + ": a few Ritz values of the " << p << "-by-" << + p << " matrix\n"; + for (int i = 0 ; i < k; i++) + os << " " << workl[iptr(5)+i-1] << "\n"; + } + + // This is a kludge, as ARPACK doesn't give its + // iteration pointer. But as workl[iptr(5)-1] is + // an output value updated at each iteration, setting + // a value in this array to NaN and testing for it + // is a way of obtaining the iteration counter. + if (ido != 99) + workl[iptr(5)-1] = octave_NaN; + } if (ido == -1 || ido == 1 || ido == 2) - { - if (have_b) - { - if (ido == -1) - { - OCTAVE_LOCAL_BUFFER (double, dtmp, n); - - vector_product (m, workd+iptr(0)-1, dtmp); - - Matrix tmp(n, 1); - - for (octave_idx_type i = 0; i < n; i++) - tmp(i,0) = dtmp[P[i]]; - - lusolve (L, U, tmp); - - double *ip2 = workd+iptr(1)-1; - for (octave_idx_type i = 0; i < n; i++) - ip2[Q[i]] = tmp(i,0); - } - else if (ido == 2) - vector_product (b, workd+iptr(0)-1, workd+iptr(1)-1); - else - { - double *ip2 = workd+iptr(2)-1; - Matrix tmp(n, 1); - - for (octave_idx_type i = 0; i < n; i++) - tmp(i,0) = ip2[P[i]]; - - lusolve (L, U, tmp); - - ip2 = workd+iptr(1)-1; - for (octave_idx_type i = 0; i < n; i++) - ip2[Q[i]] = tmp(i,0); - } - } - else - { - if (ido == 2) - { - for (octave_idx_type i = 0; i < n; i++) - workd[iptr(0) + i - 1] = workd[iptr(1) + i - 1]; - } - else - { - double *ip2 = workd+iptr(0)-1; - Matrix tmp(n, 1); - - for (octave_idx_type i = 0; i < n; i++) - tmp(i,0) = ip2[P[i]]; - - lusolve (L, U, tmp); - - ip2 = workd+iptr(1)-1; - for (octave_idx_type i = 0; i < n; i++) - ip2[Q[i]] = tmp(i,0); - } - } - } + { + if (have_b) + { + if (ido == -1) + { + OCTAVE_LOCAL_BUFFER (double, dtmp, n); + + vector_product (m, workd+iptr(0)-1, dtmp); + + Matrix tmp(n, 1); + + for (octave_idx_type i = 0; i < n; i++) + tmp(i,0) = dtmp[P[i]]; + + lusolve (L, U, tmp); + + double *ip2 = workd+iptr(1)-1; + for (octave_idx_type i = 0; i < n; i++) + ip2[Q[i]] = tmp(i,0); + } + else if (ido == 2) + vector_product (b, workd+iptr(0)-1, workd+iptr(1)-1); + else + { + double *ip2 = workd+iptr(2)-1; + Matrix tmp(n, 1); + + for (octave_idx_type i = 0; i < n; i++) + tmp(i,0) = ip2[P[i]]; + + lusolve (L, U, tmp); + + ip2 = workd+iptr(1)-1; + for (octave_idx_type i = 0; i < n; i++) + ip2[Q[i]] = tmp(i,0); + } + } + else + { + if (ido == 2) + { + for (octave_idx_type i = 0; i < n; i++) + workd[iptr(0) + i - 1] = workd[iptr(1) + i - 1]; + } + else + { + double *ip2 = workd+iptr(0)-1; + Matrix tmp(n, 1); + + for (octave_idx_type i = 0; i < n; i++) + tmp(i,0) = ip2[P[i]]; + + lusolve (L, U, tmp); + + ip2 = workd+iptr(1)-1; + for (octave_idx_type i = 0; i < n; i++) + ip2[Q[i]] = tmp(i,0); + } + } + } else - { - if (info < 0) - { - (*current_liboctave_error_handler) - ("eigs: error %d in dsaupd", info); - return -1; - } - break; - } + { + if (info < 0) + { + (*current_liboctave_error_handler) + ("eigs: error %d in dsaupd", info); + return -1; + } + break; + } } while (1); @@ -2284,7 +2284,7 @@ // avoid problems. Array<octave_idx_type> s (p); octave_idx_type *sel = s.fortran_vec (); - + Matrix eig_vec2 (n, k + 1); double *z = eig_vec2.fortran_vec (); @@ -2304,7 +2304,7 @@ if (f77_exception_encountered) { (*current_liboctave_error_handler) - ("eigs: unrecoverable exception encountered in dneupd"); + ("eigs: unrecoverable exception encountered in dneupd"); return -1; } else @@ -2313,84 +2313,84 @@ Complex *d = eig_val.fortran_vec (); if (info2 == 0) - { - octave_idx_type jj = 0; - for (octave_idx_type i = 0; i < k+1; i++) - { - if (dr[i] == 0.0 && di[i] == 0.0 && jj == 0) - jj++; - else - d [i-jj] = Complex (dr[i], di[i]); - } - if (jj == 0 && !rvec) - for (octave_idx_type i = 0; i < k; i++) - d[i] = d[i+1]; - - octave_idx_type k2 = k / 2; - for (octave_idx_type i = 0; i < k2; i++) - { - Complex dtmp = d[i]; - d[i] = d[k - i - 1]; - d[k - i - 1] = dtmp; - } - eig_val.resize(k); - - if (rvec) - { - OCTAVE_LOCAL_BUFFER (double, dtmp, n); - - for (octave_idx_type i = 0; i < k2; i++) - { - octave_idx_type off1 = i * n; - octave_idx_type off2 = (k - i - 1) * n; - - if (off1 == off2) - continue; - - for (octave_idx_type j = 0; j < n; j++) - dtmp[j] = z[off1 + j]; - - for (octave_idx_type j = 0; j < n; j++) - z[off1 + j] = z[off2 + j]; - - for (octave_idx_type j = 0; j < n; j++) - z[off2 + j] = dtmp[j]; - } - - eig_vec.resize (n, k); - octave_idx_type i = 0; - while (i < k) - { - octave_idx_type off1 = i * n; - octave_idx_type off2 = (i+1) * n; - if (std::imag(eig_val(i)) == 0) - { - for (octave_idx_type j = 0; j < n; j++) - eig_vec(j,i) = - Complex(z[j+off1],0.); - i++; - } - else - { - for (octave_idx_type j = 0; j < n; j++) - { - eig_vec(j,i) = - Complex(z[j+off1],z[j+off2]); - if (i < k - 1) - eig_vec(j,i+1) = - Complex(z[j+off1],-z[j+off2]); - } - i+=2; - } - } - } - } + { + octave_idx_type jj = 0; + for (octave_idx_type i = 0; i < k+1; i++) + { + if (dr[i] == 0.0 && di[i] == 0.0 && jj == 0) + jj++; + else + d [i-jj] = Complex (dr[i], di[i]); + } + if (jj == 0 && !rvec) + for (octave_idx_type i = 0; i < k; i++) + d[i] = d[i+1]; + + octave_idx_type k2 = k / 2; + for (octave_idx_type i = 0; i < k2; i++) + { + Complex dtmp = d[i]; + d[i] = d[k - i - 1]; + d[k - i - 1] = dtmp; + } + eig_val.resize(k); + + if (rvec) + { + OCTAVE_LOCAL_BUFFER (double, dtmp, n); + + for (octave_idx_type i = 0; i < k2; i++) + { + octave_idx_type off1 = i * n; + octave_idx_type off2 = (k - i - 1) * n; + + if (off1 == off2) + continue; + + for (octave_idx_type j = 0; j < n; j++) + dtmp[j] = z[off1 + j]; + + for (octave_idx_type j = 0; j < n; j++) + z[off1 + j] = z[off2 + j]; + + for (octave_idx_type j = 0; j < n; j++) + z[off2 + j] = dtmp[j]; + } + + eig_vec.resize (n, k); + octave_idx_type i = 0; + while (i < k) + { + octave_idx_type off1 = i * n; + octave_idx_type off2 = (i+1) * n; + if (std::imag(eig_val(i)) == 0) + { + for (octave_idx_type j = 0; j < n; j++) + eig_vec(j,i) = + Complex(z[j+off1],0.); + i++; + } + else + { + for (octave_idx_type j = 0; j < n; j++) + { + eig_vec(j,i) = + Complex(z[j+off1],z[j+off2]); + if (i < k - 1) + eig_vec(j,i+1) = + Complex(z[j+off1],-z[j+off2]); + } + i+=2; + } + } + } + } else - { - (*current_liboctave_error_handler) - ("eigs: error %d in dneupd", info2); - return -1; - } + { + (*current_liboctave_error_handler) + ("eigs: error %d in dneupd", info2); + return -1; + } } return ip(4); @@ -2398,12 +2398,12 @@ octave_idx_type EigsRealNonSymmetricFunc (EigsFunc fun, octave_idx_type n, - const std::string &_typ, double sigmar, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, ComplexMatrix &eig_vec, - ComplexColumnVector &eig_val, ColumnVector &resid, - std::ostream& os, double tol, int rvec, - bool /* cholB */, int disp, int maxit) + const std::string &_typ, double sigmar, + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, ComplexMatrix &eig_vec, + ComplexColumnVector &eig_val, ColumnVector &resid, + std::ostream& os, double tol, int rvec, + bool /* cholB */, int disp, int maxit) { std::string typ (_typ); bool have_sigma = (sigmar ? true : false); @@ -2423,7 +2423,7 @@ if (n < 3) { (*current_liboctave_error_handler) - ("eigs: n must be at least 3"); + ("eigs: n must be at least 3"); return -1; } @@ -2432,24 +2432,24 @@ p = k * 2 + 1; if (p < 20) - p = 20; + p = 20; if (p > n - 1) - p = n - 1 ; + p = n - 1 ; } if (k <= 0 || k >= n - 1) { (*current_liboctave_error_handler) - ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1).\n" - " Use 'eig(full(A))' instead"); + ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1).\n" + " Use 'eig(full(A))' instead"); return -1; } if (p <= k || p >= n) { (*current_liboctave_error_handler) - ("eigs: opts.p must be greater than k and less than n"); + ("eigs: opts.p must be greater than k and less than n"); return -1; } @@ -2457,24 +2457,24 @@ if (! have_sigma) { if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" && - typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" && - typ != "SI") - (*current_liboctave_error_handler) - ("eigs: unrecognized sigma value"); + typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" && + typ != "SI") + (*current_liboctave_error_handler) + ("eigs: unrecognized sigma value"); if (typ == "LA" || typ == "SA" || typ == "BE") - { - (*current_liboctave_error_handler) - ("eigs: invalid sigma value for unsymmetric problem"); - return -1; - } + { + (*current_liboctave_error_handler) + ("eigs: invalid sigma value for unsymmetric problem"); + return -1; + } if (typ == "SM") - { - typ = "LM"; - sigmar = 0.; - mode = 3; - } + { + typ = "LM"; + sigmar = 0.; + mode = 3; + } } else if (! std::abs (sigmar)) typ = "SM"; @@ -2515,66 +2515,66 @@ do { F77_FUNC (dnaupd, DNAUPD) - (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, - F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), - k, tol, presid, p, v, n, iparam, - ipntr, workd, workl, lwork, info - F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); + (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, + F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), + k, tol, presid, p, v, n, iparam, + ipntr, workd, workl, lwork, info + F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("eigs: unrecoverable exception encountered in dnaupd"); - return -1; - } + { + (*current_liboctave_error_handler) + ("eigs: unrecoverable exception encountered in dnaupd"); + return -1; + } if (disp > 0 && !xisnan(workl[iptr(5)-1])) - { - if (iter++) - { - os << "Iteration " << iter - 1 << - ": a few Ritz values of the " << p << "-by-" << - p << " matrix\n"; - for (int i = 0 ; i < k; i++) - os << " " << workl[iptr(5)+i-1] << "\n"; - } - - // This is a kludge, as ARPACK doesn't give its - // iteration pointer. But as workl[iptr(5)-1] is - // an output value updated at each iteration, setting - // a value in this array to NaN and testing for it - // is a way of obtaining the iteration counter. - if (ido != 99) - workl[iptr(5)-1] = octave_NaN; - } + { + if (iter++) + { + os << "Iteration " << iter - 1 << + ": a few Ritz values of the " << p << "-by-" << + p << " matrix\n"; + for (int i = 0 ; i < k; i++) + os << " " << workl[iptr(5)+i-1] << "\n"; + } + + // This is a kludge, as ARPACK doesn't give its + // iteration pointer. But as workl[iptr(5)-1] is + // an output value updated at each iteration, setting + // a value in this array to NaN and testing for it + // is a way of obtaining the iteration counter. + if (ido != 99) + workl[iptr(5)-1] = octave_NaN; + } if (ido == -1 || ido == 1 || ido == 2) - { - double *ip2 = workd + iptr(0) - 1; - ColumnVector x(n); - - for (octave_idx_type i = 0; i < n; i++) - x(i) = *ip2++; - - ColumnVector y = fun (x, err); - - if (err) - return false; - - ip2 = workd + iptr(1) - 1; - for (octave_idx_type i = 0; i < n; i++) - *ip2++ = y(i); - } + { + double *ip2 = workd + iptr(0) - 1; + ColumnVector x(n); + + for (octave_idx_type i = 0; i < n; i++) + x(i) = *ip2++; + + ColumnVector y = fun (x, err); + + if (err) + return false; + + ip2 = workd + iptr(1) - 1; + for (octave_idx_type i = 0; i < n; i++) + *ip2++ = y(i); + } else - { - if (info < 0) - { - (*current_liboctave_error_handler) - ("eigs: error %d in dsaupd", info); - return -1; - } - break; - } + { + if (info < 0) + { + (*current_liboctave_error_handler) + ("eigs: error %d in dsaupd", info); + return -1; + } + break; + } } while (1); @@ -2609,7 +2609,7 @@ if (f77_exception_encountered) { (*current_liboctave_error_handler) - ("eigs: unrecoverable exception encountered in dneupd"); + ("eigs: unrecoverable exception encountered in dneupd"); return -1; } else @@ -2618,84 +2618,84 @@ Complex *d = eig_val.fortran_vec (); if (info2 == 0) - { - octave_idx_type jj = 0; - for (octave_idx_type i = 0; i < k+1; i++) - { - if (dr[i] == 0.0 && di[i] == 0.0 && jj == 0) - jj++; - else - d [i-jj] = Complex (dr[i], di[i]); - } - if (jj == 0 && !rvec) - for (octave_idx_type i = 0; i < k; i++) - d[i] = d[i+1]; - - octave_idx_type k2 = k / 2; - for (octave_idx_type i = 0; i < k2; i++) - { - Complex dtmp = d[i]; - d[i] = d[k - i - 1]; - d[k - i - 1] = dtmp; - } - eig_val.resize(k); - - if (rvec) - { - OCTAVE_LOCAL_BUFFER (double, dtmp, n); - - for (octave_idx_type i = 0; i < k2; i++) - { - octave_idx_type off1 = i * n; - octave_idx_type off2 = (k - i - 1) * n; - - if (off1 == off2) - continue; - - for (octave_idx_type j = 0; j < n; j++) - dtmp[j] = z[off1 + j]; - - for (octave_idx_type j = 0; j < n; j++) - z[off1 + j] = z[off2 + j]; - - for (octave_idx_type j = 0; j < n; j++) - z[off2 + j] = dtmp[j]; - } - - eig_vec.resize (n, k); - octave_idx_type i = 0; - while (i < k) - { - octave_idx_type off1 = i * n; - octave_idx_type off2 = (i+1) * n; - if (std::imag(eig_val(i)) == 0) - { - for (octave_idx_type j = 0; j < n; j++) - eig_vec(j,i) = - Complex(z[j+off1],0.); - i++; - } - else - { - for (octave_idx_type j = 0; j < n; j++) - { - eig_vec(j,i) = - Complex(z[j+off1],z[j+off2]); - if (i < k - 1) - eig_vec(j,i+1) = - Complex(z[j+off1],-z[j+off2]); - } - i+=2; - } - } - } - } + { + octave_idx_type jj = 0; + for (octave_idx_type i = 0; i < k+1; i++) + { + if (dr[i] == 0.0 && di[i] == 0.0 && jj == 0) + jj++; + else + d [i-jj] = Complex (dr[i], di[i]); + } + if (jj == 0 && !rvec) + for (octave_idx_type i = 0; i < k; i++) + d[i] = d[i+1]; + + octave_idx_type k2 = k / 2; + for (octave_idx_type i = 0; i < k2; i++) + { + Complex dtmp = d[i]; + d[i] = d[k - i - 1]; + d[k - i - 1] = dtmp; + } + eig_val.resize(k); + + if (rvec) + { + OCTAVE_LOCAL_BUFFER (double, dtmp, n); + + for (octave_idx_type i = 0; i < k2; i++) + { + octave_idx_type off1 = i * n; + octave_idx_type off2 = (k - i - 1) * n; + + if (off1 == off2) + continue; + + for (octave_idx_type j = 0; j < n; j++) + dtmp[j] = z[off1 + j]; + + for (octave_idx_type j = 0; j < n; j++) + z[off1 + j] = z[off2 + j]; + + for (octave_idx_type j = 0; j < n; j++) + z[off2 + j] = dtmp[j]; + } + + eig_vec.resize (n, k); + octave_idx_type i = 0; + while (i < k) + { + octave_idx_type off1 = i * n; + octave_idx_type off2 = (i+1) * n; + if (std::imag(eig_val(i)) == 0) + { + for (octave_idx_type j = 0; j < n; j++) + eig_vec(j,i) = + Complex(z[j+off1],0.); + i++; + } + else + { + for (octave_idx_type j = 0; j < n; j++) + { + eig_vec(j,i) = + Complex(z[j+off1],z[j+off2]); + if (i < k - 1) + eig_vec(j,i+1) = + Complex(z[j+off1],-z[j+off2]); + } + i+=2; + } + } + } + } else - { - (*current_liboctave_error_handler) - ("eigs: error %d in dneupd", info2); - return -1; - } + { + (*current_liboctave_error_handler) + ("eigs: error %d in dneupd", info2); + return -1; + } } return ip(4); @@ -2704,13 +2704,13 @@ template <class M> octave_idx_type EigsComplexNonSymmetricMatrix (const M& m, const std::string typ, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, ComplexMatrix &eig_vec, - ComplexColumnVector &eig_val, const M& _b, - ColumnVector &permB, - ComplexColumnVector &cresid, - std::ostream& os, double tol, int rvec, - bool cholB, int disp, int maxit) + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, ComplexMatrix &eig_vec, + ComplexColumnVector &eig_val, const M& _b, + ColumnVector &permB, + ComplexColumnVector &cresid, + std::ostream& os, double tol, int rvec, + bool cholB, int disp, int maxit) { M b(_b); octave_idx_type n = m.cols (); @@ -2729,7 +2729,7 @@ if (have_b && (m.rows() != b.rows() || m.rows() != b.cols())) { (*current_liboctave_error_handler) - ("eigs: B must be square and the same size as A"); + ("eigs: B must be square and the same size as A"); return -1; } @@ -2741,14 +2741,14 @@ Array<double> ri (octave_rand::vector(n)); cresid = ComplexColumnVector (n); for (octave_idx_type i = 0; i < n; i++) - cresid(i) = Complex(rr(i),ri(i)); + cresid(i) = Complex(rr(i),ri(i)); octave_rand::distribution(rand_dist); } if (n < 3) { (*current_liboctave_error_handler) - ("eigs: n must be at least 3"); + ("eigs: n must be at least 3"); return -1; } @@ -2757,24 +2757,24 @@ p = k * 2 + 1; if (p < 20) - p = 20; + p = 20; if (p > n - 1) - p = n - 1 ; + p = n - 1 ; } if (k <= 0 || k >= n - 1) { (*current_liboctave_error_handler) - ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1).\n" - " Use 'eig(full(A))' instead"); + ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1).\n" + " Use 'eig(full(A))' instead"); return -1; } if (p <= k || p >= n) { (*current_liboctave_error_handler) - ("eigs: opts.p must be greater than k and less than n"); + ("eigs: opts.p must be greater than k and less than n"); return -1; } @@ -2782,27 +2782,27 @@ { // Check the we really have a permutation vector if (permB.length() != n) - { - (*current_liboctave_error_handler) - ("eigs: permB vector invalid"); - return -1; - } + { + (*current_liboctave_error_handler) + ("eigs: permB vector invalid"); + return -1; + } else - { - Array<bool> checked(n,false); - for (octave_idx_type i = 0; i < n; i++) - { - octave_idx_type bidx = - static_cast<octave_idx_type> (permB(i)); - if (checked(bidx) || bidx < 0 || - bidx >= n || D_NINT (bidx) != bidx) - { - (*current_liboctave_error_handler) - ("eigs: permB vector invalid"); - return -1; - } - } - } + { + Array<bool> checked(n,false); + for (octave_idx_type i = 0; i < n; i++) + { + octave_idx_type bidx = + static_cast<octave_idx_type> (permB(i)); + if (checked(bidx) || bidx < 0 || + bidx >= n || D_NINT (bidx) != bidx) + { + (*current_liboctave_error_handler) + ("eigs: permB vector invalid"); + return -1; + } + } + } } if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" && @@ -2810,14 +2810,14 @@ typ != "SI") { (*current_liboctave_error_handler) - ("eigs: unrecognized sigma value"); + ("eigs: unrecognized sigma value"); return -1; } if (typ == "LA" || typ == "SA" || typ == "BE") { (*current_liboctave_error_handler) - ("eigs: invalid sigma value for complex problem"); + ("eigs: invalid sigma value for complex problem"); return -1; } @@ -2826,25 +2826,25 @@ // See Note 3 dsaupd note3 = true; if (cholB) - { - bt = b; - b = b.hermitian(); - if (permB.length() == 0) - { - permB = ColumnVector(n); - for (octave_idx_type i = 0; i < n; i++) - permB(i) = i; - } - } + { + bt = b; + b = b.hermitian(); + if (permB.length() == 0) + { + permB = ColumnVector(n); + for (octave_idx_type i = 0; i < n; i++) + permB(i) = i; + } + } else - { - if (! make_cholb(b, bt, permB)) - { - (*current_liboctave_error_handler) - ("eigs: The matrix B is not positive definite"); - return -1; - } - } + { + if (! make_cholb(b, bt, permB)) + { + (*current_liboctave_error_handler) + ("eigs: The matrix B is not positive definite"); + return -1; + } + } } Array<octave_idx_type> ip (11); @@ -2869,7 +2869,7 @@ octave_idx_type ido = 0; int iter = 0; octave_idx_type lwork = p * (3 * p + 5); - + OCTAVE_LOCAL_BUFFER (Complex, v, n * p); OCTAVE_LOCAL_BUFFER (Complex, workl, lwork); OCTAVE_LOCAL_BUFFER (Complex, workd, 3 * n); @@ -2879,65 +2879,65 @@ do { F77_FUNC (znaupd, ZNAUPD) - (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, - F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), - k, tol, presid, p, v, n, iparam, - ipntr, workd, workl, lwork, rwork, info - F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); + (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, + F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), + k, tol, presid, p, v, n, iparam, + ipntr, workd, workl, lwork, rwork, info + F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("eigs: unrecoverable exception encountered in znaupd"); - return -1; - } + { + (*current_liboctave_error_handler) + ("eigs: unrecoverable exception encountered in znaupd"); + return -1; + } if (disp > 0 && !xisnan(workl[iptr(5)-1])) - { - if (iter++) - { - os << "Iteration " << iter - 1 << - ": a few Ritz values of the " << p << "-by-" << - p << " matrix\n"; - for (int i = 0 ; i < k; i++) - os << " " << workl[iptr(5)+i-1] << "\n"; - } - - // This is a kludge, as ARPACK doesn't give its - // iteration pointer. But as workl[iptr(5)-1] is - // an output value updated at each iteration, setting - // a value in this array to NaN and testing for it - // is a way of obtaining the iteration counter. - if (ido != 99) - workl[iptr(5)-1] = octave_NaN; - } + { + if (iter++) + { + os << "Iteration " << iter - 1 << + ": a few Ritz values of the " << p << "-by-" << + p << " matrix\n"; + for (int i = 0 ; i < k; i++) + os << " " << workl[iptr(5)+i-1] << "\n"; + } + + // This is a kludge, as ARPACK doesn't give its + // iteration pointer. But as workl[iptr(5)-1] is + // an output value updated at each iteration, setting + // a value in this array to NaN and testing for it + // is a way of obtaining the iteration counter. + if (ido != 99) + workl[iptr(5)-1] = octave_NaN; + } if (ido == -1 || ido == 1 || ido == 2) - { - if (have_b) - { - ComplexMatrix mtmp (n,1); - for (octave_idx_type i = 0; i < n; i++) - mtmp(i,0) = workd[i + iptr(0) - 1]; - mtmp = utsolve(bt, permB, m * ltsolve(b, permB, mtmp)); - for (octave_idx_type i = 0; i < n; i++) - workd[i+iptr(1)-1] = mtmp(i,0); - - } - else if (!vector_product (m, workd + iptr(0) - 1, - workd + iptr(1) - 1)) - break; - } + { + if (have_b) + { + ComplexMatrix mtmp (n,1); + for (octave_idx_type i = 0; i < n; i++) + mtmp(i,0) = workd[i + iptr(0) - 1]; + mtmp = utsolve(bt, permB, m * ltsolve(b, permB, mtmp)); + for (octave_idx_type i = 0; i < n; i++) + workd[i+iptr(1)-1] = mtmp(i,0); + + } + else if (!vector_product (m, workd + iptr(0) - 1, + workd + iptr(1) - 1)) + break; + } else - { - if (info < 0) - { - (*current_liboctave_error_handler) - ("eigs: error %d in znaupd", info); - return -1; - } - break; - } + { + if (info < 0) + { + (*current_liboctave_error_handler) + ("eigs: error %d in znaupd", info); + return -1; + } + break; + } } while (1); @@ -2971,7 +2971,7 @@ if (f77_exception_encountered) { (*current_liboctave_error_handler) - ("eigs: unrecoverable exception encountered in zneupd"); + ("eigs: unrecoverable exception encountered in zneupd"); return -1; } @@ -2979,43 +2979,43 @@ { octave_idx_type k2 = k / 2; for (octave_idx_type i = 0; i < k2; i++) - { - Complex ctmp = d[i]; - d[i] = d[k - i - 1]; - d[k - i - 1] = ctmp; - } + { + Complex ctmp = d[i]; + d[i] = d[k - i - 1]; + d[k - i - 1] = ctmp; + } eig_val.resize(k); if (rvec) - { - OCTAVE_LOCAL_BUFFER (Complex, ctmp, n); - - for (octave_idx_type i = 0; i < k2; i++) - { - octave_idx_type off1 = i * n; - octave_idx_type off2 = (k - i - 1) * n; - - if (off1 == off2) - continue; - - for (octave_idx_type j = 0; j < n; j++) - ctmp[j] = z[off1 + j]; - - for (octave_idx_type j = 0; j < n; j++) - z[off1 + j] = z[off2 + j]; - - for (octave_idx_type j = 0; j < n; j++) - z[off2 + j] = ctmp[j]; - } - - if (note3) - eig_vec = ltsolve(b, permB, eig_vec); - } + { + OCTAVE_LOCAL_BUFFER (Complex, ctmp, n); + + for (octave_idx_type i = 0; i < k2; i++) + { + octave_idx_type off1 = i * n; + octave_idx_type off2 = (k - i - 1) * n; + + if (off1 == off2) + continue; + + for (octave_idx_type j = 0; j < n; j++) + ctmp[j] = z[off1 + j]; + + for (octave_idx_type j = 0; j < n; j++) + z[off1 + j] = z[off2 + j]; + + for (octave_idx_type j = 0; j < n; j++) + z[off2 + j] = ctmp[j]; + } + + if (note3) + eig_vec = ltsolve(b, permB, eig_vec); + } } else { (*current_liboctave_error_handler) - ("eigs: error %d in zneupd", info2); + ("eigs: error %d in zneupd", info2); return -1; } @@ -3025,14 +3025,14 @@ template <class M> octave_idx_type EigsComplexNonSymmetricMatrixShift (const M& m, Complex sigma, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, - ComplexMatrix &eig_vec, - ComplexColumnVector &eig_val, const M& _b, - ColumnVector &permB, - ComplexColumnVector &cresid, - std::ostream& os, double tol, int rvec, - bool cholB, int disp, int maxit) + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, + ComplexMatrix &eig_vec, + ComplexColumnVector &eig_val, const M& _b, + ColumnVector &permB, + ComplexColumnVector &cresid, + std::ostream& os, double tol, int rvec, + bool cholB, int disp, int maxit) { M b(_b); octave_idx_type n = m.cols (); @@ -3048,15 +3048,15 @@ if (have_b && (m.rows() != b.rows() || m.rows() != b.cols())) { (*current_liboctave_error_handler) - ("eigs: B must be square and the same size as A"); + ("eigs: B must be square and the same size as A"); return -1; } // FIXME: The "SM" type for mode 1 seems unstable though faster!! //if (! std::abs (sigma)) // return EigsComplexNonSymmetricMatrix (m, "SM", k, p, info, eig_vec, - // eig_val, _b, permB, cresid, os, tol, - // rvec, cholB, disp, maxit); + // eig_val, _b, permB, cresid, os, tol, + // rvec, cholB, disp, maxit); if (cresid.is_empty()) { @@ -3066,14 +3066,14 @@ Array<double> ri (octave_rand::vector(n)); cresid = ComplexColumnVector (n); for (octave_idx_type i = 0; i < n; i++) - cresid(i) = Complex(rr(i),ri(i)); + cresid(i) = Complex(rr(i),ri(i)); octave_rand::distribution(rand_dist); } if (n < 3) { (*current_liboctave_error_handler) - ("eigs: n must be at least 3"); + ("eigs: n must be at least 3"); return -1; } @@ -3082,24 +3082,24 @@ p = k * 2 + 1; if (p < 20) - p = 20; + p = 20; if (p > n - 1) - p = n - 1 ; + p = n - 1 ; } if (k <= 0 || k >= n - 1) { (*current_liboctave_error_handler) - ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1).\n" - " Use 'eig(full(A))' instead"); + ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1).\n" + " Use 'eig(full(A))' instead"); return -1; } if (p <= k || p >= n) { (*current_liboctave_error_handler) - ("eigs: opts.p must be greater than k and less than n"); + ("eigs: opts.p must be greater than k and less than n"); return -1; } @@ -3107,26 +3107,26 @@ { // Check that we really have a permutation vector if (permB.length() != n) - { - (*current_liboctave_error_handler) ("eigs: permB vector invalid"); - return -1; - } + { + (*current_liboctave_error_handler) ("eigs: permB vector invalid"); + return -1; + } else - { - Array<bool> checked(n,false); - for (octave_idx_type i = 0; i < n; i++) - { - octave_idx_type bidx = - static_cast<octave_idx_type> (permB(i)); - if (checked(bidx) || bidx < 0 || - bidx >= n || D_NINT (bidx) != bidx) - { - (*current_liboctave_error_handler) - ("eigs: permB vector invalid"); - return -1; - } - } - } + { + Array<bool> checked(n,false); + for (octave_idx_type i = 0; i < n; i++) + { + octave_idx_type bidx = + static_cast<octave_idx_type> (permB(i)); + if (checked(bidx) || bidx < 0 || + bidx >= n || D_NINT (bidx) != bidx) + { + (*current_liboctave_error_handler) + ("eigs: permB vector invalid"); + return -1; + } + } + } } char bmat = 'I'; @@ -3163,7 +3163,7 @@ return -1; octave_idx_type lwork = p * (3 * p + 5); - + OCTAVE_LOCAL_BUFFER (Complex, v, n * p); OCTAVE_LOCAL_BUFFER (Complex, workl, lwork); OCTAVE_LOCAL_BUFFER (Complex, workd, 3 * n); @@ -3173,111 +3173,111 @@ do { F77_FUNC (znaupd, ZNAUPD) - (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, - F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), - k, tol, presid, p, v, n, iparam, - ipntr, workd, workl, lwork, rwork, info - F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); + (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, + F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), + k, tol, presid, p, v, n, iparam, + ipntr, workd, workl, lwork, rwork, info + F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("eigs: unrecoverable exception encountered in znaupd"); - return -1; - } + { + (*current_liboctave_error_handler) + ("eigs: unrecoverable exception encountered in znaupd"); + return -1; + } if (disp > 0 && !xisnan(workl[iptr(5)-1])) - { - if (iter++) - { - os << "Iteration " << iter - 1 << - ": a few Ritz values of the " << p << "-by-" << - p << " matrix\n"; - for (int i = 0 ; i < k; i++) - os << " " << workl[iptr(5)+i-1] << "\n"; - } - - // This is a kludge, as ARPACK doesn't give its - // iteration pointer. But as workl[iptr(5)-1] is - // an output value updated at each iteration, setting - // a value in this array to NaN and testing for it - // is a way of obtaining the iteration counter. - if (ido != 99) - workl[iptr(5)-1] = octave_NaN; - } + { + if (iter++) + { + os << "Iteration " << iter - 1 << + ": a few Ritz values of the " << p << "-by-" << + p << " matrix\n"; + for (int i = 0 ; i < k; i++) + os << " " << workl[iptr(5)+i-1] << "\n"; + } + + // This is a kludge, as ARPACK doesn't give its + // iteration pointer. But as workl[iptr(5)-1] is + // an output value updated at each iteration, setting + // a value in this array to NaN and testing for it + // is a way of obtaining the iteration counter. + if (ido != 99) + workl[iptr(5)-1] = octave_NaN; + } if (ido == -1 || ido == 1 || ido == 2) - { - if (have_b) - { - if (ido == -1) - { - OCTAVE_LOCAL_BUFFER (Complex, ctmp, n); - - vector_product (m, workd+iptr(0)-1, ctmp); - - ComplexMatrix tmp(n, 1); - - for (octave_idx_type i = 0; i < n; i++) - tmp(i,0) = ctmp[P[i]]; - - lusolve (L, U, tmp); - - Complex *ip2 = workd+iptr(1)-1; - for (octave_idx_type i = 0; i < n; i++) - ip2[Q[i]] = tmp(i,0); - } - else if (ido == 2) - vector_product (b, workd + iptr(0) - 1, workd + iptr(1) - 1); - else - { - Complex *ip2 = workd+iptr(2)-1; - ComplexMatrix tmp(n, 1); - - for (octave_idx_type i = 0; i < n; i++) - tmp(i,0) = ip2[P[i]]; - - lusolve (L, U, tmp); - - ip2 = workd+iptr(1)-1; - for (octave_idx_type i = 0; i < n; i++) - ip2[Q[i]] = tmp(i,0); - } - } - else - { - if (ido == 2) - { - for (octave_idx_type i = 0; i < n; i++) - workd[iptr(0) + i - 1] = - workd[iptr(1) + i - 1]; - } - else - { - Complex *ip2 = workd+iptr(0)-1; - ComplexMatrix tmp(n, 1); - - for (octave_idx_type i = 0; i < n; i++) - tmp(i,0) = ip2[P[i]]; - - lusolve (L, U, tmp); - - ip2 = workd+iptr(1)-1; - for (octave_idx_type i = 0; i < n; i++) - ip2[Q[i]] = tmp(i,0); - } - } - } + { + if (have_b) + { + if (ido == -1) + { + OCTAVE_LOCAL_BUFFER (Complex, ctmp, n); + + vector_product (m, workd+iptr(0)-1, ctmp); + + ComplexMatrix tmp(n, 1); + + for (octave_idx_type i = 0; i < n; i++) + tmp(i,0) = ctmp[P[i]]; + + lusolve (L, U, tmp); + + Complex *ip2 = workd+iptr(1)-1; + for (octave_idx_type i = 0; i < n; i++) + ip2[Q[i]] = tmp(i,0); + } + else if (ido == 2) + vector_product (b, workd + iptr(0) - 1, workd + iptr(1) - 1); + else + { + Complex *ip2 = workd+iptr(2)-1; + ComplexMatrix tmp(n, 1); + + for (octave_idx_type i = 0; i < n; i++) + tmp(i,0) = ip2[P[i]]; + + lusolve (L, U, tmp); + + ip2 = workd+iptr(1)-1; + for (octave_idx_type i = 0; i < n; i++) + ip2[Q[i]] = tmp(i,0); + } + } + else + { + if (ido == 2) + { + for (octave_idx_type i = 0; i < n; i++) + workd[iptr(0) + i - 1] = + workd[iptr(1) + i - 1]; + } + else + { + Complex *ip2 = workd+iptr(0)-1; + ComplexMatrix tmp(n, 1); + + for (octave_idx_type i = 0; i < n; i++) + tmp(i,0) = ip2[P[i]]; + + lusolve (L, U, tmp); + + ip2 = workd+iptr(1)-1; + for (octave_idx_type i = 0; i < n; i++) + ip2[Q[i]] = tmp(i,0); + } + } + } else - { - if (info < 0) - { - (*current_liboctave_error_handler) - ("eigs: error %d in dsaupd", info); - return -1; - } - break; - } + { + if (info < 0) + { + (*current_liboctave_error_handler) + ("eigs: error %d in dsaupd", info); + return -1; + } + break; + } } while (1); @@ -3311,7 +3311,7 @@ if (f77_exception_encountered) { (*current_liboctave_error_handler) - ("eigs: unrecoverable exception encountered in zneupd"); + ("eigs: unrecoverable exception encountered in zneupd"); return -1; } @@ -3319,40 +3319,40 @@ { octave_idx_type k2 = k / 2; for (octave_idx_type i = 0; i < k2; i++) - { - Complex ctmp = d[i]; - d[i] = d[k - i - 1]; - d[k - i - 1] = ctmp; - } + { + Complex ctmp = d[i]; + d[i] = d[k - i - 1]; + d[k - i - 1] = ctmp; + } eig_val.resize(k); if (rvec) - { - OCTAVE_LOCAL_BUFFER (Complex, ctmp, n); - - for (octave_idx_type i = 0; i < k2; i++) - { - octave_idx_type off1 = i * n; - octave_idx_type off2 = (k - i - 1) * n; - - if (off1 == off2) - continue; - - for (octave_idx_type j = 0; j < n; j++) - ctmp[j] = z[off1 + j]; - - for (octave_idx_type j = 0; j < n; j++) - z[off1 + j] = z[off2 + j]; - - for (octave_idx_type j = 0; j < n; j++) - z[off2 + j] = ctmp[j]; - } - } + { + OCTAVE_LOCAL_BUFFER (Complex, ctmp, n); + + for (octave_idx_type i = 0; i < k2; i++) + { + octave_idx_type off1 = i * n; + octave_idx_type off2 = (k - i - 1) * n; + + if (off1 == off2) + continue; + + for (octave_idx_type j = 0; j < n; j++) + ctmp[j] = z[off1 + j]; + + for (octave_idx_type j = 0; j < n; j++) + z[off1 + j] = z[off2 + j]; + + for (octave_idx_type j = 0; j < n; j++) + z[off2 + j] = ctmp[j]; + } + } } else { (*current_liboctave_error_handler) - ("eigs: error %d in zneupd", info2); + ("eigs: error %d in zneupd", info2); return -1; } @@ -3361,13 +3361,13 @@ octave_idx_type EigsComplexNonSymmetricFunc (EigsComplexFunc fun, octave_idx_type n, - const std::string &_typ, Complex sigma, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, ComplexMatrix &eig_vec, - ComplexColumnVector &eig_val, - ComplexColumnVector &cresid, std::ostream& os, - double tol, int rvec, bool /* cholB */, - int disp, int maxit) + const std::string &_typ, Complex sigma, + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, ComplexMatrix &eig_vec, + ComplexColumnVector &eig_val, + ComplexColumnVector &cresid, std::ostream& os, + double tol, int rvec, bool /* cholB */, + int disp, int maxit) { std::string typ (_typ); bool have_sigma = (std::abs(sigma) ? true : false); @@ -3383,14 +3383,14 @@ Array<double> ri (octave_rand::vector(n)); cresid = ComplexColumnVector (n); for (octave_idx_type i = 0; i < n; i++) - cresid(i) = Complex(rr(i),ri(i)); + cresid(i) = Complex(rr(i),ri(i)); octave_rand::distribution(rand_dist); } if (n < 3) { (*current_liboctave_error_handler) - ("eigs: n must be at least 3"); + ("eigs: n must be at least 3"); return -1; } @@ -3399,48 +3399,48 @@ p = k * 2 + 1; if (p < 20) - p = 20; + p = 20; if (p > n - 1) - p = n - 1 ; + p = n - 1 ; } if (k <= 0 || k >= n - 1) { (*current_liboctave_error_handler) - ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1).\n" - " Use 'eig(full(A))' instead"); + ("eigs: Invalid number of eigenvalues to extract (must be 0 < k < n-1).\n" + " Use 'eig(full(A))' instead"); return -1; } if (p <= k || p >= n) { (*current_liboctave_error_handler) - ("eigs: opts.p must be greater than k and less than n"); + ("eigs: opts.p must be greater than k and less than n"); return -1; } if (! have_sigma) { if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" && - typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" && - typ != "SI") - (*current_liboctave_error_handler) - ("eigs: unrecognized sigma value"); + typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" && + typ != "SI") + (*current_liboctave_error_handler) + ("eigs: unrecognized sigma value"); if (typ == "LA" || typ == "SA" || typ == "BE") - { - (*current_liboctave_error_handler) - ("eigs: invalid sigma value for complex problem"); - return -1; - } + { + (*current_liboctave_error_handler) + ("eigs: invalid sigma value for complex problem"); + return -1; + } if (typ == "SM") - { - typ = "LM"; - sigma = 0.; - mode = 3; - } + { + typ = "LM"; + sigma = 0.; + mode = 3; + } } else if (! std::abs (sigma)) typ = "SM"; @@ -3472,7 +3472,7 @@ octave_idx_type ido = 0; int iter = 0; octave_idx_type lwork = p * (3 * p + 5); - + OCTAVE_LOCAL_BUFFER (Complex, v, n * p); OCTAVE_LOCAL_BUFFER (Complex, workl, lwork); OCTAVE_LOCAL_BUFFER (Complex, workd, 3 * n); @@ -3482,66 +3482,66 @@ do { F77_FUNC (znaupd, ZNAUPD) - (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, - F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), - k, tol, presid, p, v, n, iparam, - ipntr, workd, workl, lwork, rwork, info - F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); + (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, + F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), + k, tol, presid, p, v, n, iparam, + ipntr, workd, workl, lwork, rwork, info + F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("eigs: unrecoverable exception encountered in znaupd"); - return -1; - } + { + (*current_liboctave_error_handler) + ("eigs: unrecoverable exception encountered in znaupd"); + return -1; + } if (disp > 0 && !xisnan(workl[iptr(5)-1])) - { - if (iter++) - { - os << "Iteration " << iter - 1 << - ": a few Ritz values of the " << p << "-by-" << - p << " matrix\n"; - for (int i = 0 ; i < k; i++) - os << " " << workl[iptr(5)+i-1] << "\n"; - } - - // This is a kludge, as ARPACK doesn't give its - // iteration pointer. But as workl[iptr(5)-1] is - // an output value updated at each iteration, setting - // a value in this array to NaN and testing for it - // is a way of obtaining the iteration counter. - if (ido != 99) - workl[iptr(5)-1] = octave_NaN; - } + { + if (iter++) + { + os << "Iteration " << iter - 1 << + ": a few Ritz values of the " << p << "-by-" << + p << " matrix\n"; + for (int i = 0 ; i < k; i++) + os << " " << workl[iptr(5)+i-1] << "\n"; + } + + // This is a kludge, as ARPACK doesn't give its + // iteration pointer. But as workl[iptr(5)-1] is + // an output value updated at each iteration, setting + // a value in this array to NaN and testing for it + // is a way of obtaining the iteration counter. + if (ido != 99) + workl[iptr(5)-1] = octave_NaN; + } if (ido == -1 || ido == 1 || ido == 2) - { - Complex *ip2 = workd + iptr(0) - 1; - ComplexColumnVector x(n); - - for (octave_idx_type i = 0; i < n; i++) - x(i) = *ip2++; - - ComplexColumnVector y = fun (x, err); - - if (err) - return false; - - ip2 = workd + iptr(1) - 1; - for (octave_idx_type i = 0; i < n; i++) - *ip2++ = y(i); - } + { + Complex *ip2 = workd + iptr(0) - 1; + ComplexColumnVector x(n); + + for (octave_idx_type i = 0; i < n; i++) + x(i) = *ip2++; + + ComplexColumnVector y = fun (x, err); + + if (err) + return false; + + ip2 = workd + iptr(1) - 1; + for (octave_idx_type i = 0; i < n; i++) + *ip2++ = y(i); + } else - { - if (info < 0) - { - (*current_liboctave_error_handler) - ("eigs: error %d in dsaupd", info); - return -1; - } - break; - } + { + if (info < 0) + { + (*current_liboctave_error_handler) + ("eigs: error %d in dsaupd", info); + return -1; + } + break; + } } while (1); @@ -3575,7 +3575,7 @@ if (f77_exception_encountered) { (*current_liboctave_error_handler) - ("eigs: unrecoverable exception encountered in zneupd"); + ("eigs: unrecoverable exception encountered in zneupd"); return -1; } @@ -3583,40 +3583,40 @@ { octave_idx_type k2 = k / 2; for (octave_idx_type i = 0; i < k2; i++) - { - Complex ctmp = d[i]; - d[i] = d[k - i - 1]; - d[k - i - 1] = ctmp; - } + { + Complex ctmp = d[i]; + d[i] = d[k - i - 1]; + d[k - i - 1] = ctmp; + } eig_val.resize(k); if (rvec) - { - OCTAVE_LOCAL_BUFFER (Complex, ctmp, n); - - for (octave_idx_type i = 0; i < k2; i++) - { - octave_idx_type off1 = i * n; - octave_idx_type off2 = (k - i - 1) * n; - - if (off1 == off2) - continue; - - for (octave_idx_type j = 0; j < n; j++) - ctmp[j] = z[off1 + j]; - - for (octave_idx_type j = 0; j < n; j++) - z[off1 + j] = z[off2 + j]; - - for (octave_idx_type j = 0; j < n; j++) - z[off2 + j] = ctmp[j]; - } - } + { + OCTAVE_LOCAL_BUFFER (Complex, ctmp, n); + + for (octave_idx_type i = 0; i < k2; i++) + { + octave_idx_type off1 = i * n; + octave_idx_type off2 = (k - i - 1) * n; + + if (off1 == off2) + continue; + + for (octave_idx_type j = 0; j < n; j++) + ctmp[j] = z[off1 + j]; + + for (octave_idx_type j = 0; j < n; j++) + z[off1 + j] = z[off2 + j]; + + for (octave_idx_type j = 0; j < n; j++) + z[off2 + j] = ctmp[j]; + } + } } else { (*current_liboctave_error_handler) - ("eigs: error %d in zneupd", info2); + ("eigs: error %d in zneupd", info2); return -1; } @@ -3626,168 +3626,168 @@ #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL) extern octave_idx_type EigsRealSymmetricMatrix (const Matrix& m, const std::string typ, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, Matrix &eig_vec, - ColumnVector &eig_val, const Matrix& b, - ColumnVector &permB, ColumnVector &resid, - std::ostream &os, double tol = DBL_EPSILON, - int rvec = 0, bool cholB = 0, int disp = 0, - int maxit = 300); + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, Matrix &eig_vec, + ColumnVector &eig_val, const Matrix& b, + ColumnVector &permB, ColumnVector &resid, + std::ostream &os, double tol = DBL_EPSILON, + int rvec = 0, bool cholB = 0, int disp = 0, + int maxit = 300); extern octave_idx_type EigsRealSymmetricMatrix (const SparseMatrix& m, const std::string typ, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, Matrix &eig_vec, - ColumnVector &eig_val, const SparseMatrix& b, - ColumnVector &permB, ColumnVector &resid, - std::ostream& os, double tol = DBL_EPSILON, - int rvec = 0, bool cholB = 0, int disp = 0, - int maxit = 300); + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, Matrix &eig_vec, + ColumnVector &eig_val, const SparseMatrix& b, + ColumnVector &permB, ColumnVector &resid, + std::ostream& os, double tol = DBL_EPSILON, + int rvec = 0, bool cholB = 0, int disp = 0, + int maxit = 300); extern octave_idx_type EigsRealSymmetricMatrixShift (const Matrix& m, double sigma, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, Matrix &eig_vec, - ColumnVector &eig_val, const Matrix& b, - ColumnVector &permB, ColumnVector &resid, - std::ostream &os, double tol = DBL_EPSILON, - int rvec = 0, bool cholB = 0, int disp = 0, - int maxit = 300); + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, Matrix &eig_vec, + ColumnVector &eig_val, const Matrix& b, + ColumnVector &permB, ColumnVector &resid, + std::ostream &os, double tol = DBL_EPSILON, + int rvec = 0, bool cholB = 0, int disp = 0, + int maxit = 300); extern octave_idx_type EigsRealSymmetricMatrixShift (const SparseMatrix& m, double sigma, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, Matrix &eig_vec, - ColumnVector &eig_val, const SparseMatrix& b, - ColumnVector &permB, ColumnVector &resid, - std::ostream &os, double tol = DBL_EPSILON, - int rvec = 0, bool cholB = 0, int disp = 0, - int maxit = 300); + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, Matrix &eig_vec, + ColumnVector &eig_val, const SparseMatrix& b, + ColumnVector &permB, ColumnVector &resid, + std::ostream &os, double tol = DBL_EPSILON, + int rvec = 0, bool cholB = 0, int disp = 0, + int maxit = 300); extern octave_idx_type EigsRealSymmetricFunc (EigsFunc fun, octave_idx_type n, - const std::string &typ, double sigma, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, - Matrix &eig_vec, ColumnVector &eig_val, - ColumnVector &resid, std::ostream &os, - double tol = DBL_EPSILON, int rvec = 0, - bool cholB = 0, int disp = 0, int maxit = 300); + const std::string &typ, double sigma, + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, + Matrix &eig_vec, ColumnVector &eig_val, + ColumnVector &resid, std::ostream &os, + double tol = DBL_EPSILON, int rvec = 0, + bool cholB = 0, int disp = 0, int maxit = 300); extern octave_idx_type EigsRealNonSymmetricMatrix (const Matrix& m, const std::string typ, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, ComplexMatrix &eig_vec, - ComplexColumnVector &eig_val, const Matrix& b, - ColumnVector &permB, ColumnVector &resid, - std::ostream &os, double tol = DBL_EPSILON, - int rvec = 0, bool cholB = 0, int disp = 0, - int maxit = 300); + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, ComplexMatrix &eig_vec, + ComplexColumnVector &eig_val, const Matrix& b, + ColumnVector &permB, ColumnVector &resid, + std::ostream &os, double tol = DBL_EPSILON, + int rvec = 0, bool cholB = 0, int disp = 0, + int maxit = 300); extern octave_idx_type EigsRealNonSymmetricMatrix (const SparseMatrix& m, const std::string typ, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, ComplexMatrix &eig_vec, - ComplexColumnVector &eig_val, - const SparseMatrix& b, - ColumnVector &permB, ColumnVector &resid, - std::ostream &os, double tol = DBL_EPSILON, - int rvec = 0, bool cholB = 0, int disp = 0, - int maxit = 300); + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, ComplexMatrix &eig_vec, + ComplexColumnVector &eig_val, + const SparseMatrix& b, + ColumnVector &permB, ColumnVector &resid, + std::ostream &os, double tol = DBL_EPSILON, + int rvec = 0, bool cholB = 0, int disp = 0, + int maxit = 300); extern octave_idx_type EigsRealNonSymmetricMatrixShift (const Matrix& m, double sigma, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, - ComplexMatrix &eig_vec, - ComplexColumnVector &eig_val, const Matrix& b, - ColumnVector &permB, ColumnVector &resid, - std::ostream &os, double tol = DBL_EPSILON, - int rvec = 0, bool cholB = 0, int disp = 0, - int maxit = 300); + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, + ComplexMatrix &eig_vec, + ComplexColumnVector &eig_val, const Matrix& b, + ColumnVector &permB, ColumnVector &resid, + std::ostream &os, double tol = DBL_EPSILON, + int rvec = 0, bool cholB = 0, int disp = 0, + int maxit = 300); extern octave_idx_type EigsRealNonSymmetricMatrixShift (const SparseMatrix& m, double sigma, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, - ComplexMatrix &eig_vec, - ComplexColumnVector &eig_val, - const SparseMatrix& b, - ColumnVector &permB, ColumnVector &resid, - std::ostream &os, double tol = DBL_EPSILON, - int rvec = 0, bool cholB = 0, int disp = 0, - int maxit = 300); + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, + ComplexMatrix &eig_vec, + ComplexColumnVector &eig_val, + const SparseMatrix& b, + ColumnVector &permB, ColumnVector &resid, + std::ostream &os, double tol = DBL_EPSILON, + int rvec = 0, bool cholB = 0, int disp = 0, + int maxit = 300); extern octave_idx_type EigsRealNonSymmetricFunc (EigsFunc fun, octave_idx_type n, - const std::string &_typ, double sigma, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, ComplexMatrix &eig_vec, - ComplexColumnVector &eig_val, - ColumnVector &resid, std::ostream& os, - double tol = DBL_EPSILON, int rvec = 0, - bool cholB = 0, int disp = 0, int maxit = 300); + const std::string &_typ, double sigma, + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, ComplexMatrix &eig_vec, + ComplexColumnVector &eig_val, + ColumnVector &resid, std::ostream& os, + double tol = DBL_EPSILON, int rvec = 0, + bool cholB = 0, int disp = 0, int maxit = 300); extern octave_idx_type EigsComplexNonSymmetricMatrix (const ComplexMatrix& m, const std::string typ, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, ComplexMatrix &eig_vec, - ComplexColumnVector &eig_val, - const ComplexMatrix& b, ColumnVector &permB, - ComplexColumnVector &resid, - std::ostream &os, double tol = DBL_EPSILON, - int rvec = 0, bool cholB = 0, int disp = 0, - int maxit = 300); + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, ComplexMatrix &eig_vec, + ComplexColumnVector &eig_val, + const ComplexMatrix& b, ColumnVector &permB, + ComplexColumnVector &resid, + std::ostream &os, double tol = DBL_EPSILON, + int rvec = 0, bool cholB = 0, int disp = 0, + int maxit = 300); extern octave_idx_type EigsComplexNonSymmetricMatrix (const SparseComplexMatrix& m, - const std::string typ, octave_idx_type k, - octave_idx_type p, octave_idx_type &info, - ComplexMatrix &eig_vec, - ComplexColumnVector &eig_val, - const SparseComplexMatrix& b, - ColumnVector &permB, - ComplexColumnVector &resid, - std::ostream &os, double tol = DBL_EPSILON, - int rvec = 0, bool cholB = 0, int disp = 0, - int maxit = 300); + const std::string typ, octave_idx_type k, + octave_idx_type p, octave_idx_type &info, + ComplexMatrix &eig_vec, + ComplexColumnVector &eig_val, + const SparseComplexMatrix& b, + ColumnVector &permB, + ComplexColumnVector &resid, + std::ostream &os, double tol = DBL_EPSILON, + int rvec = 0, bool cholB = 0, int disp = 0, + int maxit = 300); extern octave_idx_type EigsComplexNonSymmetricMatrixShift (const ComplexMatrix& m, Complex sigma, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, - ComplexMatrix &eig_vec, - ComplexColumnVector &eig_val, - const ComplexMatrix& b, - ColumnVector &permB, - ComplexColumnVector &resid, - std::ostream &os, double tol = DBL_EPSILON, - int rvec = 0, bool cholB = 0, - int disp = 0, int maxit = 300); + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, + ComplexMatrix &eig_vec, + ComplexColumnVector &eig_val, + const ComplexMatrix& b, + ColumnVector &permB, + ComplexColumnVector &resid, + std::ostream &os, double tol = DBL_EPSILON, + int rvec = 0, bool cholB = 0, + int disp = 0, int maxit = 300); extern octave_idx_type EigsComplexNonSymmetricMatrixShift (const SparseComplexMatrix& m, - Complex sigma, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, - ComplexMatrix &eig_vec, - ComplexColumnVector &eig_val, - const SparseComplexMatrix& b, - ColumnVector &permB, - ComplexColumnVector &resid, - std::ostream &os, double tol = DBL_EPSILON, - int rvec = 0, bool cholB = 0, - int disp = 0, int maxit = 300); + Complex sigma, + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, + ComplexMatrix &eig_vec, + ComplexColumnVector &eig_val, + const SparseComplexMatrix& b, + ColumnVector &permB, + ComplexColumnVector &resid, + std::ostream &os, double tol = DBL_EPSILON, + int rvec = 0, bool cholB = 0, + int disp = 0, int maxit = 300); extern octave_idx_type EigsComplexNonSymmetricFunc (EigsComplexFunc fun, octave_idx_type n, - const std::string &_typ, Complex sigma, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, ComplexMatrix &eig_vec, - ComplexColumnVector &eig_val, - ComplexColumnVector &resid, std::ostream& os, - double tol = DBL_EPSILON, int rvec = 0, - bool cholB = 0, int disp = 0, int maxit = 300); + const std::string &_typ, Complex sigma, + octave_idx_type k, octave_idx_type p, + octave_idx_type &info, ComplexMatrix &eig_vec, + ComplexColumnVector &eig_val, + ComplexColumnVector &resid, std::ostream& os, + double tol = DBL_EPSILON, int rvec = 0, + bool cholB = 0, int disp = 0, int maxit = 300); #endif #ifndef _MSC_VER @@ -3796,7 +3796,7 @@ template static octave_idx_type lusolve (const SparseComplexMatrix&, const SparseComplexMatrix&, - ComplexMatrix&); + ComplexMatrix&); template static octave_idx_type lusolve (const Matrix&, const Matrix&, Matrix&); @@ -3806,7 +3806,7 @@ template static ComplexMatrix ltsolve (const SparseComplexMatrix&, const ColumnVector&, - const ComplexMatrix&); + const ComplexMatrix&); template static Matrix ltsolve (const SparseMatrix&, const ColumnVector&, const Matrix&); @@ -3819,7 +3819,7 @@ template static ComplexMatrix utsolve (const SparseComplexMatrix&, const ColumnVector&, - const ComplexMatrix&); + const ComplexMatrix&); template static Matrix utsolve (const SparseMatrix&, const ColumnVector&, const Matrix&);