Mercurial > octave
changeset 22953:fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
* EIG.cc, fEIG.cc: Use F77_INT instead of octave_idx_type for integer
data passed to Fortran subroutines.
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Mon, 26 Dec 2016 21:44:04 -0500 |
| parents | 7c9492d3b421 |
| children | 6cd3e9acf443 |
| files | liboctave/numeric/EIG.cc liboctave/numeric/fEIG.cc |
| diffstat | 2 files changed, 179 insertions(+), 147 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/numeric/EIG.cc Mon Dec 26 21:28:03 2016 -0500 +++ b/liboctave/numeric/EIG.cc Mon Dec 26 21:44:04 2016 -0500 @@ -40,12 +40,13 @@ if (a.is_symmetric ()) return symmetric_init (a, calc_rev, calc_lev); - octave_idx_type n = a.rows (); + F77_INT n = to_f77_int (a.rows ()); + F77_INT a_nc = to_f77_int (a.cols ()); - if (n != a.cols ()) + if (n != a_nc) (*current_liboctave_error_handler) ("EIG requires square matrix"); - octave_idx_type info = 0; + F77_INT info = 0; Matrix atmp = a; double *tmp_data = atmp.fortran_vec (); @@ -56,19 +57,19 @@ Array<double> wi (dim_vector (n, 1)); double *pwi = wi.fortran_vec (); - octave_idx_type tnvr = calc_rev ? n : 0; + F77_INT tnvr = calc_rev ? n : 0; Matrix vr (tnvr, tnvr); double *pvr = vr.fortran_vec (); - octave_idx_type tnvl = calc_lev ? n : 0; + F77_INT tnvl = calc_lev ? n : 0; Matrix vl (tnvl, tnvl); double *pvl = vl.fortran_vec (); - octave_idx_type lwork = -1; + F77_INT lwork = -1; double dummy_work; - octave_idx_type ilo; - octave_idx_type ihi; + F77_INT ilo; + F77_INT ihi; Array<double> scale (dim_vector (n, 1)); double *pscale = scale.fortran_vec (); @@ -81,7 +82,7 @@ Array<double> rcondv (dim_vector (n, 1)); double *prcondv = rcondv.fortran_vec (); - octave_idx_type dummy_iwork; + F77_INT dummy_iwork; F77_XFCN (dgeevx, DGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1), F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1), @@ -99,7 +100,7 @@ if (info != 0) (*current_liboctave_error_handler) ("dgeevx workspace query failed"); - lwork = static_cast<octave_idx_type> (dummy_work); + lwork = static_cast<F77_INT> (dummy_work); Array<double> work (dim_vector (lwork, 1)); double *pwork = work.fortran_vec (); @@ -123,20 +124,20 @@ (*current_liboctave_error_handler) ("dgeevx failed to converge"); lambda.resize (n); - octave_idx_type nvr = calc_rev ? n : 0; + F77_INT nvr = calc_rev ? n : 0; v.resize (nvr, nvr); - octave_idx_type nvl = calc_lev ? n : 0; + F77_INT nvl = calc_lev ? n : 0; w.resize (nvl, nvl); - for (octave_idx_type j = 0; j < n; j++) + for (F77_INT j = 0; j < n; j++) { if (wi.elem (j) == 0.0) { lambda.elem (j) = Complex (wr.elem (j)); - for (octave_idx_type i = 0; i < nvr; i++) + for (F77_INT i = 0; i < nvr; i++) v.elem (i, j) = vr.elem (i, j); - for (octave_idx_type i = 0; i < nvl; i++) + for (F77_INT i = 0; i < nvl; i++) w.elem (i, j) = vl.elem (i, j); } else @@ -147,7 +148,7 @@ lambda.elem (j) = Complex (wr.elem (j), wi.elem (j)); lambda.elem (j+1) = Complex (wr.elem (j+1), wi.elem (j+1)); - for (octave_idx_type i = 0; i < nvr; i++) + for (F77_INT i = 0; i < nvr; i++) { double real_part = vr.elem (i, j); double imag_part = vr.elem (i, j+1); @@ -155,7 +156,7 @@ v.elem (i, j+1) = Complex (real_part, -imag_part); } - for (octave_idx_type i = 0; i < nvl; i++) + for (F77_INT i = 0; i < nvl; i++) { double real_part = vl.elem (i, j); double imag_part = vl.elem (i, j+1); @@ -172,12 +173,13 @@ octave_idx_type EIG::symmetric_init (const Matrix& a, bool calc_rev, bool calc_lev) { - octave_idx_type n = a.rows (); + F77_INT n = to_f77_int (a.rows ()); + F77_INT a_nc = to_f77_int (a.cols ()); - if (n != a.cols ()) + if (n != a_nc) (*current_liboctave_error_handler) ("EIG requires square matrix"); - octave_idx_type info = 0; + F77_INT info = 0; Matrix atmp = a; double *tmp_data = atmp.fortran_vec (); @@ -185,7 +187,7 @@ ColumnVector wr (n); double *pwr = wr.fortran_vec (); - octave_idx_type lwork = -1; + F77_INT lwork = -1; double dummy_work; F77_XFCN (dsyev, DSYEV, (F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), @@ -197,7 +199,7 @@ if (info != 0) (*current_liboctave_error_handler) ("dsyev workspace query failed"); - lwork = static_cast<octave_idx_type> (dummy_work); + lwork = static_cast<F77_INT> (dummy_work); Array<double> work (dim_vector (lwork, 1)); double *pwork = work.fortran_vec (); @@ -230,12 +232,13 @@ if (a.is_hermitian ()) return hermitian_init (a, calc_rev, calc_lev); - octave_idx_type n = a.rows (); + F77_INT n = to_f77_int (a.rows ()); + F77_INT a_nc = to_f77_int (a.cols ()); - if (n != a.cols ()) + if (n != a_nc) (*current_liboctave_error_handler) ("EIG requires square matrix"); - octave_idx_type info = 0; + F77_INT info = 0; ComplexMatrix atmp = a; Complex *tmp_data = atmp.fortran_vec (); @@ -243,23 +246,23 @@ ComplexColumnVector wr (n); Complex *pw = wr.fortran_vec (); - octave_idx_type nvr = calc_rev ? n : 0; + F77_INT nvr = calc_rev ? n : 0; ComplexMatrix vrtmp (nvr, nvr); Complex *pvr = vrtmp.fortran_vec (); - octave_idx_type nvl = calc_lev ? n : 0; + F77_INT nvl = calc_lev ? n : 0; ComplexMatrix vltmp (nvl, nvl); Complex *pvl = vltmp.fortran_vec (); - octave_idx_type lwork = -1; + F77_INT lwork = -1; Complex dummy_work; - octave_idx_type lrwork = 2*n; + F77_INT lrwork = 2*n; Array<double> rwork (dim_vector (lrwork, 1)); double *prwork = rwork.fortran_vec (); - octave_idx_type ilo; - octave_idx_type ihi; + F77_INT ilo; + F77_INT ihi; Array<double> scale (dim_vector (n, 1)); double *pscale = scale.fortran_vec (); @@ -289,7 +292,7 @@ if (info != 0) (*current_liboctave_error_handler) ("zgeevx workspace query failed"); - lwork = static_cast<octave_idx_type> (dummy_work.real ()); + lwork = static_cast<F77_INT> (dummy_work.real ()); Array<Complex> work (dim_vector (lwork, 1)); Complex *pwork = work.fortran_vec (); @@ -323,12 +326,13 @@ octave_idx_type EIG::hermitian_init (const ComplexMatrix& a, bool calc_rev, bool calc_lev) { - octave_idx_type n = a.rows (); + F77_INT n = to_f77_int (a.rows ()); + F77_INT a_nc = to_f77_int (a.cols ()); - if (n != a.cols ()) + if (n != a_nc) (*current_liboctave_error_handler) ("EIG requires square matrix"); - octave_idx_type info = 0; + F77_INT info = 0; ComplexMatrix atmp = a; Complex *tmp_data = atmp.fortran_vec (); @@ -336,10 +340,10 @@ ColumnVector wr (n); double *pwr = wr.fortran_vec (); - octave_idx_type lwork = -1; + F77_INT lwork = -1; Complex dummy_work; - octave_idx_type lrwork = 3*n; + F77_INT lrwork = 3*n; Array<double> rwork (dim_vector (lrwork, 1)); double *prwork = rwork.fortran_vec (); @@ -354,7 +358,7 @@ if (info != 0) (*current_liboctave_error_handler) ("zheev workspace query failed"); - lwork = static_cast<octave_idx_type> (dummy_work.real ()); + lwork = static_cast<F77_INT> (dummy_work.real ()); Array<Complex> work (dim_vector (lwork, 1)); Complex *pwork = work.fortran_vec (); @@ -386,16 +390,19 @@ (*current_liboctave_error_handler) ("EIG: matrix contains Inf or NaN values"); - octave_idx_type n = a.rows (); - octave_idx_type nb = b.rows (); + F77_INT n = to_f77_int (a.rows ()); + F77_INT nb = to_f77_int (b.rows ()); - if (n != a.cols () || nb != b.cols ()) + F77_INT a_nc = to_f77_int (a.cols ()); + F77_INT b_nc = to_f77_int (b.cols ()); + + if (n != a_nc || nb != b_nc) (*current_liboctave_error_handler) ("EIG requires square matrix"); if (n != nb) (*current_liboctave_error_handler) ("EIG requires same size matrices"); - octave_idx_type info = 0; + F77_INT info = 0; Matrix tmp = b; double *tmp_data = tmp.fortran_vec (); @@ -426,15 +433,15 @@ Array<double> beta (dim_vector (n, 1)); double *pbeta = beta.fortran_vec (); - octave_idx_type tnvr = calc_rev ? n : 0; + F77_INT tnvr = calc_rev ? n : 0; Matrix vr (tnvr, tnvr); double *pvr = vr.fortran_vec (); - octave_idx_type tnvl = calc_lev ? n : 0; + F77_INT tnvl = calc_lev ? n : 0; Matrix vl (tnvl, tnvl); double *pvl = vl.fortran_vec (); - octave_idx_type lwork = -1; + F77_INT lwork = -1; double dummy_work; F77_XFCN (dggev, DGGEV, (F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1), @@ -449,7 +456,7 @@ if (info != 0) (*current_liboctave_error_handler) ("dggev workspace query failed"); - lwork = static_cast<octave_idx_type> (dummy_work); + lwork = static_cast<F77_INT> (dummy_work); Array<double> work (dim_vector (lwork, 1)); double *pwork = work.fortran_vec (); @@ -469,20 +476,20 @@ (*current_liboctave_error_handler) ("dggev failed to converge"); lambda.resize (n); - octave_idx_type nvr = calc_rev ? n : 0; + F77_INT nvr = calc_rev ? n : 0; v.resize (nvr, nvr); - octave_idx_type nvl = calc_lev ? n : 0; + F77_INT nvl = calc_lev ? n : 0; w.resize (nvl, nvl); - for (octave_idx_type j = 0; j < n; j++) + for (F77_INT j = 0; j < n; j++) { if (ai.elem (j) == 0.0) { lambda.elem (j) = Complex (ar.elem (j) / beta.elem (j)); - for (octave_idx_type i = 0; i < nvr; i++) + for (F77_INT i = 0; i < nvr; i++) v.elem (i, j) = vr.elem (i, j); - for (octave_idx_type i = 0; i < nvl; i++) + for (F77_INT i = 0; i < nvl; i++) w.elem (i, j) = vl.elem (i, j); } else @@ -495,14 +502,14 @@ lambda.elem (j+1) = Complex (ar.elem (j+1) / beta.elem (j+1), ai.elem (j+1) / beta.elem (j+1)); - for (octave_idx_type i = 0; i < nvr; i++) + for (F77_INT i = 0; i < nvr; i++) { double real_part = vr.elem (i, j); double imag_part = vr.elem (i, j+1); v.elem (i, j) = Complex (real_part, imag_part); v.elem (i, j+1) = Complex (real_part, -imag_part); } - for (octave_idx_type i = 0; i < nvl; i++) + for (F77_INT i = 0; i < nvl; i++) { double real_part = vl.elem (i, j); double imag_part = vl.elem (i, j+1); @@ -520,16 +527,19 @@ EIG::symmetric_init (const Matrix& a, const Matrix& b, bool calc_rev, bool calc_lev) { - octave_idx_type n = a.rows (); - octave_idx_type nb = b.rows (); + F77_INT n = to_f77_int (a.rows ()); + F77_INT nb = to_f77_int (b.rows ()); - if (n != a.cols () || nb != b.cols ()) + F77_INT a_nc = to_f77_int (a.cols ()); + F77_INT b_nc = to_f77_int (b.cols ()); + + if (n != a_nc || nb != b_nc) (*current_liboctave_error_handler) ("EIG requires square matrix"); if (n != nb) (*current_liboctave_error_handler) ("EIG requires same size matrices"); - octave_idx_type info = 0; + F77_INT info = 0; Matrix atmp = a; double *atmp_data = atmp.fortran_vec (); @@ -540,7 +550,7 @@ ColumnVector wr (n); double *pwr = wr.fortran_vec (); - octave_idx_type lwork = -1; + F77_INT lwork = -1; double dummy_work; F77_XFCN (dsygv, DSYGV, (1, F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), @@ -554,7 +564,7 @@ if (info != 0) (*current_liboctave_error_handler) ("dsygv workspace query failed"); - lwork = static_cast<octave_idx_type> (dummy_work); + lwork = static_cast<F77_INT> (dummy_work); Array<double> work (dim_vector (lwork, 1)); double *pwork = work.fortran_vec (); @@ -587,16 +597,19 @@ (*current_liboctave_error_handler) ("EIG: matrix contains Inf or NaN values"); - octave_idx_type n = a.rows (); - octave_idx_type nb = b.rows (); + F77_INT n = to_f77_int (a.rows ()); + F77_INT nb = to_f77_int (b.rows ()); - if (n != a.cols () || nb != b.cols ()) + F77_INT a_nc = to_f77_int (a.cols ()); + F77_INT b_nc = to_f77_int (b.cols ()); + + if (n != a_nc || nb != b_nc) (*current_liboctave_error_handler) ("EIG requires square matrix"); if (n != nb) (*current_liboctave_error_handler) ("EIG requires same size matrices"); - octave_idx_type info = 0; + F77_INT info = 0; ComplexMatrix tmp = b; Complex*tmp_data = tmp.fortran_vec (); @@ -624,18 +637,18 @@ ComplexColumnVector beta (n); Complex *pbeta = beta.fortran_vec (); - octave_idx_type nvr = calc_rev ? n : 0; + F77_INT nvr = calc_rev ? n : 0; ComplexMatrix vrtmp (nvr, nvr); Complex *pvr = vrtmp.fortran_vec (); - octave_idx_type nvl = calc_lev ? n : 0; + F77_INT nvl = calc_lev ? n : 0; ComplexMatrix vltmp (nvl, nvl); Complex *pvl = vltmp.fortran_vec (); - octave_idx_type lwork = -1; + F77_INT lwork = -1; Complex dummy_work; - octave_idx_type lrwork = 8*n; + F77_INT lrwork = 8*n; Array<double> rwork (dim_vector (lrwork, 1)); double *prwork = rwork.fortran_vec (); @@ -654,7 +667,7 @@ if (info != 0) (*current_liboctave_error_handler) ("zggev workspace query failed"); - lwork = static_cast<octave_idx_type> (dummy_work.real ()); + lwork = static_cast<F77_INT> (dummy_work.real ()); Array<Complex> work (dim_vector (lwork, 1)); Complex *pwork = work.fortran_vec (); @@ -678,7 +691,7 @@ lambda.resize (n); - for (octave_idx_type j = 0; j < n; j++) + for (F77_INT j = 0; j < n; j++) lambda.elem (j) = alpha.elem (j) / beta.elem (j); v = vrtmp; @@ -691,16 +704,19 @@ EIG::hermitian_init (const ComplexMatrix& a, const ComplexMatrix& b, bool calc_rev, bool calc_lev) { - octave_idx_type n = a.rows (); - octave_idx_type nb = b.rows (); + F77_INT n = to_f77_int (a.rows ()); + F77_INT nb = to_f77_int (b.rows ()); - if (n != a.cols () || nb != b.cols ()) + F77_INT a_nc = to_f77_int (a.cols ()); + F77_INT b_nc = to_f77_int (b.cols ()); + + if (n != a_nc || nb != b_nc) (*current_liboctave_error_handler) ("EIG requires square matrix"); if (n != nb) (*current_liboctave_error_handler) ("EIG requires same size matrices"); - octave_idx_type info = 0; + F77_INT info = 0; ComplexMatrix atmp = a; Complex *atmp_data = atmp.fortran_vec (); @@ -711,10 +727,10 @@ ColumnVector wr (n); double *pwr = wr.fortran_vec (); - octave_idx_type lwork = -1; + F77_INT lwork = -1; Complex dummy_work; - octave_idx_type lrwork = 3*n; + F77_INT lrwork = 3*n; Array<double> rwork (dim_vector (lrwork, 1)); double *prwork = rwork.fortran_vec (); @@ -730,7 +746,7 @@ if (info != 0) (*current_liboctave_error_handler) ("zhegv workspace query failed"); - lwork = static_cast<octave_idx_type> (dummy_work.real ()); + lwork = static_cast<F77_INT> (dummy_work.real ()); Array<Complex> work (dim_vector (lwork, 1)); Complex *pwork = work.fortran_vec ();
--- a/liboctave/numeric/fEIG.cc Mon Dec 26 21:28:03 2016 -0500 +++ b/liboctave/numeric/fEIG.cc Mon Dec 26 21:44:04 2016 -0500 @@ -41,12 +41,13 @@ if (a.is_symmetric ()) return symmetric_init (a, calc_rev, calc_lev); - octave_idx_type n = a.rows (); + F77_INT n = to_f77_int (a.rows ()); + F77_INT a_nc = to_f77_int (a.cols ()); - if (n != a.cols ()) + if (n != a_nc) (*current_liboctave_error_handler) ("EIG requires square matrix"); - octave_idx_type info = 0; + F77_INT info = 0; FloatMatrix atmp = a; float *tmp_data = atmp.fortran_vec (); @@ -57,19 +58,19 @@ Array<float> wi (dim_vector (n, 1)); float *pwi = wi.fortran_vec (); - volatile octave_idx_type nvr = calc_rev ? n : 0; + volatile F77_INT nvr = calc_rev ? n : 0; FloatMatrix vr (nvr, nvr); float *pvr = vr.fortran_vec (); - volatile octave_idx_type nvl = calc_lev ? n : 0; + volatile F77_INT nvl = calc_lev ? n : 0; FloatMatrix vl (nvl, nvl); float *pvl = vl.fortran_vec (); - octave_idx_type lwork = -1; + F77_INT lwork = -1; float dummy_work; - octave_idx_type ilo; - octave_idx_type ihi; + F77_INT ilo; + F77_INT ihi; Array<float> scale (dim_vector (n, 1)); float *pscale = scale.fortran_vec (); @@ -82,7 +83,7 @@ Array<float> rcondv (dim_vector (n, 1)); float *prcondv = rcondv.fortran_vec (); - octave_idx_type dummy_iwork; + F77_INT dummy_iwork; F77_XFCN (sgeevx, SGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1), F77_CONST_CHAR_ARG2 ("N", 1), @@ -100,7 +101,7 @@ if (info != 0) (*current_liboctave_error_handler) ("sgeevx workspace query failed"); - lwork = static_cast<octave_idx_type> (dummy_work); + lwork = static_cast<F77_INT> (dummy_work); Array<float> work (dim_vector (lwork, 1)); float *pwork = work.fortran_vec (); @@ -127,7 +128,7 @@ v.resize (nvr, nvr); w.resize (nvl, nvl); - for (octave_idx_type j = 0; j < n; j++) + for (F77_INT j = 0; j < n; j++) { if (wi.elem (j) == 0.0) { @@ -135,7 +136,7 @@ for (octave_idx_type i = 0; i < nvr; i++) v.elem (i, j) = vr.elem (i, j); - for (octave_idx_type i = 0; i < nvl; i++) + for (F77_INT i = 0; i < nvl; i++) w.elem (i, j) = vl.elem (i, j); } else @@ -146,14 +147,14 @@ lambda.elem (j) = FloatComplex (wr.elem (j), wi.elem (j)); lambda.elem (j+1) = FloatComplex (wr.elem (j+1), wi.elem (j+1)); - for (octave_idx_type i = 0; i < nvr; i++) + for (F77_INT i = 0; i < nvr; i++) { float real_part = vr.elem (i, j); float imag_part = vr.elem (i, j+1); v.elem (i, j) = FloatComplex (real_part, imag_part); v.elem (i, j+1) = FloatComplex (real_part, -imag_part); } - for (octave_idx_type i = 0; i < nvl; i++) + for (F77_INT i = 0; i < nvl; i++) { float real_part = vl.elem (i, j); float imag_part = vl.elem (i, j+1); @@ -170,12 +171,13 @@ octave_idx_type FloatEIG::symmetric_init (const FloatMatrix& a, bool calc_rev, bool calc_lev) { - octave_idx_type n = a.rows (); + F77_INT n = to_f77_int (a.rows ()); + F77_INT a_nc = to_f77_int (a.cols ()); - if (n != a.cols ()) + if (n != a_nc) (*current_liboctave_error_handler) ("EIG requires square matrix"); - octave_idx_type info = 0; + F77_INT info = 0; FloatMatrix atmp = a; float *tmp_data = atmp.fortran_vec (); @@ -183,7 +185,7 @@ FloatColumnVector wr (n); float *pwr = wr.fortran_vec (); - octave_idx_type lwork = -1; + F77_INT lwork = -1; float dummy_work; F77_XFCN (ssyev, SSYEV, (F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), @@ -195,7 +197,7 @@ if (info != 0) (*current_liboctave_error_handler) ("ssyev workspace query failed"); - lwork = static_cast<octave_idx_type> (dummy_work); + lwork = static_cast<F77_INT> (dummy_work); Array<float> work (dim_vector (lwork, 1)); float *pwork = work.fortran_vec (); @@ -229,12 +231,13 @@ if (a.is_hermitian ()) return hermitian_init (a, calc_rev, calc_lev); - octave_idx_type n = a.rows (); + F77_INT n = to_f77_int (a.rows ()); + F77_INT a_nc = to_f77_int (a.cols ()); - if (n != a.cols ()) + if (n != a_nc) (*current_liboctave_error_handler) ("EIG requires square matrix"); - octave_idx_type info = 0; + F77_INT info = 0; FloatComplexMatrix atmp = a; FloatComplex *tmp_data = atmp.fortran_vec (); @@ -242,23 +245,23 @@ FloatComplexColumnVector wr (n); FloatComplex *pw = wr.fortran_vec (); - octave_idx_type nvr = calc_rev ? n : 0; + F77_INT nvr = calc_rev ? n : 0; FloatComplexMatrix vrtmp (nvr, nvr); FloatComplex *pvr = vrtmp.fortran_vec (); - octave_idx_type nvl = calc_lev ? n : 0; + F77_INT nvl = calc_lev ? n : 0; FloatComplexMatrix vltmp (nvl, nvl); FloatComplex *pvl = vltmp.fortran_vec (); - octave_idx_type lwork = -1; + F77_INT lwork = -1; FloatComplex dummy_work; - octave_idx_type lrwork = 2*n; + F77_INT lrwork = 2*n; Array<float> rwork (dim_vector (lrwork, 1)); float *prwork = rwork.fortran_vec (); - octave_idx_type ilo; - octave_idx_type ihi; + F77_INT ilo; + F77_INT ihi; Array<float> scale (dim_vector (n, 1)); float *pscale = scale.fortran_vec (); @@ -287,7 +290,7 @@ if (info != 0) (*current_liboctave_error_handler) ("cgeevx workspace query failed"); - lwork = static_cast<octave_idx_type> (dummy_work.real ()); + lwork = static_cast<F77_INT> (dummy_work.real ()); Array<FloatComplex> work (dim_vector (lwork, 1)); FloatComplex *pwork = work.fortran_vec (); @@ -321,12 +324,13 @@ FloatEIG::hermitian_init (const FloatComplexMatrix& a, bool calc_rev, bool calc_lev) { - octave_idx_type n = a.rows (); + F77_INT n = to_f77_int (a.rows ()); + F77_INT a_nc = to_f77_int (a.cols ()); - if (n != a.cols ()) + if (n != a_nc) (*current_liboctave_error_handler) ("EIG requires square matrix"); - octave_idx_type info = 0; + F77_INT info = 0; FloatComplexMatrix atmp = a; FloatComplex *tmp_data = atmp.fortran_vec (); @@ -334,10 +338,10 @@ FloatColumnVector wr (n); float *pwr = wr.fortran_vec (); - octave_idx_type lwork = -1; + F77_INT lwork = -1; FloatComplex dummy_work; - octave_idx_type lrwork = 3*n; + F77_INT lrwork = 3*n; Array<float> rwork (dim_vector (lrwork, 1)); float *prwork = rwork.fortran_vec (); @@ -352,7 +356,7 @@ if (info != 0) (*current_liboctave_error_handler) ("cheev workspace query failed"); - lwork = static_cast<octave_idx_type> (dummy_work.real ()); + lwork = static_cast<F77_INT> (dummy_work.real ()); Array<FloatComplex> work (dim_vector (lwork, 1)); FloatComplex *pwork = work.fortran_vec (); @@ -384,16 +388,19 @@ (*current_liboctave_error_handler) ("EIG: matrix contains Inf or NaN values"); - octave_idx_type n = a.rows (); - octave_idx_type nb = b.rows (); + F77_INT n = to_f77_int (a.rows ()); + F77_INT nb = to_f77_int (b.rows ()); - if (n != a.cols () || nb != b.cols ()) + F77_INT a_nc = to_f77_int (a.cols ()); + F77_INT b_nc = to_f77_int (b.cols ()); + + if (n != a_nc || nb != b_nc) (*current_liboctave_error_handler) ("EIG requires square matrix"); if (n != nb) (*current_liboctave_error_handler) ("EIG requires same size matrices"); - octave_idx_type info = 0; + F77_INT info = 0; FloatMatrix tmp = b; float *tmp_data = tmp.fortran_vec (); @@ -423,15 +430,15 @@ Array<float> beta (dim_vector (n, 1)); float *pbeta = beta.fortran_vec (); - volatile octave_idx_type nvr = calc_rev ? n : 0; + volatile F77_INT nvr = calc_rev ? n : 0; FloatMatrix vr (nvr, nvr); float *pvr = vr.fortran_vec (); - volatile octave_idx_type nvl = calc_lev ? n : 0; + volatile F77_INT nvl = calc_lev ? n : 0; FloatMatrix vl (nvl, nvl); float *pvl = vl.fortran_vec (); - octave_idx_type lwork = -1; + F77_INT lwork = -1; float dummy_work; F77_XFCN (sggev, SGGEV, (F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1), @@ -446,7 +453,7 @@ if (info != 0) (*current_liboctave_error_handler) ("sggev workspace query failed"); - lwork = static_cast<octave_idx_type> (dummy_work); + lwork = static_cast<F77_INT> (dummy_work); Array<float> work (dim_vector (lwork, 1)); float *pwork = work.fortran_vec (); @@ -470,15 +477,15 @@ w.resize (nvl, nvl); - for (octave_idx_type j = 0; j < n; j++) + for (F77_INT j = 0; j < n; j++) { if (ai.elem (j) == 0.0) { lambda.elem (j) = FloatComplex (ar.elem (j) / beta.elem (j)); - for (octave_idx_type i = 0; i < nvr; i++) + for (F77_INT i = 0; i < nvr; i++) v.elem (i, j) = vr.elem (i, j); - for (octave_idx_type i = 0; i < nvl; i++) + for (F77_INT i = 0; i < nvl; i++) w.elem (i, j) = vl.elem (i, j); } else @@ -491,14 +498,14 @@ lambda.elem (j+1) = FloatComplex (ar.elem (j+1) / beta.elem (j+1), ai.elem (j+1) / beta.elem (j+1)); - for (octave_idx_type i = 0; i < nvr; i++) + for (F77_INT i = 0; i < nvr; i++) { float real_part = vr.elem (i, j); float imag_part = vr.elem (i, j+1); v.elem (i, j) = FloatComplex (real_part, imag_part); v.elem (i, j+1) = FloatComplex (real_part, -imag_part); } - for (octave_idx_type i = 0; i < nvl; i++) + for (F77_INT i = 0; i < nvl; i++) { float real_part = vl.elem (i, j); float imag_part = vl.elem (i, j+1); @@ -516,16 +523,19 @@ FloatEIG::symmetric_init (const FloatMatrix& a, const FloatMatrix& b, bool calc_rev, bool calc_lev) { - octave_idx_type n = a.rows (); - octave_idx_type nb = b.rows (); + F77_INT n = to_f77_int (a.rows ()); + F77_INT nb = to_f77_int (b.rows ()); - if (n != a.cols () || nb != b.cols ()) + F77_INT a_nc = to_f77_int (a.cols ()); + F77_INT b_nc = to_f77_int (b.cols ()); + + if (n != a_nc || nb != b_nc) (*current_liboctave_error_handler) ("EIG requires square matrix"); if (n != nb) (*current_liboctave_error_handler) ("EIG requires same size matrices"); - octave_idx_type info = 0; + F77_INT info = 0; FloatMatrix atmp = a; float *atmp_data = atmp.fortran_vec (); @@ -536,7 +546,7 @@ FloatColumnVector wr (n); float *pwr = wr.fortran_vec (); - octave_idx_type lwork = -1; + F77_INT lwork = -1; float dummy_work; F77_XFCN (ssygv, SSYGV, (1, F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), @@ -550,7 +560,7 @@ if (info != 0) (*current_liboctave_error_handler) ("ssygv workspace query failed"); - lwork = static_cast<octave_idx_type> (dummy_work); + lwork = static_cast<F77_INT> (dummy_work); Array<float> work (dim_vector (lwork, 1)); float *pwork = work.fortran_vec (); @@ -583,16 +593,19 @@ (*current_liboctave_error_handler) ("EIG: matrix contains Inf or NaN values"); - octave_idx_type n = a.rows (); - octave_idx_type nb = b.rows (); + F77_INT n = to_f77_int (a.rows ()); + F77_INT nb = to_f77_int (b.rows ()); - if (n != a.cols () || nb != b.cols ()) + F77_INT a_nc = to_f77_int (a.cols ()); + F77_INT b_nc = to_f77_int (b.cols ()); + + if (n != a_nc || nb != b_nc) (*current_liboctave_error_handler) ("EIG requires square matrix"); if (n != nb) (*current_liboctave_error_handler) ("EIG requires same size matrices"); - octave_idx_type info = 0; + F77_INT info = 0; FloatComplexMatrix tmp = b; FloatComplex *tmp_data = tmp.fortran_vec (); @@ -620,18 +633,18 @@ FloatComplexColumnVector beta (n); FloatComplex *pbeta = beta.fortran_vec (); - octave_idx_type nvr = calc_rev ? n : 0; + F77_INT nvr = calc_rev ? n : 0; FloatComplexMatrix vrtmp (nvr, nvr); FloatComplex *pvr = vrtmp.fortran_vec (); - octave_idx_type nvl = calc_lev ? n : 0; + F77_INT nvl = calc_lev ? n : 0; FloatComplexMatrix vltmp (nvl, nvl); FloatComplex *pvl = vltmp.fortran_vec (); - octave_idx_type lwork = -1; + F77_INT lwork = -1; FloatComplex dummy_work; - octave_idx_type lrwork = 8*n; + F77_INT lrwork = 8*n; Array<float> rwork (dim_vector (lrwork, 1)); float *prwork = rwork.fortran_vec (); @@ -648,7 +661,7 @@ if (info != 0) (*current_liboctave_error_handler) ("cggev workspace query failed"); - lwork = static_cast<octave_idx_type> (dummy_work.real ()); + lwork = static_cast<F77_INT> (dummy_work.real ()); Array<FloatComplex> work (dim_vector (lwork, 1)); FloatComplex *pwork = work.fortran_vec (); @@ -670,7 +683,7 @@ lambda.resize (n); - for (octave_idx_type j = 0; j < n; j++) + for (F77_INT j = 0; j < n; j++) lambda.elem (j) = alpha.elem (j) / beta.elem (j); v = vrtmp; @@ -684,16 +697,19 @@ const FloatComplexMatrix& b, bool calc_rev, bool calc_lev) { - octave_idx_type n = a.rows (); - octave_idx_type nb = b.rows (); + F77_INT n = to_f77_int (a.rows ()); + F77_INT nb = to_f77_int (b.rows ()); - if (n != a.cols () || nb != b.cols ()) + F77_INT a_nc = to_f77_int (a.cols ()); + F77_INT b_nc = to_f77_int (b.cols ()); + + if (n != a_nc || nb != b_nc) (*current_liboctave_error_handler) ("EIG requires square matrix"); if (n != nb) (*current_liboctave_error_handler) ("EIG requires same size matrices"); - octave_idx_type info = 0; + F77_INT info = 0; FloatComplexMatrix atmp = a; FloatComplex *atmp_data = atmp.fortran_vec (); @@ -704,10 +720,10 @@ FloatColumnVector wr (n); float *pwr = wr.fortran_vec (); - octave_idx_type lwork = -1; + F77_INT lwork = -1; FloatComplex dummy_work; - octave_idx_type lrwork = 3*n; + F77_INT lrwork = 3*n; Array<float> rwork (dim_vector (lrwork, 1)); float *prwork = rwork.fortran_vec (); @@ -723,7 +739,7 @@ if (info != 0) (*current_liboctave_error_handler) ("zhegv workspace query failed"); - lwork = static_cast<octave_idx_type> (dummy_work.real ()); + lwork = static_cast<F77_INT> (dummy_work.real ()); Array<FloatComplex> work (dim_vector (lwork, 1)); FloatComplex *pwork = work.fortran_vec ();