/* * eig.c * * Code generation for function 'eig' * * C source code generated on: Wed Aug 26 14:59:33 2015 * */ /* Include files */ #include "rt_nonfinite.h" #include "Select_Ground_Motions.h" #include "eig.h" #include "inv.h" #include "Select_Ground_Motions_emxutil.h" #include "SP_1996.h" #include "sqrt.h" #include "Select_Ground_Motions_mexutil.h" #include "Select_Ground_Motions_data.h" /* Variable Definitions */ static emlrtRSInfo oi_emlrtRSI = { 40, "eig", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/matfun/eig.m" }; static emlrtRSInfo pi_emlrtRSI = { 41, "eig", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/matfun/eig.m" }; static emlrtRSInfo qi_emlrtRSI = { 55, "eig", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/matfun/eig.m" }; static emlrtRSInfo ri_emlrtRSI = { 59, "eig", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/matfun/eig.m" }; static emlrtRSInfo si_emlrtRSI = { 9, "eml_xgeev", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/eml_xgeev.m" }; static emlrtRSInfo ti_emlrtRSI = { 8, "eml_lapack_xgeev", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/internal/eml_lapack_xgeev.m" }; static emlrtRSInfo ui_emlrtRSI = { 15, "eml_lapack_xgeev", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/internal/eml_lapack_xgeev.m" }; static emlrtRSInfo vi_emlrtRSI = { 16, "eml_lapack_xgeev", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/internal/eml_lapack_xgeev.m" }; static emlrtRSInfo wi_emlrtRSI = { 17, "eml_lapack_xgeev", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/internal/eml_lapack_xgeev.m" }; static emlrtRSInfo xi_emlrtRSI = { 84, "eml_matlab_zggev", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zggev.m" }; static emlrtRSInfo yi_emlrtRSI = { 100, "eml_matlab_zggev", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zggev.m" }; static emlrtRSInfo aj_emlrtRSI = { 101, "eml_matlab_zggev", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zggev.m" }; static emlrtRSInfo bj_emlrtRSI = { 106, "eml_matlab_zggev", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zggev.m" }; static emlrtRSInfo cj_emlrtRSI = { 108, "eml_matlab_zggev", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zggev.m" }; static emlrtRSInfo dj_emlrtRSI = { 23, "eml_matlab_zggbal", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zggbal.m" }; static emlrtRSInfo ej_emlrtRSI = { 27, "eml_matlab_zggbal", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zggbal.m" }; static emlrtRSInfo fj_emlrtRSI = { 36, "eml_matlab_zggbal", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zggbal.m" }; static emlrtRSInfo gj_emlrtRSI = { 40, "eml_matlab_zggbal", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zggbal.m" }; static emlrtRSInfo hj_emlrtRSI = { 140, "eml_matlab_zggbal", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zggbal.m" }; static emlrtRSInfo ij_emlrtRSI = { 60, "eml_matlab_zggbal", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zggbal.m" }; static emlrtRSInfo jj_emlrtRSI = { 77, "eml_matlab_zggbal", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zggbal.m" }; static emlrtRSInfo kj_emlrtRSI = { 102, "eml_matlab_zggbal", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zggbal.m" }; static emlrtRSInfo lj_emlrtRSI = { 106, "eml_matlab_zggbal", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zggbal.m" }; static emlrtRSInfo mj_emlrtRSI = { 30, "eml_matlab_zgghrd", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zgghrd.m" }; static emlrtRSInfo nj_emlrtRSI = { 62, "eml_matlab_zgghrd", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zgghrd.m" }; static emlrtRSInfo oj_emlrtRSI = { 64, "eml_matlab_zgghrd", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zgghrd.m" }; static emlrtRSInfo pj_emlrtRSI = { 77, "eml_matlab_zgghrd", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zgghrd.m" }; static emlrtRSInfo qj_emlrtRSI = { 79, "eml_matlab_zgghrd", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zgghrd.m" }; static emlrtRSInfo rj_emlrtRSI = { 25, "eye", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/elmat/eye.m" }; static emlrtRSInfo sj_emlrtRSI = { 50, "eye", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/elmat/eye.m" }; static emlrtRSInfo tj_emlrtRSI = { 59, "eye", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/elmat/eye.m" }; static emlrtRSInfo vj_emlrtRSI = { 67, "eml_matlab_zlartg", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zlartg.m" }; static emlrtRSInfo wj_emlrtRSI = { 93, "eml_matlab_zlartg", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zlartg.m" }; static emlrtRSInfo xj_emlrtRSI = { 102, "eml_matlab_zlartg", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zlartg.m" }; static emlrtRSInfo yj_emlrtRSI = { 106, "eml_matlab_zlartg", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zlartg.m" }; static emlrtRSInfo ak_emlrtRSI = { 10, "eml_zrot_rows", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_zrot_rows.m" }; static emlrtRSInfo bk_emlrtRSI = { 10, "eml_zrot_cols", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_zrot_cols.m" }; static emlrtRSInfo ck_emlrtRSI = { 444, "eml_matlab_zhgeqz", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zhgeqz.m" }; static emlrtRSInfo dk_emlrtRSI = { 432, "eml_matlab_zhgeqz", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zhgeqz.m" }; static emlrtRSInfo ek_emlrtRSI = { 420, "eml_matlab_zhgeqz", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zhgeqz.m" }; static emlrtRSInfo fk_emlrtRSI = { 418, "eml_matlab_zhgeqz", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zhgeqz.m" }; static emlrtRSInfo gk_emlrtRSI = { 409, "eml_matlab_zhgeqz", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zhgeqz.m" }; static emlrtRSInfo hk_emlrtRSI = { 406, "eml_matlab_zhgeqz", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zhgeqz.m" }; static emlrtRSInfo ik_emlrtRSI = { 399, "eml_matlab_zhgeqz", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zhgeqz.m" }; static emlrtRSInfo jk_emlrtRSI = { 107, "eml_matlab_zhgeqz", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zhgeqz.m" }; static emlrtRSInfo kk_emlrtRSI = { 47, "eml_matlab_zhgeqz", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zhgeqz.m" }; static emlrtRSInfo lk_emlrtRSI = { 34, "eml_matlab_zhgeqz", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zhgeqz.m" }; static emlrtRSInfo mk_emlrtRSI = { 20, "eml_matlab_zlanhs", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zlanhs.m" }; static emlrtRSInfo nk_emlrtRSI = { 21, "eml_matlab_zlanhs", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zlanhs.m" }; static emlrtRSInfo ok_emlrtRSI = { 57, "eml_matlab_zlanhs", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zlanhs.m" }; static emlrtRSInfo pk_emlrtRSI = { 22, "eml_matlab_zggbak", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zggbak.m" }; static emlrtRSInfo qk_emlrtRSI = { 32, "eml_matlab_zggbak", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zggbak.m" }; static emlrtRSInfo rk_emlrtRSI = { 36, "eml_matlab_zggbak", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zggbak.m" }; static emlrtRTEInfo wb_emlrtRTEI = { 1, 18, "eig", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/matfun/eig.m" }; static emlrtRTEInfo xb_emlrtRTEI = { 1, 34, "eml_matlab_zggev", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zggev.m" }; static emlrtRTEInfo yb_emlrtRTEI = { 1, 33, "eml_matlab_zggbal", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zggbal.m" }; static emlrtRTEInfo ac_emlrtRTEI = { 1, 38, "eml_matlab_zhgeqz", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_zhgeqz.m" }; static emlrtRTEInfo bc_emlrtRTEI = { 1, 14, "eml_matlab_ztgevc", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_ztgevc.m" }; static emlrtRTEInfo wd_emlrtRTEI = { 13, 1, "eml_matlab_ztgevc", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_ztgevc.m" }; static emlrtRTEInfo xd_emlrtRTEI = { 14, 1, "eml_matlab_ztgevc", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_ztgevc.m" }; static emlrtRTEInfo yd_emlrtRTEI = { 23, 1, "eml_matlab_ztgevc", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_ztgevc.m" }; static emlrtRTEInfo he_emlrtRTEI = { 113, 5, "eml_matlab_ztgevc", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_ztgevc.m" }; static emlrtRTEInfo ie_emlrtRTEI = { 48, 1, "eml_matlab_ztgevc", "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/eml/lapack/matlab/eml_matlab_ztgevc.m" }; /* Function Declarations */ static void b_eml_matlab_zlartg(const creal_T f, const creal_T g, real_T *cs, creal_T *sn); static void c_eml_warning(void); static void d_eml_warning(void); static int32_T div_s32(int32_T numerator, int32_T denominator); static creal_T eml_div(const creal_T x, real_T y); static void eml_matlab_zggbal(emxArray_creal_T *A, int32_T *ilo, int32_T *ihi, emxArray_int32_T *rscale); static void eml_matlab_zggev(emxArray_creal_T *A, real_T *info, emxArray_creal_T *alpha1, emxArray_creal_T *beta1, emxArray_creal_T *V); static void eml_matlab_zhgeqz(emxArray_creal_T *A, int32_T ilo, int32_T ihi, emxArray_creal_T *Z, real_T *info, emxArray_creal_T *alpha1, emxArray_creal_T * beta1); static real_T eml_matlab_zlanhs(const emxArray_creal_T *A, int32_T ilo, int32_T ihi); static void eml_matlab_zlartg(const creal_T f, const creal_T g, real_T *cs, creal_T *sn, creal_T *r); static void eml_matlab_ztgevc(const emxArray_creal_T *A, emxArray_creal_T *V); /* Function Definitions */ static void b_eml_matlab_zlartg(const creal_T f, const creal_T g, real_T *cs, creal_T *sn) { real_T scale; real_T f2s; real_T g2; real_T fs_re; real_T fs_im; real_T gs_re; real_T gs_im; int32_T count; int32_T rescaledir; boolean_T guard1 = FALSE; real_T g2s; boolean_T b8; boolean_T b9; scale = muDoubleScalarAbs(f.re); f2s = muDoubleScalarAbs(f.im); if (f2s > scale) { scale = f2s; } f2s = muDoubleScalarAbs(g.re); g2 = muDoubleScalarAbs(g.im); if (g2 > f2s) { f2s = g2; } if (f2s > scale) { scale = f2s; } fs_re = f.re; fs_im = f.im; gs_re = g.re; gs_im = g.im; count = 0; rescaledir = 0; guard1 = FALSE; if (scale >= 7.4428285367870146E+137) { do { count++; fs_re *= 1.3435752215134178E-138; fs_im *= 1.3435752215134178E-138; gs_re *= 1.3435752215134178E-138; gs_im *= 1.3435752215134178E-138; scale *= 1.3435752215134178E-138; } while (!(scale < 7.4428285367870146E+137)); rescaledir = 1; guard1 = TRUE; } else if (scale <= 1.3435752215134178E-138) { if ((g.re == 0.0) && (g.im == 0.0)) { *cs = 1.0; sn->re = 0.0; sn->im = 0.0; } else { do { count++; fs_re *= 7.4428285367870146E+137; fs_im *= 7.4428285367870146E+137; gs_re *= 7.4428285367870146E+137; gs_im *= 7.4428285367870146E+137; scale *= 7.4428285367870146E+137; } while (!(scale > 1.3435752215134178E-138)); rescaledir = -1; guard1 = TRUE; } } else { guard1 = TRUE; } if (guard1 == TRUE) { scale = fs_re * fs_re + fs_im * fs_im; g2 = gs_re * gs_re + gs_im * gs_im; f2s = g2; if (1.0 > g2) { f2s = 1.0; } if (scale <= f2s * 2.0041683600089728E-292) { if ((f.re == 0.0) && (f.im == 0.0)) { *cs = 0.0; scale = muDoubleScalarHypot(gs_re, gs_im); sn->re = gs_re / scale; sn->im = -gs_im / scale; } else { emlrtPushRtStackR2012b(&vj_emlrtRSI, emlrtRootTLSGlobal); if (g2 < 0.0) { emlrtPushRtStackR2012b(&md_emlrtRSI, emlrtRootTLSGlobal); e_eml_error(); emlrtPopRtStackR2012b(&md_emlrtRSI, emlrtRootTLSGlobal); } g2s = muDoubleScalarSqrt(g2); emlrtPopRtStackR2012b(&vj_emlrtRSI, emlrtRootTLSGlobal); *cs = muDoubleScalarHypot(fs_re, fs_im) / g2s; f2s = muDoubleScalarAbs(f.re); g2 = muDoubleScalarAbs(f.im); if (g2 > f2s) { f2s = g2; } if (f2s > 1.0) { scale = muDoubleScalarHypot(f.re, f.im); fs_re = f.re / scale; fs_im = f.im / scale; } else { f2s = 7.4428285367870146E+137 * f.re; g2 = 7.4428285367870146E+137 * f.im; scale = muDoubleScalarHypot(f2s, g2); fs_re = f2s / scale; fs_im = g2 / scale; } gs_re /= g2s; gs_im = -gs_im / g2s; sn->re = fs_re * gs_re - fs_im * gs_im; sn->im = fs_re * gs_im + fs_im * gs_re; } } else { emlrtPushRtStackR2012b(&wj_emlrtRSI, emlrtRootTLSGlobal); f2s = 1.0 + g2 / scale; if (f2s < 0.0) { emlrtPushRtStackR2012b(&md_emlrtRSI, emlrtRootTLSGlobal); e_eml_error(); emlrtPopRtStackR2012b(&md_emlrtRSI, emlrtRootTLSGlobal); } f2s = muDoubleScalarSqrt(f2s); emlrtPopRtStackR2012b(&wj_emlrtRSI, emlrtRootTLSGlobal); *cs = 1.0 / f2s; scale += g2; fs_re = f2s * fs_re / scale; fs_im = f2s * fs_im / scale; sn->re = fs_re * gs_re - fs_im * -gs_im; sn->im = fs_re * -gs_im + fs_im * gs_re; if (rescaledir > 0) { emlrtPushRtStackR2012b(&xj_emlrtRSI, emlrtRootTLSGlobal); if (1 > count) { b8 = FALSE; } else { b8 = (count > 2147483646); } if (b8) { emlrtPushRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); check_forloop_overflow_error(); emlrtPopRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); } emlrtPopRtStackR2012b(&xj_emlrtRSI, emlrtRootTLSGlobal); } else { if (rescaledir < 0) { emlrtPushRtStackR2012b(&yj_emlrtRSI, emlrtRootTLSGlobal); if (1 > count) { b9 = FALSE; } else { b9 = (count > 2147483646); } if (b9) { emlrtPushRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); check_forloop_overflow_error(); emlrtPopRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); } emlrtPopRtStackR2012b(&yj_emlrtRSI, emlrtRootTLSGlobal); } } } } } static void c_eml_warning(void) { const mxArray *y; static const int32_T iv54[2] = { 1, 26 }; const mxArray *m14; char_T cv67[26]; int32_T i; static const char_T cv68[26] = { 'C', 'o', 'd', 'e', 'r', ':', 't', 'o', 'o', 'l', 'b', 'o', 'x', ':', 'e', 'i', 'g', '_', 'Q', 'Z', 'f', 'a', 'i', 'l', 'e', 'd' }; emlrtPushRtStackR2012b(&ee_emlrtRSI, emlrtRootTLSGlobal); y = NULL; m14 = mxCreateCharArray(2, iv54); for (i = 0; i < 26; i++) { cv67[i] = cv68[i]; } emlrtInitCharArrayR2013a(emlrtRootTLSGlobal, 26, m14, cv67); emlrtAssign(&y, m14); warning(message(y, &x_emlrtMCI), &y_emlrtMCI); emlrtPopRtStackR2012b(&ee_emlrtRSI, emlrtRootTLSGlobal); } static void d_eml_warning(void) { const mxArray *y; static const int32_T iv55[2] = { 1, 34 }; const mxArray *m15; char_T cv69[34]; int32_T i; static const char_T cv70[34] = { 'C', 'o', 'd', 'e', 'r', ':', 't', 'o', 'o', 'l', 'b', 'o', 'x', ':', 'e', 'i', 'g', '_', 'Q', 'Z', 'n', 'o', 'n', 'c', 'o', 'n', 'v', 'e', 'r', 'g', 'e', 'n', 'c', 'e' }; emlrtPushRtStackR2012b(&ee_emlrtRSI, emlrtRootTLSGlobal); y = NULL; m15 = mxCreateCharArray(2, iv55); for (i = 0; i < 34; i++) { cv69[i] = cv70[i]; } emlrtInitCharArrayR2013a(emlrtRootTLSGlobal, 34, m15, cv69); emlrtAssign(&y, m15); warning(message(y, &x_emlrtMCI), &y_emlrtMCI); emlrtPopRtStackR2012b(&ee_emlrtRSI, emlrtRootTLSGlobal); } static int32_T div_s32(int32_T numerator, int32_T denominator) { int32_T quotient; uint32_T absNumerator; uint32_T absDenominator; int32_T quotientNeedsNegation; if (denominator == 0) { if (numerator >= 0) { quotient = MAX_int32_T; } else { quotient = MIN_int32_T; } emlrtDivisionByZeroErrorR2012b(0, emlrtRootTLSGlobal); } else { if (numerator >= 0) { absNumerator = (uint32_T)numerator; } else { absNumerator = (uint32_T)-numerator; } if (denominator >= 0) { absDenominator = (uint32_T)denominator; } else { absDenominator = (uint32_T)-denominator; } quotientNeedsNegation = ((numerator < 0) != (denominator < 0)); absNumerator /= absDenominator; if ((uint32_T)quotientNeedsNegation) { quotient = -(int32_T)absNumerator; } else { quotient = (int32_T)absNumerator; } } return quotient; } static creal_T eml_div(const creal_T x, real_T y) { creal_T z; if (x.im == 0.0) { z.re = x.re / y; z.im = 0.0; } else if (x.re == 0.0) { z.re = 0.0; z.im = x.im / y; } else { z.re = x.re / y; z.im = x.im / y; } return z; } static void eml_matlab_zggbal(emxArray_creal_T *A, int32_T *ilo, int32_T *ihi, emxArray_int32_T *rscale) { int32_T ii; int32_T nzcount; emxArray_creal_T *b_A; int32_T exitg2; int32_T i; int32_T j; boolean_T found; boolean_T exitg5; boolean_T overflow; int32_T jj; boolean_T exitg6; boolean_T guard2 = FALSE; int32_T i29; real_T atmp_re; real_T atmp_im; int32_T exitg1; boolean_T exitg3; boolean_T exitg4; boolean_T guard1 = FALSE; emlrtHeapReferenceStackEnterFcnR2012b(emlrtRootTLSGlobal); ii = rscale->size[0]; rscale->size[0] = A->size[0]; emxEnsureCapacity((emxArray__common *)rscale, ii, (int32_T)sizeof(int32_T), &yb_emlrtRTEI); nzcount = A->size[0]; for (ii = 0; ii < nzcount; ii++) { rscale->data[ii] = 0; } *ilo = 1; *ihi = A->size[0]; if (A->size[0] <= 1) { *ihi = 1; rscale->data[0] = 1; } else { emxInit_creal_T(&b_A, 2, &yb_emlrtRTEI, TRUE); do { exitg2 = 0; emlrtPushRtStackR2012b(&dj_emlrtRSI, emlrtRootTLSGlobal); i = 0; j = 0; found = FALSE; ii = *ihi; exitg5 = FALSE; while ((exitg5 == FALSE) && (ii > 0)) { nzcount = 0; i = ii; j = *ihi; emlrtPushRtStackR2012b(&hj_emlrtRSI, emlrtRootTLSGlobal); if (1 > *ihi) { overflow = FALSE; } else { overflow = (*ihi > 2147483646); } if (overflow) { emlrtPushRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); check_forloop_overflow_error(); emlrtPopRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); } emlrtPopRtStackR2012b(&hj_emlrtRSI, emlrtRootTLSGlobal); jj = 1; exitg6 = FALSE; while ((exitg6 == FALSE) && (jj <= *ihi)) { overflow = ((A->data[(ii + A->size[0] * (jj - 1)) - 1].re != 0.0) || (A->data[(ii + A->size[0] * (jj - 1)) - 1].im != 0.0)); guard2 = FALSE; if (overflow || (ii == jj)) { if (nzcount == 0) { j = jj; nzcount = 1; guard2 = TRUE; } else { nzcount = 2; exitg6 = TRUE; } } else { guard2 = TRUE; } if (guard2 == TRUE) { jj++; } } if (nzcount < 2) { found = TRUE; exitg5 = TRUE; } else { ii--; } } emlrtPopRtStackR2012b(&dj_emlrtRSI, emlrtRootTLSGlobal); if (!found) { exitg2 = 2; } else { emlrtPushRtStackR2012b(&ej_emlrtRSI, emlrtRootTLSGlobal); ii = b_A->size[0] * b_A->size[1]; b_A->size[0] = A->size[0]; b_A->size[1] = A->size[1]; emxEnsureCapacity((emxArray__common *)b_A, ii, (int32_T)sizeof(creal_T), &yb_emlrtRTEI); nzcount = A->size[1]; for (ii = 0; ii < nzcount; ii++) { jj = A->size[0]; for (i29 = 0; i29 < jj; i29++) { b_A->data[i29 + b_A->size[0] * ii] = A->data[i29 + A->size[0] * ii]; } } if (i != *ihi) { emlrtPushRtStackR2012b(&ij_emlrtRSI, emlrtRootTLSGlobal); if (1 > A->size[0]) { overflow = FALSE; } else { overflow = (A->size[0] > 2147483646); } if (overflow) { emlrtPushRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); check_forloop_overflow_error(); emlrtPopRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); } emlrtPopRtStackR2012b(&ij_emlrtRSI, emlrtRootTLSGlobal); for (ii = 0; ii + 1 <= A->size[0]; ii++) { atmp_re = b_A->data[(i + b_A->size[0] * ii) - 1].re; atmp_im = b_A->data[(i + b_A->size[0] * ii) - 1].im; b_A->data[(i + b_A->size[0] * ii) - 1] = b_A->data[(*ihi + b_A-> size[0] * ii) - 1]; b_A->data[(*ihi + b_A->size[0] * ii) - 1].re = atmp_re; b_A->data[(*ihi + b_A->size[0] * ii) - 1].im = atmp_im; } } if (j != *ihi) { emlrtPushRtStackR2012b(&jj_emlrtRSI, emlrtRootTLSGlobal); if (1 > *ihi) { overflow = FALSE; } else { overflow = (*ihi > 2147483646); } if (overflow) { emlrtPushRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); check_forloop_overflow_error(); emlrtPopRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); } emlrtPopRtStackR2012b(&jj_emlrtRSI, emlrtRootTLSGlobal); for (ii = 0; ii + 1 <= *ihi; ii++) { atmp_re = b_A->data[ii + b_A->size[0] * (j - 1)].re; atmp_im = b_A->data[ii + b_A->size[0] * (j - 1)].im; b_A->data[ii + b_A->size[0] * (j - 1)] = b_A->data[ii + b_A->size[0] * (*ihi - 1)]; b_A->data[ii + b_A->size[0] * (*ihi - 1)].re = atmp_re; b_A->data[ii + b_A->size[0] * (*ihi - 1)].im = atmp_im; } } emlrtPopRtStackR2012b(&ej_emlrtRSI, emlrtRootTLSGlobal); ii = A->size[0] * A->size[1]; A->size[0] = b_A->size[0]; A->size[1] = b_A->size[1]; emxEnsureCapacity((emxArray__common *)A, ii, (int32_T)sizeof(creal_T), &yb_emlrtRTEI); nzcount = b_A->size[1]; for (ii = 0; ii < nzcount; ii++) { jj = b_A->size[0]; for (i29 = 0; i29 < jj; i29++) { A->data[i29 + A->size[0] * ii] = b_A->data[i29 + b_A->size[0] * ii]; } } rscale->data[*ihi - 1] = j; (*ihi)--; if (*ihi == 1) { rscale->data[0] = 1; exitg2 = 1; } } } while (exitg2 == 0); if (exitg2 == 1) { } else { do { exitg1 = 0; emlrtPushRtStackR2012b(&fj_emlrtRSI, emlrtRootTLSGlobal); i = 0; j = 0; found = FALSE; emlrtPushRtStackR2012b(&kj_emlrtRSI, emlrtRootTLSGlobal); if (*ilo > *ihi) { overflow = FALSE; } else { overflow = (*ihi > 2147483646); } if (overflow) { emlrtPushRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); check_forloop_overflow_error(); emlrtPopRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); } emlrtPopRtStackR2012b(&kj_emlrtRSI, emlrtRootTLSGlobal); jj = *ilo; exitg3 = FALSE; while ((exitg3 == FALSE) && (jj <= *ihi)) { nzcount = 0; i = *ihi; j = jj; emlrtPushRtStackR2012b(&lj_emlrtRSI, emlrtRootTLSGlobal); if (*ilo > *ihi) { overflow = FALSE; } else { overflow = (*ihi > 2147483646); } if (overflow) { emlrtPushRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); check_forloop_overflow_error(); emlrtPopRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); } emlrtPopRtStackR2012b(&lj_emlrtRSI, emlrtRootTLSGlobal); ii = *ilo; exitg4 = FALSE; while ((exitg4 == FALSE) && (ii <= *ihi)) { overflow = ((A->data[(ii + A->size[0] * (jj - 1)) - 1].re != 0.0) || (A->data[(ii + A->size[0] * (jj - 1)) - 1].im != 0.0)); guard1 = FALSE; if (overflow || (ii == jj)) { if (nzcount == 0) { i = ii; nzcount = 1; guard1 = TRUE; } else { nzcount = 2; exitg4 = TRUE; } } else { guard1 = TRUE; } if (guard1 == TRUE) { ii++; } } if (nzcount < 2) { found = TRUE; exitg3 = TRUE; } else { jj++; } } emlrtPopRtStackR2012b(&fj_emlrtRSI, emlrtRootTLSGlobal); if (!found) { exitg1 = 1; } else { emlrtPushRtStackR2012b(&gj_emlrtRSI, emlrtRootTLSGlobal); ii = b_A->size[0] * b_A->size[1]; b_A->size[0] = A->size[0]; b_A->size[1] = A->size[1]; emxEnsureCapacity((emxArray__common *)b_A, ii, (int32_T)sizeof(creal_T), &yb_emlrtRTEI); nzcount = A->size[1]; for (ii = 0; ii < nzcount; ii++) { jj = A->size[0]; for (i29 = 0; i29 < jj; i29++) { b_A->data[i29 + b_A->size[0] * ii] = A->data[i29 + A->size[0] * ii]; } } if (i != *ilo) { emlrtPushRtStackR2012b(&ij_emlrtRSI, emlrtRootTLSGlobal); if (*ilo > A->size[0]) { overflow = FALSE; } else { overflow = (A->size[0] > 2147483646); } if (overflow) { emlrtPushRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); check_forloop_overflow_error(); emlrtPopRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); } emlrtPopRtStackR2012b(&ij_emlrtRSI, emlrtRootTLSGlobal); for (ii = *ilo - 1; ii + 1 <= A->size[0]; ii++) { atmp_re = b_A->data[(i + b_A->size[0] * ii) - 1].re; atmp_im = b_A->data[(i + b_A->size[0] * ii) - 1].im; b_A->data[(i + b_A->size[0] * ii) - 1] = b_A->data[(*ilo + b_A->size[0] * ii) - 1]; b_A->data[(*ilo + b_A->size[0] * ii) - 1].re = atmp_re; b_A->data[(*ilo + b_A->size[0] * ii) - 1].im = atmp_im; } } if (j != *ilo) { emlrtPushRtStackR2012b(&jj_emlrtRSI, emlrtRootTLSGlobal); if (1 > *ihi) { overflow = FALSE; } else { overflow = (*ihi > 2147483646); } if (overflow) { emlrtPushRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); check_forloop_overflow_error(); emlrtPopRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); } emlrtPopRtStackR2012b(&jj_emlrtRSI, emlrtRootTLSGlobal); for (ii = 0; ii + 1 <= *ihi; ii++) { atmp_re = b_A->data[ii + b_A->size[0] * (j - 1)].re; atmp_im = b_A->data[ii + b_A->size[0] * (j - 1)].im; b_A->data[ii + b_A->size[0] * (j - 1)] = b_A->data[ii + b_A->size [0] * (*ilo - 1)]; b_A->data[ii + b_A->size[0] * (*ilo - 1)].re = atmp_re; b_A->data[ii + b_A->size[0] * (*ilo - 1)].im = atmp_im; } } emlrtPopRtStackR2012b(&gj_emlrtRSI, emlrtRootTLSGlobal); ii = A->size[0] * A->size[1]; A->size[0] = b_A->size[0]; A->size[1] = b_A->size[1]; emxEnsureCapacity((emxArray__common *)A, ii, (int32_T)sizeof(creal_T), &yb_emlrtRTEI); nzcount = b_A->size[1]; for (ii = 0; ii < nzcount; ii++) { jj = b_A->size[0]; for (i29 = 0; i29 < jj; i29++) { A->data[i29 + A->size[0] * ii] = b_A->data[i29 + b_A->size[0] * ii]; } } rscale->data[*ilo - 1] = j; (*ilo)++; if (*ilo == *ihi) { rscale->data[*ilo - 1] = *ilo; exitg1 = 1; } } } while (exitg1 == 0); } emxFree_creal_T(&b_A); } emlrtHeapReferenceStackLeaveFcnR2012b(emlrtRootTLSGlobal); } static void eml_matlab_zggev(emxArray_creal_T *A, real_T *info, emxArray_creal_T *alpha1, emxArray_creal_T *beta1, emxArray_creal_T *V) { int32_T n; int32_T m; int32_T q; int32_T jcol; real_T anrm; boolean_T exitg1; real_T absxk; boolean_T ilascl; real_T anrmto; real_T ctoc; boolean_T notdone; real_T cfrom1; real_T cto1; real_T mul; emxArray_creal_T *b_A; emxArray_int32_T *rscale; emxArray_creal_T *c_A; int32_T ihi; int32_T ilo; const mxArray *y; static const int32_T iv53[2] = { 1, 21 }; const mxArray *m13; char_T cv65[21]; int32_T i; static const char_T cv66[21] = { 'C', 'o', 'd', 'e', 'r', ':', 'M', 'A', 'T', 'L', 'A', 'B', ':', 'p', 'm', 'a', 'x', 's', 'i', 'z', 'e' }; emxArray_real_T *I; creal_T d_A; creal_T e_A; creal_T tmp; creal_T s; boolean_T b_jcol; int32_T j; real_T y_re; real_T y_im; real_T b_info; boolean_T b4; boolean_T b_ihi; boolean_T b5; emlrtHeapReferenceStackEnterFcnR2012b(emlrtRootTLSGlobal); *info = 0.0; n = A->size[0]; m = alpha1->size[0]; alpha1->size[0] = A->size[0]; emxEnsureCapacity((emxArray__common *)alpha1, m, (int32_T)sizeof(creal_T), &xb_emlrtRTEI); q = A->size[0]; for (m = 0; m < q; m++) { alpha1->data[m].re = 0.0; alpha1->data[m].im = 0.0; } m = beta1->size[0]; beta1->size[0] = A->size[0]; emxEnsureCapacity((emxArray__common *)beta1, m, (int32_T)sizeof(creal_T), &xb_emlrtRTEI); q = A->size[0]; for (m = 0; m < q; m++) { beta1->data[m].re = 0.0; beta1->data[m].im = 0.0; } jcol = A->size[0]; m = V->size[0] * V->size[1]; V->size[0] = jcol; emxEnsureCapacity((emxArray__common *)V, m, (int32_T)sizeof(creal_T), &xb_emlrtRTEI); q = A->size[0]; m = V->size[0] * V->size[1]; V->size[1] = q; emxEnsureCapacity((emxArray__common *)V, m, (int32_T)sizeof(creal_T), &xb_emlrtRTEI); q = A->size[0] * A->size[0]; for (m = 0; m < q; m++) { V->data[m].re = 0.0; V->data[m].im = 0.0; } if ((A->size[0] == 0) || (A->size[1] == 0)) { } else { anrm = 0.0; jcol = 0; exitg1 = FALSE; while ((exitg1 == FALSE) && (jcol <= A->size[0] * A->size[1] - 1)) { absxk = muDoubleScalarHypot(A->data[(int32_T)(1.0 + (real_T)jcol) - 1].re, A->data[(int32_T)(1.0 + (real_T)jcol) - 1].im); if (muDoubleScalarIsNaN(absxk)) { anrm = rtNaN; exitg1 = TRUE; } else { if (absxk > anrm) { anrm = absxk; } jcol++; } } if (!((!muDoubleScalarIsInf(anrm)) && (!muDoubleScalarIsNaN(anrm)))) { m = alpha1->size[0]; alpha1->size[0] = A->size[0]; emxEnsureCapacity((emxArray__common *)alpha1, m, (int32_T)sizeof(creal_T), &xb_emlrtRTEI); q = A->size[0]; for (m = 0; m < q; m++) { alpha1->data[m].re = rtNaN; alpha1->data[m].im = 0.0; } m = beta1->size[0]; beta1->size[0] = A->size[0]; emxEnsureCapacity((emxArray__common *)beta1, m, (int32_T)sizeof(creal_T), &xb_emlrtRTEI); q = A->size[0]; for (m = 0; m < q; m++) { beta1->data[m].re = rtNaN; beta1->data[m].im = 0.0; } jcol = A->size[0]; m = V->size[0] * V->size[1]; V->size[0] = jcol; emxEnsureCapacity((emxArray__common *)V, m, (int32_T)sizeof(creal_T), &xb_emlrtRTEI); q = A->size[0]; m = V->size[0] * V->size[1]; V->size[1] = q; emxEnsureCapacity((emxArray__common *)V, m, (int32_T)sizeof(creal_T), &xb_emlrtRTEI); q = A->size[0] * A->size[0]; for (m = 0; m < q; m++) { V->data[m].re = rtNaN; V->data[m].im = 0.0; } } else { ilascl = FALSE; anrmto = anrm; if ((anrm > 0.0) && (anrm < 6.7178761075670888E-139)) { anrmto = 6.7178761075670888E-139; ilascl = TRUE; } else { if (anrm > 1.4885657073574029E+138) { anrmto = 1.4885657073574029E+138; ilascl = TRUE; } } if (ilascl) { absxk = anrm; ctoc = anrmto; notdone = TRUE; while (notdone) { cfrom1 = absxk * 2.0041683600089728E-292; cto1 = ctoc / 4.9896007738368E+291; if ((cfrom1 > ctoc) && (ctoc != 0.0)) { mul = 2.0041683600089728E-292; absxk = cfrom1; } else if (cto1 > absxk) { mul = 4.9896007738368E+291; ctoc = cto1; } else { mul = ctoc / absxk; notdone = FALSE; } m = A->size[0] * A->size[1]; emxEnsureCapacity((emxArray__common *)A, m, (int32_T)sizeof(creal_T), &xb_emlrtRTEI); jcol = A->size[0]; q = A->size[1]; q *= jcol; for (m = 0; m < q; m++) { A->data[m].re *= mul; A->data[m].im *= mul; } } } emxInit_creal_T(&b_A, 2, &xb_emlrtRTEI, TRUE); emlrtPushRtStackR2012b(&xi_emlrtRSI, emlrtRootTLSGlobal); m = b_A->size[0] * b_A->size[1]; b_A->size[0] = A->size[0]; b_A->size[1] = A->size[1]; emxEnsureCapacity((emxArray__common *)b_A, m, (int32_T)sizeof(creal_T), &xb_emlrtRTEI); q = A->size[0] * A->size[1]; for (m = 0; m < q; m++) { b_A->data[m] = A->data[m]; } emxInit_int32_T(&rscale, 1, &xb_emlrtRTEI, TRUE); emxInit_creal_T(&c_A, 2, &xb_emlrtRTEI, TRUE); eml_matlab_zggbal(b_A, &ilo, &ihi, rscale); emlrtPopRtStackR2012b(&xi_emlrtRSI, emlrtRootTLSGlobal); emlrtPushRtStackR2012b(&yi_emlrtRSI, emlrtRootTLSGlobal); m = c_A->size[0] * c_A->size[1]; c_A->size[0] = b_A->size[0]; c_A->size[1] = b_A->size[1]; emxEnsureCapacity((emxArray__common *)c_A, m, (int32_T)sizeof(creal_T), &xb_emlrtRTEI); q = b_A->size[0] * b_A->size[1]; for (m = 0; m < q; m++) { c_A->data[m] = b_A->data[m]; } emlrtPushRtStackR2012b(&mj_emlrtRSI, emlrtRootTLSGlobal); emlrtPushRtStackR2012b(&rj_emlrtRSI, emlrtRootTLSGlobal); emlrtPushRtStackR2012b(&sj_emlrtRSI, emlrtRootTLSGlobal); if (b_A->size[0] <= 0) { absxk = 0.0; } else { absxk = b_A->size[0]; } if (b_A->size[0] <= 0) { absxk = 0.0; } else { absxk *= (real_T)b_A->size[0]; } if (2.147483647E+9 >= absxk) { } else { emlrtPushRtStackR2012b(&uj_emlrtRSI, emlrtRootTLSGlobal); y = NULL; m13 = mxCreateCharArray(2, iv53); for (i = 0; i < 21; i++) { cv65[i] = cv66[i]; } emlrtInitCharArrayR2013a(emlrtRootTLSGlobal, 21, m13, cv65); emlrtAssign(&y, m13); error(message(y, &rb_emlrtMCI), &sb_emlrtMCI); emlrtPopRtStackR2012b(&uj_emlrtRSI, emlrtRootTLSGlobal); } b_emxInit_real_T(&I, 2, &xb_emlrtRTEI, TRUE); emlrtPopRtStackR2012b(&sj_emlrtRSI, emlrtRootTLSGlobal); jcol = b_A->size[0]; m = I->size[0] * I->size[1]; I->size[0] = jcol; emxEnsureCapacity((emxArray__common *)I, m, (int32_T)sizeof(real_T), &xb_emlrtRTEI); q = b_A->size[0]; m = I->size[0] * I->size[1]; I->size[1] = q; emxEnsureCapacity((emxArray__common *)I, m, (int32_T)sizeof(real_T), &xb_emlrtRTEI); q = b_A->size[0] * b_A->size[0]; for (m = 0; m < q; m++) { I->data[m] = 0.0; } q = muIntScalarMin_sint32(b_A->size[0], b_A->size[0]); if (q > 0) { emlrtPushRtStackR2012b(&tj_emlrtRSI, emlrtRootTLSGlobal); if (q > 2147483646) { emlrtPushRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); check_forloop_overflow_error(); emlrtPopRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); } emlrtPopRtStackR2012b(&tj_emlrtRSI, emlrtRootTLSGlobal); for (i = 0; i + 1 <= q; i++) { I->data[i + I->size[0] * i] = 1.0; } } emlrtPopRtStackR2012b(&rj_emlrtRSI, emlrtRootTLSGlobal); m = V->size[0] * V->size[1]; V->size[0] = I->size[0]; V->size[1] = I->size[1]; emxEnsureCapacity((emxArray__common *)V, m, (int32_T)sizeof(creal_T), &xb_emlrtRTEI); q = I->size[0] * I->size[1]; for (m = 0; m < q; m++) { V->data[m].re = I->data[m]; V->data[m].im = 0.0; } emxFree_real_T(&I); emlrtPopRtStackR2012b(&mj_emlrtRSI, emlrtRootTLSGlobal); if ((b_A->size[0] <= 1) || (ihi < ilo + 2)) { } else { for (jcol = ilo - 1; jcol + 1 < ihi - 1; jcol++) { for (q = ihi - 1; q + 1 > jcol + 2; q--) { emlrtPushRtStackR2012b(&nj_emlrtRSI, emlrtRootTLSGlobal); d_A = c_A->data[(q + c_A->size[0] * jcol) - 1]; e_A = c_A->data[q + c_A->size[0] * jcol]; eml_matlab_zlartg(d_A, e_A, &mul, &s, &tmp); emlrtPopRtStackR2012b(&nj_emlrtRSI, emlrtRootTLSGlobal); c_A->data[(q + c_A->size[0] * jcol) - 1] = tmp; c_A->data[q + c_A->size[0] * jcol].re = 0.0; c_A->data[q + c_A->size[0] * jcol].im = 0.0; emlrtPushRtStackR2012b(&oj_emlrtRSI, emlrtRootTLSGlobal); emlrtPushRtStackR2012b(&ak_emlrtRSI, emlrtRootTLSGlobal); if (jcol + 2 > ihi) { b_jcol = FALSE; } else { b_jcol = (ihi > 2147483646); } if (b_jcol) { emlrtPushRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); check_forloop_overflow_error(); emlrtPopRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); } emlrtPopRtStackR2012b(&ak_emlrtRSI, emlrtRootTLSGlobal); for (j = jcol + 1; j + 1 <= ihi; j++) { tmp.re = mul * c_A->data[(q + c_A->size[0] * j) - 1].re; tmp.im = mul * c_A->data[(q + c_A->size[0] * j) - 1].im; y_re = s.re * c_A->data[q + c_A->size[0] * j].re - s.im * c_A->data[q + c_A->size[0] * j].im; y_im = s.re * c_A->data[q + c_A->size[0] * j].im + s.im * c_A->data[q + c_A->size[0] * j].re; absxk = c_A->data[(q + c_A->size[0] * j) - 1].re; ctoc = c_A->data[(q + c_A->size[0] * j) - 1].im; cfrom1 = c_A->data[(q + c_A->size[0] * j) - 1].im; cto1 = c_A->data[(q + c_A->size[0] * j) - 1].re; c_A->data[q + c_A->size[0] * j].re = mul * c_A->data[q + c_A-> size[0] * j].re - (s.re * absxk + s.im * ctoc); c_A->data[q + c_A->size[0] * j].im = mul * c_A->data[q + c_A-> size[0] * j].im - (s.re * cfrom1 - s.im * cto1); c_A->data[(q + c_A->size[0] * j) - 1].re = tmp.re + y_re; c_A->data[(q + c_A->size[0] * j) - 1].im = tmp.im + y_im; } emlrtPopRtStackR2012b(&oj_emlrtRSI, emlrtRootTLSGlobal); s.re = -s.re; s.im = -s.im; emlrtPushRtStackR2012b(&pj_emlrtRSI, emlrtRootTLSGlobal); emlrtPushRtStackR2012b(&bk_emlrtRSI, emlrtRootTLSGlobal); emlrtPopRtStackR2012b(&bk_emlrtRSI, emlrtRootTLSGlobal); for (i = ilo - 1; i + 1 <= ihi; i++) { tmp.re = mul * c_A->data[i + c_A->size[0] * q].re; tmp.im = mul * c_A->data[i + c_A->size[0] * q].im; y_re = s.re * c_A->data[i + c_A->size[0] * (q - 1)].re - s.im * c_A->data[i + c_A->size[0] * (q - 1)].im; y_im = s.re * c_A->data[i + c_A->size[0] * (q - 1)].im + s.im * c_A->data[i + c_A->size[0] * (q - 1)].re; absxk = c_A->data[i + c_A->size[0] * q].re; ctoc = c_A->data[i + c_A->size[0] * q].im; cfrom1 = c_A->data[i + c_A->size[0] * q].im; cto1 = c_A->data[i + c_A->size[0] * q].re; c_A->data[i + c_A->size[0] * (q - 1)].re = mul * c_A->data[i + c_A->size[0] * (q - 1)].re - (s.re * absxk + s.im * ctoc); c_A->data[i + c_A->size[0] * (q - 1)].im = mul * c_A->data[i + c_A->size[0] * (q - 1)].im - (s.re * cfrom1 - s.im * cto1); c_A->data[i + c_A->size[0] * q].re = tmp.re + y_re; c_A->data[i + c_A->size[0] * q].im = tmp.im + y_im; } emlrtPopRtStackR2012b(&pj_emlrtRSI, emlrtRootTLSGlobal); emlrtPushRtStackR2012b(&qj_emlrtRSI, emlrtRootTLSGlobal); emlrtPushRtStackR2012b(&bk_emlrtRSI, emlrtRootTLSGlobal); if (1 > b_A->size[0]) { notdone = FALSE; } else { notdone = (b_A->size[0] > 2147483646); } if (notdone) { emlrtPushRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); check_forloop_overflow_error(); emlrtPopRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); } emlrtPopRtStackR2012b(&bk_emlrtRSI, emlrtRootTLSGlobal); for (i = 0; i + 1 <= b_A->size[0]; i++) { tmp.re = mul * V->data[i + V->size[0] * q].re; tmp.im = mul * V->data[i + V->size[0] * q].im; y_re = s.re * V->data[i + V->size[0] * (q - 1)].re - s.im * V->data[i + V->size[0] * (q - 1)].im; y_im = s.re * V->data[i + V->size[0] * (q - 1)].im + s.im * V->data[i + V->size[0] * (q - 1)].re; absxk = V->data[i + V->size[0] * q].re; ctoc = V->data[i + V->size[0] * q].im; cfrom1 = V->data[i + V->size[0] * q].im; cto1 = V->data[i + V->size[0] * q].re; V->data[i + V->size[0] * (q - 1)].re = mul * V->data[i + V->size[0] * (q - 1)].re - (s.re * absxk + s.im * ctoc); V->data[i + V->size[0] * (q - 1)].im = mul * V->data[i + V->size[0] * (q - 1)].im - (s.re * cfrom1 - s.im * cto1); V->data[i + V->size[0] * q].re = tmp.re + y_re; V->data[i + V->size[0] * q].im = tmp.im + y_im; } emlrtPopRtStackR2012b(&qj_emlrtRSI, emlrtRootTLSGlobal); } } } emxFree_creal_T(&b_A); emlrtPopRtStackR2012b(&yi_emlrtRSI, emlrtRootTLSGlobal); emlrtPushRtStackR2012b(&aj_emlrtRSI, emlrtRootTLSGlobal); eml_matlab_zhgeqz(c_A, ilo, ihi, V, &b_info, alpha1, beta1); emlrtPopRtStackR2012b(&aj_emlrtRSI, emlrtRootTLSGlobal); *info = b_info; if (b_info != 0.0) { } else { emlrtPushRtStackR2012b(&bj_emlrtRSI, emlrtRootTLSGlobal); eml_matlab_ztgevc(c_A, V); emlrtPopRtStackR2012b(&bj_emlrtRSI, emlrtRootTLSGlobal); emlrtPushRtStackR2012b(&cj_emlrtRSI, emlrtRootTLSGlobal); q = V->size[0]; m = V->size[1]; if (ilo > 1) { for (i = ilo - 2; i + 1 >= 1; i--) { if (rscale->data[i] != i + 1) { emlrtPushRtStackR2012b(&pk_emlrtRSI, emlrtRootTLSGlobal); if (1 > m) { b4 = FALSE; } else { b4 = (m > 2147483646); } if (b4) { emlrtPushRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); check_forloop_overflow_error(); emlrtPopRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); } emlrtPopRtStackR2012b(&pk_emlrtRSI, emlrtRootTLSGlobal); for (j = 0; j + 1 <= m; j++) { tmp = V->data[i + V->size[0] * j]; V->data[i + V->size[0] * j] = V->data[(rscale->data[i] + V-> size[0] * j) - 1]; V->data[(rscale->data[i] + V->size[0] * j) - 1] = tmp; } } } } if (ihi < q) { emlrtPushRtStackR2012b(&qk_emlrtRSI, emlrtRootTLSGlobal); if (ihi + 1 > q) { b_ihi = FALSE; } else { b_ihi = (q > 2147483646); } if (b_ihi) { emlrtPushRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); check_forloop_overflow_error(); emlrtPopRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); } emlrtPopRtStackR2012b(&qk_emlrtRSI, emlrtRootTLSGlobal); while (ihi + 1 <= q) { if (rscale->data[ihi] != ihi + 1) { emlrtPushRtStackR2012b(&rk_emlrtRSI, emlrtRootTLSGlobal); if (1 > m) { b5 = FALSE; } else { b5 = (m > 2147483646); } if (b5) { emlrtPushRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); check_forloop_overflow_error(); emlrtPopRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); } emlrtPopRtStackR2012b(&rk_emlrtRSI, emlrtRootTLSGlobal); for (j = 0; j + 1 <= m; j++) { tmp = V->data[ihi + V->size[0] * j]; V->data[ihi + V->size[0] * j] = V->data[(rscale->data[ihi] + V->size[0] * j) - 1]; V->data[(rscale->data[ihi] + V->size[0] * j) - 1] = tmp; } } ihi++; } } emlrtPopRtStackR2012b(&cj_emlrtRSI, emlrtRootTLSGlobal); for (jcol = 0; jcol < n; jcol++) { absxk = muDoubleScalarAbs(V->data[V->size[0] * jcol].re) + muDoubleScalarAbs(V->data[V->size[0] * jcol].im); if (n > 1) { for (q = 1; q - 1 <= n - 2; q++) { ctoc = muDoubleScalarAbs(V->data[q + V->size[0] * jcol].re) + muDoubleScalarAbs(V->data[q + V->size[0] * jcol].im); if (ctoc > absxk) { absxk = ctoc; } } } if (absxk >= 6.7178761075670888E-139) { absxk = 1.0 / absxk; for (q = 0; q < n; q++) { V->data[q + V->size[0] * jcol].re *= absxk; V->data[q + V->size[0] * jcol].im *= absxk; } } } if (ilascl) { notdone = TRUE; while (notdone) { cfrom1 = anrmto * 2.0041683600089728E-292; cto1 = anrm / 4.9896007738368E+291; if ((cfrom1 > anrm) && (anrm != 0.0)) { mul = 2.0041683600089728E-292; anrmto = cfrom1; } else if (cto1 > anrmto) { mul = 4.9896007738368E+291; anrm = cto1; } else { mul = anrm / anrmto; notdone = FALSE; } m = alpha1->size[0]; emxEnsureCapacity((emxArray__common *)alpha1, m, (int32_T)sizeof (creal_T), &xb_emlrtRTEI); q = alpha1->size[0]; for (m = 0; m < q; m++) { alpha1->data[m].re *= mul; alpha1->data[m].im *= mul; } } } } emxFree_creal_T(&c_A); emxFree_int32_T(&rscale); } } emlrtHeapReferenceStackLeaveFcnR2012b(emlrtRootTLSGlobal); } static void eml_matlab_zhgeqz(emxArray_creal_T *A, int32_T ilo, int32_T ihi, emxArray_creal_T *Z, real_T *info, emxArray_creal_T *alpha1, emxArray_creal_T * beta1) { boolean_T b18; int32_T n; int32_T i; int32_T loop_ub; real_T eshift_re; real_T eshift_im; creal_T ctemp; real_T rho_re; real_T rho_im; real_T anorm; real_T temp; real_T b_atol; real_T ascale; boolean_T failed; boolean_T overflow; int32_T j; boolean_T guard1 = FALSE; boolean_T guard2 = FALSE; int32_T ifirst; int32_T istart; int32_T ilast; int32_T ilastm1; int32_T ifrstm; int32_T ilastm; int32_T iiter; int32_T maxit; boolean_T goto60; boolean_T goto70; boolean_T goto90; boolean_T b19; int32_T jiter; int32_T exitg1; boolean_T exitg3; int32_T jp1; int32_T jm1; boolean_T b_guard1 = FALSE; creal_T b_A; creal_T t1; creal_T d; creal_T sigma1; real_T sigma2_re; real_T sigma2_im; boolean_T exitg2; real_T tempr; creal_T c_A; creal_T d_A; boolean_T b_j; boolean_T b_ifrstm; boolean_T b20; boolean_T b21; boolean_T b22; b18 = !((Z->size[0] == 0) || (Z->size[1] == 0)); n = A->size[0]; i = alpha1->size[0]; alpha1->size[0] = A->size[0]; emxEnsureCapacity((emxArray__common *)alpha1, i, (int32_T)sizeof(creal_T), &ac_emlrtRTEI); loop_ub = A->size[0]; for (i = 0; i < loop_ub; i++) { alpha1->data[i].re = 0.0; alpha1->data[i].im = 0.0; } i = beta1->size[0]; beta1->size[0] = A->size[0]; emxEnsureCapacity((emxArray__common *)beta1, i, (int32_T)sizeof(creal_T), &ac_emlrtRTEI); loop_ub = A->size[0]; for (i = 0; i < loop_ub; i++) { beta1->data[i].re = 1.0; beta1->data[i].im = 0.0; } eshift_re = 0.0; eshift_im = 0.0; ctemp.re = 0.0; ctemp.im = 0.0; rho_re = 0.0; rho_im = 0.0; emlrtPushRtStackR2012b(&lk_emlrtRSI, emlrtRootTLSGlobal); anorm = eml_matlab_zlanhs(A, ilo, ihi); emlrtPopRtStackR2012b(&lk_emlrtRSI, emlrtRootTLSGlobal); temp = 2.2204460492503131E-16 * anorm; b_atol = 2.2250738585072014E-308; if (temp > 2.2250738585072014E-308) { b_atol = temp; } temp = 2.2250738585072014E-308; if (anorm > 2.2250738585072014E-308) { temp = anorm; } ascale = 1.0 / temp; failed = TRUE; emlrtPushRtStackR2012b(&kk_emlrtRSI, emlrtRootTLSGlobal); if (ihi + 1 > A->size[0]) { overflow = FALSE; } else { overflow = (A->size[0] > 2147483646); } if (overflow) { emlrtPushRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); check_forloop_overflow_error(); emlrtPopRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); } emlrtPopRtStackR2012b(&kk_emlrtRSI, emlrtRootTLSGlobal); for (j = ihi; j + 1 <= A->size[0]; j++) { alpha1->data[j] = A->data[j + A->size[0] * j]; } guard1 = FALSE; guard2 = FALSE; if (ihi >= ilo) { ifirst = ilo; istart = ilo; ilast = ihi - 1; ilastm1 = ihi - 2; if (b18) { ifrstm = 1; ilastm = A->size[0]; } else { ifrstm = ilo; ilastm = ihi; } iiter = 0; maxit = 30 * ((ihi - ilo) + 1); goto60 = FALSE; goto70 = FALSE; goto90 = FALSE; emlrtPushRtStackR2012b(&jk_emlrtRSI, emlrtRootTLSGlobal); if (1 > maxit) { b19 = FALSE; } else { b19 = (maxit > 2147483646); } if (b19) { emlrtPushRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); check_forloop_overflow_error(); emlrtPopRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); } emlrtPopRtStackR2012b(&jk_emlrtRSI, emlrtRootTLSGlobal); jiter = 1; do { exitg1 = 0; if (jiter <= maxit) { if (ilast + 1 == ilo) { goto60 = TRUE; } else if (muDoubleScalarAbs(A->data[ilast + A->size[0] * ilastm1].re) + muDoubleScalarAbs(A->data[ilast + A->size[0] * ilastm1].im) <= b_atol) { A->data[ilast + A->size[0] * ilastm1].re = 0.0; A->data[ilast + A->size[0] * ilastm1].im = 0.0; goto60 = TRUE; } else { j = ilastm1; exitg3 = FALSE; while ((exitg3 == FALSE) && (j + 1 >= ilo)) { if (j + 1 == ilo) { overflow = TRUE; } else if (muDoubleScalarAbs(A->data[j + A->size[0] * (j - 1)].re) + muDoubleScalarAbs(A->data[j + A->size[0] * (j - 1)].im) <= b_atol) { A->data[j + A->size[0] * (j - 1)].re = 0.0; A->data[j + A->size[0] * (j - 1)].im = 0.0; overflow = TRUE; } else { overflow = FALSE; } if (overflow) { ifirst = j + 1; goto70 = TRUE; exitg3 = TRUE; } else { j--; } } } if (goto60 || goto70) { overflow = TRUE; } else { overflow = FALSE; } if (!overflow) { jp1 = alpha1->size[0]; i = alpha1->size[0]; alpha1->size[0] = jp1; emxEnsureCapacity((emxArray__common *)alpha1, i, (int32_T)sizeof (creal_T), &ac_emlrtRTEI); for (i = 0; i < jp1; i++) { alpha1->data[i].re = rtNaN; alpha1->data[i].im = 0.0; } jp1 = beta1->size[0]; i = beta1->size[0]; beta1->size[0] = jp1; emxEnsureCapacity((emxArray__common *)beta1, i, (int32_T)sizeof (creal_T), &ac_emlrtRTEI); for (i = 0; i < jp1; i++) { beta1->data[i].re = rtNaN; beta1->data[i].im = 0.0; } if (b18) { i = Z->size[0] * Z->size[1]; emxEnsureCapacity((emxArray__common *)Z, i, (int32_T)sizeof(creal_T), &ac_emlrtRTEI); loop_ub = Z->size[1]; for (i = 0; i < loop_ub; i++) { jp1 = Z->size[0]; for (jm1 = 0; jm1 < jp1; jm1++) { Z->data[jm1 + Z->size[0] * i].re = rtNaN; Z->data[jm1 + Z->size[0] * i].im = 0.0; } } } *info = -1.0; exitg1 = 1; } else { b_guard1 = FALSE; if (goto60) { goto60 = FALSE; alpha1->data[ilast] = A->data[ilast + A->size[0] * ilast]; ilast = ilastm1; ilastm1--; if (ilast + 1 < ilo) { failed = FALSE; guard2 = TRUE; exitg1 = 1; } else { iiter = 0; eshift_re = 0.0; eshift_im = 0.0; if (!b18) { ilastm = ilast + 1; if (ifrstm > ilast + 1) { ifrstm = ilo; } } b_guard1 = TRUE; } } else { if (goto70) { goto70 = FALSE; iiter++; if (!b18) { ifrstm = ifirst; } if (iiter - div_s32(iiter, 10) * 10 != 0) { b_A.re = -(A->data[ilast + A->size[0] * ilast].re - A-> data[ilastm1 + A->size[0] * ilastm1].re); b_A.im = -(A->data[ilast + A->size[0] * ilast].im - A-> data[ilastm1 + A->size[0] * ilastm1].im); t1 = eml_div(b_A, 2.0); temp = A->data[ilastm1 + A->size[0] * ilast].re * A->data[ilast + A->size[0] * ilastm1].re - A->data[ilastm1 + A->size[0] * ilast].im * A->data[ilast + A->size[0] * ilastm1].im; anorm = A->data[ilastm1 + A->size[0] * ilast].re * A->data[ilast + A->size[0] * ilastm1].im + A->data[ilastm1 + A->size[0] * ilast].im * A->data[ilast + A->size[0] * ilastm1].re; d.re = (t1.re * t1.re - t1.im * t1.im) + temp; d.im = (t1.re * t1.im + t1.im * t1.re) + anorm; e_sqrt(&d); sigma1.re = A->data[ilastm1 + A->size[0] * ilastm1].re - (t1.re - d.re); sigma1.im = A->data[ilastm1 + A->size[0] * ilastm1].im - (t1.im - d.im); sigma2_re = A->data[ilastm1 + A->size[0] * ilastm1].re - (t1.re + d.re); sigma2_im = A->data[ilastm1 + A->size[0] * ilastm1].im - (t1.im + d.im); rho_re = sigma1.re - A->data[ilast + A->size[0] * ilast].re; rho_im = sigma1.im - A->data[ilast + A->size[0] * ilast].im; temp = sigma2_re - A->data[ilast + A->size[0] * ilast].re; anorm = sigma2_im - A->data[ilast + A->size[0] * ilast].im; if (muDoubleScalarHypot(rho_re, rho_im) <= muDoubleScalarHypot (temp, anorm)) { sigma2_re = sigma1.re; sigma2_im = sigma1.im; rho_re = t1.re - d.re; rho_im = t1.im - d.im; } else { rho_re = t1.re + d.re; rho_im = t1.im + d.im; } } else { eshift_re += A->data[ilast + A->size[0] * ilastm1].re; eshift_im += A->data[ilast + A->size[0] * ilastm1].im; sigma2_re = eshift_re; sigma2_im = eshift_im; } j = ilastm1; jp1 = ilastm1 + 1; exitg2 = FALSE; while ((exitg2 == FALSE) && (j + 1 > ifirst)) { istart = j + 1; ctemp.re = A->data[j + A->size[0] * j].re - sigma2_re; ctemp.im = A->data[j + A->size[0] * j].im - sigma2_im; temp = ascale * (muDoubleScalarAbs(ctemp.re) + muDoubleScalarAbs (ctemp.im)); anorm = ascale * (muDoubleScalarAbs(A->data[jp1 + A->size[0] * j] .re) + muDoubleScalarAbs(A->data[jp1 + A->size[0] * j].im)); tempr = temp; if (anorm > temp) { tempr = anorm; } if ((tempr < 1.0) && (tempr != 0.0)) { temp /= tempr; anorm /= tempr; } if ((muDoubleScalarAbs(A->data[j + A->size[0] * (j - 1)].re) + muDoubleScalarAbs(A->data[j + A->size[0] * (j - 1)].im)) * anorm <= temp * b_atol) { goto90 = TRUE; exitg2 = TRUE; } else { jp1 = j; j--; } } if (!goto90) { istart = ifirst; if (ifirst == ilastm1 + 1) { ctemp.re = rho_re; ctemp.im = rho_im; } else { ctemp.re = A->data[(ifirst + A->size[0] * (ifirst - 1)) - 1]. re - sigma2_re; ctemp.im = A->data[(ifirst + A->size[0] * (ifirst - 1)) - 1]. im - sigma2_im; } goto90 = TRUE; } } if (goto90) { goto90 = FALSE; emlrtPushRtStackR2012b(&ik_emlrtRSI, emlrtRootTLSGlobal); c_A = A->data[istart + A->size[0] * (istart - 1)]; b_eml_matlab_zlartg(ctemp, c_A, &sigma2_im, &sigma1); emlrtPopRtStackR2012b(&ik_emlrtRSI, emlrtRootTLSGlobal); j = istart; jm1 = istart - 2; while (j < ilast + 1) { if (j > istart) { emlrtPushRtStackR2012b(&hk_emlrtRSI, emlrtRootTLSGlobal); c_A = A->data[(j + A->size[0] * jm1) - 1]; d_A = A->data[j + A->size[0] * jm1]; eml_matlab_zlartg(c_A, d_A, &sigma2_im, &sigma1, &t1); emlrtPopRtStackR2012b(&hk_emlrtRSI, emlrtRootTLSGlobal); A->data[(j + A->size[0] * jm1) - 1] = t1; A->data[j + A->size[0] * jm1].re = 0.0; A->data[j + A->size[0] * jm1].im = 0.0; } emlrtPushRtStackR2012b(&gk_emlrtRSI, emlrtRootTLSGlobal); emlrtPushRtStackR2012b(&ak_emlrtRSI, emlrtRootTLSGlobal); if (j > ilastm) { b_j = FALSE; } else { b_j = (ilastm > 2147483646); } if (b_j) { emlrtPushRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); check_forloop_overflow_error(); emlrtPopRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); } emlrtPopRtStackR2012b(&ak_emlrtRSI, emlrtRootTLSGlobal); for (loop_ub = j - 1; loop_ub + 1 <= ilastm; loop_ub++) { t1.re = sigma2_im * A->data[(j + A->size[0] * loop_ub) - 1].re; t1.im = sigma2_im * A->data[(j + A->size[0] * loop_ub) - 1].im; d.re = sigma1.re * A->data[j + A->size[0] * loop_ub].re - sigma1.im * A->data[j + A->size[0] * loop_ub].im; d.im = sigma1.re * A->data[j + A->size[0] * loop_ub].im + sigma1.im * A->data[j + A->size[0] * loop_ub].re; temp = A->data[(j + A->size[0] * loop_ub) - 1].re; anorm = A->data[(j + A->size[0] * loop_ub) - 1].im; tempr = A->data[(j + A->size[0] * loop_ub) - 1].im; sigma2_re = A->data[(j + A->size[0] * loop_ub) - 1].re; A->data[j + A->size[0] * loop_ub].re = sigma2_im * A->data[j + A->size[0] * loop_ub].re - (sigma1.re * temp + sigma1.im * anorm); A->data[j + A->size[0] * loop_ub].im = sigma2_im * A->data[j + A->size[0] * loop_ub].im - (sigma1.re * tempr - sigma1.im * sigma2_re); A->data[(j + A->size[0] * loop_ub) - 1].re = t1.re + d.re; A->data[(j + A->size[0] * loop_ub) - 1].im = t1.im + d.im; } emlrtPopRtStackR2012b(&gk_emlrtRSI, emlrtRootTLSGlobal); sigma1.re = -sigma1.re; sigma1.im = -sigma1.im; emlrtPushRtStackR2012b(&fk_emlrtRSI, emlrtRootTLSGlobal); loop_ub = j + 2; if (ilast + 1 < j + 2) { loop_ub = ilast + 1; } emlrtPushRtStackR2012b(&bk_emlrtRSI, emlrtRootTLSGlobal); if (ifrstm > loop_ub) { b_ifrstm = FALSE; } else { b_ifrstm = (loop_ub > 2147483646); } if (b_ifrstm) { emlrtPushRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); check_forloop_overflow_error(); emlrtPopRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); } emlrtPopRtStackR2012b(&bk_emlrtRSI, emlrtRootTLSGlobal); for (i = ifrstm - 1; i + 1 <= loop_ub; i++) { t1.re = sigma2_im * A->data[i + A->size[0] * j].re; t1.im = sigma2_im * A->data[i + A->size[0] * j].im; d.re = sigma1.re * A->data[i + A->size[0] * (j - 1)].re - sigma1.im * A->data[i + A->size[0] * (j - 1)].im; d.im = sigma1.re * A->data[i + A->size[0] * (j - 1)].im + sigma1.im * A->data[i + A->size[0] * (j - 1)].re; temp = A->data[i + A->size[0] * j].re; anorm = A->data[i + A->size[0] * j].im; tempr = A->data[i + A->size[0] * j].im; sigma2_re = A->data[i + A->size[0] * j].re; A->data[i + A->size[0] * (j - 1)].re = sigma2_im * A->data[i + A->size[0] * (j - 1)].re - (sigma1.re * temp + sigma1.im * anorm); A->data[i + A->size[0] * (j - 1)].im = sigma2_im * A->data[i + A->size[0] * (j - 1)].im - (sigma1.re * tempr - sigma1.im * sigma2_re); A->data[i + A->size[0] * j].re = t1.re + d.re; A->data[i + A->size[0] * j].im = t1.im + d.im; } emlrtPopRtStackR2012b(&fk_emlrtRSI, emlrtRootTLSGlobal); if (b18) { emlrtPushRtStackR2012b(&ek_emlrtRSI, emlrtRootTLSGlobal); emlrtPushRtStackR2012b(&bk_emlrtRSI, emlrtRootTLSGlobal); if (1 > n) { b20 = FALSE; } else { b20 = (n > 2147483646); } if (b20) { emlrtPushRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); check_forloop_overflow_error(); emlrtPopRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); } emlrtPopRtStackR2012b(&bk_emlrtRSI, emlrtRootTLSGlobal); for (i = 0; i + 1 <= n; i++) { t1.re = sigma2_im * Z->data[i + Z->size[0] * j].re; t1.im = sigma2_im * Z->data[i + Z->size[0] * j].im; d.re = sigma1.re * Z->data[i + Z->size[0] * (j - 1)].re - sigma1.im * Z->data[i + Z->size[0] * (j - 1)].im; d.im = sigma1.re * Z->data[i + Z->size[0] * (j - 1)].im + sigma1.im * Z->data[i + Z->size[0] * (j - 1)].re; anorm = Z->data[i + Z->size[0] * j].re; temp = Z->data[i + Z->size[0] * j].im; tempr = Z->data[i + Z->size[0] * j].im; sigma2_re = Z->data[i + Z->size[0] * j].re; Z->data[i + Z->size[0] * (j - 1)].re = sigma2_im * Z->data[i + Z->size[0] * (j - 1)].re - (sigma1.re * anorm + sigma1.im * temp); Z->data[i + Z->size[0] * (j - 1)].im = sigma2_im * Z->data[i + Z->size[0] * (j - 1)].im - (sigma1.re * tempr - sigma1.im * sigma2_re); Z->data[i + Z->size[0] * j].re = t1.re + d.re; Z->data[i + Z->size[0] * j].im = t1.im + d.im; } emlrtPopRtStackR2012b(&ek_emlrtRSI, emlrtRootTLSGlobal); } jm1 = j - 1; j++; } } b_guard1 = TRUE; } if (b_guard1 == TRUE) { jiter++; } } } else { guard2 = TRUE; exitg1 = 1; } } while (exitg1 == 0); } else { guard1 = TRUE; } if (guard2 == TRUE) { if (failed) { *info = ilast + 1; emlrtPushRtStackR2012b(&dk_emlrtRSI, emlrtRootTLSGlobal); if (1 > ilast + 1) { b21 = FALSE; } else { b21 = (ilast + 1 > 2147483646); } if (b21) { emlrtPushRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); check_forloop_overflow_error(); emlrtPopRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); } emlrtPopRtStackR2012b(&dk_emlrtRSI, emlrtRootTLSGlobal); for (jp1 = 0; jp1 + 1 <= ilast + 1; jp1++) { alpha1->data[jp1].re = rtNaN; alpha1->data[jp1].im = 0.0; beta1->data[jp1].re = rtNaN; beta1->data[jp1].im = 0.0; } if (b18) { i = Z->size[0] * Z->size[1]; emxEnsureCapacity((emxArray__common *)Z, i, (int32_T)sizeof(creal_T), &ac_emlrtRTEI); loop_ub = Z->size[1]; for (i = 0; i < loop_ub; i++) { jp1 = Z->size[0]; for (jm1 = 0; jm1 < jp1; jm1++) { Z->data[jm1 + Z->size[0] * i].re = rtNaN; Z->data[jm1 + Z->size[0] * i].im = 0.0; } } } } else { guard1 = TRUE; } } if (guard1 == TRUE) { emlrtPushRtStackR2012b(&ck_emlrtRSI, emlrtRootTLSGlobal); if (1 > ilo - 1) { b22 = FALSE; } else { b22 = (ilo - 1 > 2147483646); } if (b22) { emlrtPushRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); check_forloop_overflow_error(); emlrtPopRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); } emlrtPopRtStackR2012b(&ck_emlrtRSI, emlrtRootTLSGlobal); for (j = 0; j + 1 < ilo; j++) { alpha1->data[j] = A->data[j + A->size[0] * j]; } *info = 0.0; } } static real_T eml_matlab_zlanhs(const emxArray_creal_T *A, int32_T ilo, int32_T ihi) { real_T f; real_T scale; real_T sumsq; boolean_T firstNonZero; boolean_T b_ilo; int32_T j; int32_T b; boolean_T c_ilo; int32_T i; real_T reAij; real_T imAij; real_T temp2; f = 0.0; if (ilo > ihi) { } else { scale = 0.0; sumsq = 0.0; firstNonZero = TRUE; emlrtPushRtStackR2012b(&mk_emlrtRSI, emlrtRootTLSGlobal); if (ilo > ihi) { b_ilo = FALSE; } else { b_ilo = (ihi > 2147483646); } if (b_ilo) { emlrtPushRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); check_forloop_overflow_error(); emlrtPopRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); } emlrtPopRtStackR2012b(&mk_emlrtRSI, emlrtRootTLSGlobal); for (j = ilo; j <= ihi; j++) { b = j + 1; if (ihi < j + 1) { b = ihi; } emlrtPushRtStackR2012b(&nk_emlrtRSI, emlrtRootTLSGlobal); if (ilo > b) { c_ilo = FALSE; } else { c_ilo = (b > 2147483646); } if (c_ilo) { emlrtPushRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); check_forloop_overflow_error(); emlrtPopRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); } emlrtPopRtStackR2012b(&nk_emlrtRSI, emlrtRootTLSGlobal); for (i = ilo; i <= b; i++) { reAij = A->data[(i + A->size[0] * (j - 1)) - 1].re; imAij = A->data[(i + A->size[0] * (j - 1)) - 1].im; if (reAij != 0.0) { reAij = muDoubleScalarAbs(reAij); if (firstNonZero) { sumsq = 1.0; scale = reAij; firstNonZero = FALSE; } else if (scale < reAij) { temp2 = scale / reAij; sumsq = 1.0 + sumsq * temp2 * temp2; scale = reAij; } else { temp2 = reAij / scale; sumsq += temp2 * temp2; } } if (imAij != 0.0) { reAij = muDoubleScalarAbs(imAij); if (firstNonZero) { sumsq = 1.0; scale = reAij; firstNonZero = FALSE; } else if (scale < reAij) { temp2 = scale / reAij; sumsq = 1.0 + sumsq * temp2 * temp2; scale = reAij; } else { temp2 = reAij / scale; sumsq += temp2 * temp2; } } } } emlrtPushRtStackR2012b(&ok_emlrtRSI, emlrtRootTLSGlobal); if (sumsq < 0.0) { emlrtPushRtStackR2012b(&md_emlrtRSI, emlrtRootTLSGlobal); e_eml_error(); emlrtPopRtStackR2012b(&md_emlrtRSI, emlrtRootTLSGlobal); } f = scale * muDoubleScalarSqrt(sumsq); emlrtPopRtStackR2012b(&ok_emlrtRSI, emlrtRootTLSGlobal); } return f; } static void eml_matlab_zlartg(const creal_T f, const creal_T g, real_T *cs, creal_T *sn, creal_T *r) { real_T scale; real_T f2s; real_T g2; real_T fs_re; real_T fs_im; real_T gs_re; real_T gs_im; int32_T count; int32_T rescaledir; boolean_T guard1 = FALSE; real_T g2s; boolean_T b6; boolean_T b7; scale = muDoubleScalarAbs(f.re); f2s = muDoubleScalarAbs(f.im); if (f2s > scale) { scale = f2s; } f2s = muDoubleScalarAbs(g.re); g2 = muDoubleScalarAbs(g.im); if (g2 > f2s) { f2s = g2; } if (f2s > scale) { scale = f2s; } fs_re = f.re; fs_im = f.im; gs_re = g.re; gs_im = g.im; count = 0; rescaledir = 0; guard1 = FALSE; if (scale >= 7.4428285367870146E+137) { do { count++; fs_re *= 1.3435752215134178E-138; fs_im *= 1.3435752215134178E-138; gs_re *= 1.3435752215134178E-138; gs_im *= 1.3435752215134178E-138; scale *= 1.3435752215134178E-138; } while (!(scale < 7.4428285367870146E+137)); rescaledir = 1; guard1 = TRUE; } else if (scale <= 1.3435752215134178E-138) { if ((g.re == 0.0) && (g.im == 0.0)) { *cs = 1.0; sn->re = 0.0; sn->im = 0.0; *r = f; } else { do { count++; fs_re *= 7.4428285367870146E+137; fs_im *= 7.4428285367870146E+137; gs_re *= 7.4428285367870146E+137; gs_im *= 7.4428285367870146E+137; scale *= 7.4428285367870146E+137; } while (!(scale > 1.3435752215134178E-138)); rescaledir = -1; guard1 = TRUE; } } else { guard1 = TRUE; } if (guard1 == TRUE) { scale = fs_re * fs_re + fs_im * fs_im; g2 = gs_re * gs_re + gs_im * gs_im; f2s = g2; if (1.0 > g2) { f2s = 1.0; } if (scale <= f2s * 2.0041683600089728E-292) { if ((f.re == 0.0) && (f.im == 0.0)) { *cs = 0.0; r->re = muDoubleScalarHypot(g.re, g.im); r->im = 0.0; f2s = muDoubleScalarHypot(gs_re, gs_im); sn->re = gs_re / f2s; sn->im = -gs_im / f2s; } else { emlrtPushRtStackR2012b(&vj_emlrtRSI, emlrtRootTLSGlobal); if (g2 < 0.0) { emlrtPushRtStackR2012b(&md_emlrtRSI, emlrtRootTLSGlobal); e_eml_error(); emlrtPopRtStackR2012b(&md_emlrtRSI, emlrtRootTLSGlobal); } g2s = muDoubleScalarSqrt(g2); emlrtPopRtStackR2012b(&vj_emlrtRSI, emlrtRootTLSGlobal); *cs = muDoubleScalarHypot(fs_re, fs_im) / g2s; f2s = muDoubleScalarAbs(f.re); g2 = muDoubleScalarAbs(f.im); if (g2 > f2s) { f2s = g2; } if (f2s > 1.0) { f2s = muDoubleScalarHypot(f.re, f.im); fs_re = f.re / f2s; fs_im = f.im / f2s; } else { g2 = 7.4428285367870146E+137 * f.re; scale = 7.4428285367870146E+137 * f.im; f2s = muDoubleScalarHypot(g2, scale); fs_re = g2 / f2s; fs_im = scale / f2s; } gs_re /= g2s; gs_im = -gs_im / g2s; sn->re = fs_re * gs_re - fs_im * gs_im; sn->im = fs_re * gs_im + fs_im * gs_re; r->re = *cs * f.re + (sn->re * g.re - sn->im * g.im); r->im = *cs * f.im + (sn->re * g.im + sn->im * g.re); } } else { emlrtPushRtStackR2012b(&wj_emlrtRSI, emlrtRootTLSGlobal); f2s = 1.0 + g2 / scale; if (f2s < 0.0) { emlrtPushRtStackR2012b(&md_emlrtRSI, emlrtRootTLSGlobal); e_eml_error(); emlrtPopRtStackR2012b(&md_emlrtRSI, emlrtRootTLSGlobal); } f2s = muDoubleScalarSqrt(f2s); emlrtPopRtStackR2012b(&wj_emlrtRSI, emlrtRootTLSGlobal); r->re = f2s * fs_re; r->im = f2s * fs_im; *cs = 1.0 / f2s; f2s = scale + g2; g2 = r->re / f2s; f2s = r->im / f2s; sn->re = g2 * gs_re - f2s * -gs_im; sn->im = g2 * -gs_im + f2s * gs_re; if (rescaledir > 0) { emlrtPushRtStackR2012b(&xj_emlrtRSI, emlrtRootTLSGlobal); if (1 > count) { b6 = FALSE; } else { b6 = (count > 2147483646); } if (b6) { emlrtPushRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); check_forloop_overflow_error(); emlrtPopRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); } emlrtPopRtStackR2012b(&xj_emlrtRSI, emlrtRootTLSGlobal); for (rescaledir = 1; rescaledir <= count; rescaledir++) { r->re *= 7.4428285367870146E+137; r->im *= 7.4428285367870146E+137; } } else { if (rescaledir < 0) { emlrtPushRtStackR2012b(&yj_emlrtRSI, emlrtRootTLSGlobal); if (1 > count) { b7 = FALSE; } else { b7 = (count > 2147483646); } if (b7) { emlrtPushRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); check_forloop_overflow_error(); emlrtPopRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); } emlrtPopRtStackR2012b(&yj_emlrtRSI, emlrtRootTLSGlobal); for (rescaledir = 1; rescaledir <= count; rescaledir++) { r->re *= 1.3435752215134178E-138; r->im *= 1.3435752215134178E-138; } } } } } } static void eml_matlab_ztgevc(const emxArray_creal_T *A, emxArray_creal_T *V) { emxArray_creal_T *work1; int32_T je; int32_T loop_ub; emxArray_creal_T *work2; emxArray_real_T *rworka; real_T SMALL; real_T BIG; real_T BIGNUM; real_T anorm; real_T y; real_T xmx; real_T ascale; int32_T b_je; real_T temp; real_T scale; real_T salpha_re; real_T salpha_im; real_T acoeff; boolean_T b23; boolean_T b24; real_T acoefa; int32_T jr; real_T dmin; real_T j; real_T d_re; real_T d_im; real_T work1_im; emlrtHeapReferenceStackEnterFcnR2012b(emlrtRootTLSGlobal); b_emxInit_creal_T(&work1, 1, &wd_emlrtRTEI, TRUE); je = work1->size[0]; work1->size[0] = A->size[0]; emxEnsureCapacity((emxArray__common *)work1, je, (int32_T)sizeof(creal_T), &bc_emlrtRTEI); loop_ub = A->size[0]; for (je = 0; je < loop_ub; je++) { work1->data[je].re = 0.0; work1->data[je].im = 0.0; } b_emxInit_creal_T(&work2, 1, &xd_emlrtRTEI, TRUE); je = work2->size[0]; work2->size[0] = A->size[0]; emxEnsureCapacity((emxArray__common *)work2, je, (int32_T)sizeof(creal_T), &bc_emlrtRTEI); loop_ub = A->size[0]; for (je = 0; je < loop_ub; je++) { work2->data[je].re = 0.0; work2->data[je].im = 0.0; } emxInit_real_T(&rworka, 1, &yd_emlrtRTEI, TRUE); SMALL = 2.2250738585072014E-308 * (real_T)A->size[0] / 2.2204460492503131E-16; BIG = 1.0 / SMALL; BIGNUM = 1.0 / (2.2250738585072014E-308 * (real_T)A->size[0]); je = rworka->size[0]; rworka->size[0] = A->size[0]; emxEnsureCapacity((emxArray__common *)rworka, je, (int32_T)sizeof(real_T), &bc_emlrtRTEI); loop_ub = A->size[0]; for (je = 0; je < loop_ub; je++) { rworka->data[je] = 0.0; } anorm = muDoubleScalarAbs(A->data[0].re) + muDoubleScalarAbs(A->data[0].im); for (loop_ub = 1; loop_ub - 1 <= A->size[0] - 2; loop_ub++) { for (je = 0; je < loop_ub; je++) { rworka->data[loop_ub] += muDoubleScalarAbs(A->data[je + A->size[0] * loop_ub].re) + muDoubleScalarAbs(A->data[je + A->size[0] * loop_ub].im); } y = rworka->data[loop_ub] + (muDoubleScalarAbs(A->data[loop_ub + A->size[0] * loop_ub].re) + muDoubleScalarAbs(A->data[loop_ub + A->size[0] * loop_ub]. im)); if (y > anorm) { anorm = y; } } xmx = anorm; if (2.2250738585072014E-308 > anorm) { xmx = 2.2250738585072014E-308; } ascale = 1.0 / xmx; emlrtForLoopVectorCheckR2012b(A->size[0], -1.0, 1.0, mxDOUBLE_CLASS, A->size[0], &ie_emlrtRTEI, emlrtRootTLSGlobal); for (je = 1; je - 1 < A->size[0]; je++) { b_je = A->size[0] - je; y = (muDoubleScalarAbs(A->data[b_je + A->size[0] * b_je].re) + muDoubleScalarAbs(A->data[b_je + A->size[0] * b_je].im)) * ascale; if (1.0 > y) { y = 1.0; } temp = 1.0 / y; xmx = temp * A->data[b_je + A->size[0] * b_je].re; scale = temp * A->data[b_je + A->size[0] * b_je].im; salpha_re = ascale * xmx; salpha_im = ascale * scale; acoeff = temp * ascale; if ((muDoubleScalarAbs(temp) >= 2.2250738585072014E-308) && (muDoubleScalarAbs(acoeff) < SMALL)) { b23 = TRUE; } else { b23 = FALSE; } if ((muDoubleScalarAbs(salpha_re) + muDoubleScalarAbs(salpha_im) >= 2.2250738585072014E-308) && (muDoubleScalarAbs(salpha_re) + muDoubleScalarAbs(salpha_im) < SMALL)) { b24 = TRUE; } else { b24 = FALSE; } scale = 1.0; if (b23) { xmx = anorm; if (BIG < anorm) { xmx = BIG; } scale = SMALL / muDoubleScalarAbs(temp) * xmx; } if (b24) { xmx = 1.0; if (BIG < 1.0) { xmx = BIG; } y = SMALL / (muDoubleScalarAbs(salpha_re) + muDoubleScalarAbs(salpha_im)) * xmx; if (y > scale) { scale = y; } } if (b23 || b24) { y = muDoubleScalarAbs(acoeff); xmx = muDoubleScalarAbs(salpha_re) + muDoubleScalarAbs(salpha_im); if (1.0 > y) { y = 1.0; } if (xmx > y) { y = xmx; } y = 1.0 / (2.2250738585072014E-308 * y); if (y < scale) { scale = y; } if (b23) { acoeff = ascale * (scale * temp); } else { acoeff *= scale; } if (b24) { salpha_re *= scale; salpha_im *= scale; } else { salpha_re *= scale; salpha_im *= scale; } } acoefa = muDoubleScalarAbs(acoeff); for (jr = 0; jr < A->size[0]; jr++) { work1->data[jr].re = 0.0; work1->data[jr].im = 0.0; } work1->data[b_je].re = 1.0; work1->data[b_je].im = 0.0; dmin = 2.2204460492503131E-16 * acoefa * anorm; y = 2.2204460492503131E-16 * (muDoubleScalarAbs(salpha_re) + muDoubleScalarAbs(salpha_im)); if (y > dmin) { dmin = y; } if (2.2250738585072014E-308 > dmin) { dmin = 2.2250738585072014E-308; } for (jr = 0; jr < b_je; jr++) { work1->data[jr].re = acoeff * A->data[jr + A->size[0] * b_je].re; work1->data[jr].im = acoeff * A->data[jr + A->size[0] * b_je].im; } work1->data[b_je].re = 1.0; work1->data[b_je].im = 0.0; emlrtForLoopVectorCheckR2012b((real_T)(b_je + 1) - 1.0, -1.0, 1.0, mxDOUBLE_CLASS, b_je, &he_emlrtRTEI, emlrtRootTLSGlobal); for (loop_ub = 0; loop_ub < b_je; loop_ub++) { j = ((real_T)(b_je + 1) - 1.0) + -(real_T)loop_ub; d_re = acoeff * A->data[((int32_T)j + A->size[0] * ((int32_T)j - 1)) - 1]. re - salpha_re; d_im = acoeff * A->data[((int32_T)j + A->size[0] * ((int32_T)j - 1)) - 1]. im - salpha_im; if (muDoubleScalarAbs(d_re) + muDoubleScalarAbs(d_im) <= dmin) { d_re = dmin; d_im = 0.0; } if ((muDoubleScalarAbs(d_re) + muDoubleScalarAbs(d_im) < 1.0) && (muDoubleScalarAbs(work1->data[(int32_T)j - 1].re) + muDoubleScalarAbs (work1->data[(int32_T)j - 1].im) >= BIGNUM * (muDoubleScalarAbs(d_re) + muDoubleScalarAbs(d_im)))) { temp = 1.0 / (muDoubleScalarAbs(work1->data[(int32_T)j - 1].re) + muDoubleScalarAbs(work1->data[(int32_T)j - 1].im)); for (jr = 0; jr <= b_je; jr++) { work1->data[jr].re *= temp; work1->data[jr].im *= temp; } } temp = -work1->data[(int32_T)j - 1].re; work1_im = -work1->data[(int32_T)j - 1].im; if (d_im == 0.0) { if (work1_im == 0.0) { work1->data[(int32_T)j - 1].re = temp / d_re; work1->data[(int32_T)j - 1].im = 0.0; } else if (temp == 0.0) { work1->data[(int32_T)j - 1].re = 0.0; work1->data[(int32_T)j - 1].im = work1_im / d_re; } else { work1->data[(int32_T)j - 1].re = temp / d_re; work1->data[(int32_T)j - 1].im = work1_im / d_re; } } else if (d_re == 0.0) { if (temp == 0.0) { work1->data[(int32_T)j - 1].re = work1_im / d_im; work1->data[(int32_T)j - 1].im = 0.0; } else if (work1_im == 0.0) { work1->data[(int32_T)j - 1].re = 0.0; work1->data[(int32_T)j - 1].im = -(temp / d_im); } else { work1->data[(int32_T)j - 1].re = work1_im / d_im; work1->data[(int32_T)j - 1].im = -(temp / d_im); } } else { y = muDoubleScalarAbs(d_re); xmx = muDoubleScalarAbs(d_im); if (y > xmx) { scale = d_im / d_re; xmx = d_re + scale * d_im; work1->data[(int32_T)j - 1].re = (temp + scale * work1_im) / xmx; work1->data[(int32_T)j - 1].im = (work1_im - scale * temp) / xmx; } else if (xmx == y) { if (d_re > 0.0) { scale = 0.5; } else { scale = -0.5; } if (d_im > 0.0) { xmx = 0.5; } else { xmx = -0.5; } work1->data[(int32_T)j - 1].re = (temp * scale + work1_im * xmx) / y; work1->data[(int32_T)j - 1].im = (work1_im * scale - temp * xmx) / y; } else { scale = d_re / d_im; xmx = d_im + scale * d_re; work1->data[(int32_T)j - 1].re = (scale * temp + work1_im) / xmx; work1->data[(int32_T)j - 1].im = (scale * work1_im - temp) / xmx; } } if (j > 1.0) { if (muDoubleScalarAbs(work1->data[(int32_T)j - 1].re) + muDoubleScalarAbs(work1->data[(int32_T)j - 1].im) > 1.0) { temp = 1.0 / (muDoubleScalarAbs(work1->data[(int32_T)j - 1].re) + muDoubleScalarAbs(work1->data[(int32_T)j - 1].im)); if (acoefa * rworka->data[(int32_T)j - 1] >= BIGNUM * temp) { for (jr = 0; jr <= b_je; jr++) { work1->data[jr].re *= temp; work1->data[jr].im *= temp; } } } d_re = acoeff * work1->data[(int32_T)j - 1].re; d_im = acoeff * work1->data[(int32_T)j - 1].im; for (jr = 0; jr <= (int32_T)j - 2; jr++) { xmx = d_re * A->data[jr + A->size[0] * ((int32_T)j - 1)].re - d_im * A->data[jr + A->size[0] * ((int32_T)j - 1)].im; scale = d_re * A->data[jr + A->size[0] * ((int32_T)j - 1)].im + d_im * A->data[jr + A->size[0] * ((int32_T)j - 1)].re; work1->data[jr].re += xmx; work1->data[jr].im += scale; } } } for (jr = 0; jr < A->size[0]; jr++) { work2->data[jr].re = 0.0; work2->data[jr].im = 0.0; } for (loop_ub = 0; loop_ub <= b_je; loop_ub++) { for (jr = 0; jr < A->size[0]; jr++) { xmx = V->data[jr + V->size[0] * loop_ub].re * work1->data[loop_ub].re - V->data[jr + V->size[0] * loop_ub].im * work1->data[loop_ub].im; scale = V->data[jr + V->size[0] * loop_ub].re * work1->data[loop_ub].im + V->data[jr + V->size[0] * loop_ub].im * work1->data[loop_ub].re; work2->data[jr].re += xmx; work2->data[jr].im += scale; } } xmx = muDoubleScalarAbs(work2->data[0].re) + muDoubleScalarAbs(work2->data[0] .im); if (A->size[0] > 1) { for (jr = 1; jr - 1 <= A->size[0] - 2; jr++) { y = muDoubleScalarAbs(work2->data[jr].re) + muDoubleScalarAbs (work2->data[jr].im); if (y > xmx) { xmx = y; } } } if (xmx > 2.2250738585072014E-308) { temp = 1.0 / xmx; for (jr = 0; jr < A->size[0]; jr++) { V->data[jr + V->size[0] * b_je].re = temp * work2->data[jr].re; V->data[jr + V->size[0] * b_je].im = temp * work2->data[jr].im; } } else { for (jr = 0; jr < A->size[0]; jr++) { V->data[jr + V->size[0] * b_je].re = 0.0; V->data[jr + V->size[0] * b_je].im = 0.0; } } } emxFree_real_T(&rworka); emxFree_creal_T(&work2); emxFree_creal_T(&work1); emlrtHeapReferenceStackLeaveFcnR2012b(emlrtRootTLSGlobal); } void eig(const emxArray_real_T *A, emxArray_creal_T *V, emxArray_creal_T *D) { emxArray_creal_T *b_A; int32_T j; int32_T loop_ub; emxArray_creal_T *c_A; emxArray_creal_T *alpha1; emxArray_creal_T *beta1; real_T info; int32_T coltop; ptrdiff_t n_t; ptrdiff_t incx_t; double * xix0_t; real_T colnorm; int32_T b; boolean_T b_coltop; creal_T b_V; real_T alpha1_re; real_T alpha1_im; real_T beta1_re; real_T beta1_im; real_T brm; real_T s; boolean_T overflow; emlrtHeapReferenceStackEnterFcnR2012b(emlrtRootTLSGlobal); emxInit_creal_T(&b_A, 2, &wb_emlrtRTEI, TRUE); emlrtPushRtStackR2012b(&oi_emlrtRSI, emlrtRootTLSGlobal); j = b_A->size[0] * b_A->size[1]; b_A->size[0] = A->size[0]; b_A->size[1] = A->size[1]; emxEnsureCapacity((emxArray__common *)b_A, j, (int32_T)sizeof(creal_T), &wb_emlrtRTEI); loop_ub = A->size[0] * A->size[1]; for (j = 0; j < loop_ub; j++) { b_A->data[j].re = A->data[j]; b_A->data[j].im = 0.0; } emxInit_creal_T(&c_A, 2, &wb_emlrtRTEI, TRUE); emlrtPushRtStackR2012b(&si_emlrtRSI, emlrtRootTLSGlobal); emlrtPushRtStackR2012b(&ti_emlrtRSI, emlrtRootTLSGlobal); j = c_A->size[0] * c_A->size[1]; c_A->size[0] = b_A->size[0]; c_A->size[1] = b_A->size[1]; emxEnsureCapacity((emxArray__common *)c_A, j, (int32_T)sizeof(creal_T), &wb_emlrtRTEI); loop_ub = b_A->size[0] * b_A->size[1]; for (j = 0; j < loop_ub; j++) { c_A->data[j] = b_A->data[j]; } b_emxInit_creal_T(&alpha1, 1, &wb_emlrtRTEI, TRUE); b_emxInit_creal_T(&beta1, 1, &wb_emlrtRTEI, TRUE); eml_matlab_zggev(c_A, &info, alpha1, beta1, V); emlrtPopRtStackR2012b(&ti_emlrtRSI, emlrtRootTLSGlobal); emxFree_creal_T(&c_A); if (b_A->size[0] > 0) { loop_ub = (b_A->size[0] - 1) * b_A->size[0] + 1; emlrtPushRtStackR2012b(&ui_emlrtRSI, emlrtRootTLSGlobal); if (loop_ub > MAX_int32_T - b_A->size[0]) { emlrtPushRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); check_forloop_overflow_error(); emlrtPopRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); } emlrtPopRtStackR2012b(&ui_emlrtRSI, emlrtRootTLSGlobal); for (coltop = 0; coltop + 1 <= loop_ub; coltop += b_A->size[0]) { emlrtPushRtStackR2012b(&vi_emlrtRSI, emlrtRootTLSGlobal); emlrtPushRtStackR2012b(&sk_emlrtRSI, emlrtRootTLSGlobal); emlrtPushRtStackR2012b(&tk_emlrtRSI, emlrtRootTLSGlobal); emlrtPushRtStackR2012b(&wk_emlrtRSI, emlrtRootTLSGlobal); emlrtPushRtStackR2012b(&af_emlrtRSI, emlrtRootTLSGlobal); n_t = (ptrdiff_t)(b_A->size[0]); emlrtPopRtStackR2012b(&af_emlrtRSI, emlrtRootTLSGlobal); emlrtPopRtStackR2012b(&wk_emlrtRSI, emlrtRootTLSGlobal); emlrtPushRtStackR2012b(&xk_emlrtRSI, emlrtRootTLSGlobal); emlrtPushRtStackR2012b(&af_emlrtRSI, emlrtRootTLSGlobal); incx_t = (ptrdiff_t)(1); emlrtPopRtStackR2012b(&af_emlrtRSI, emlrtRootTLSGlobal); emlrtPopRtStackR2012b(&xk_emlrtRSI, emlrtRootTLSGlobal); emlrtPushRtStackR2012b(&yk_emlrtRSI, emlrtRootTLSGlobal); xix0_t = (double *)(&V->data[coltop]); emlrtPopRtStackR2012b(&yk_emlrtRSI, emlrtRootTLSGlobal); colnorm = dznrm2(&n_t, xix0_t, &incx_t); emlrtPopRtStackR2012b(&tk_emlrtRSI, emlrtRootTLSGlobal); emlrtPopRtStackR2012b(&sk_emlrtRSI, emlrtRootTLSGlobal); emlrtPopRtStackR2012b(&vi_emlrtRSI, emlrtRootTLSGlobal); b = coltop + b_A->size[0]; emlrtPushRtStackR2012b(&wi_emlrtRSI, emlrtRootTLSGlobal); if (coltop + 1 > b) { b_coltop = FALSE; } else { b_coltop = (b > 2147483646); } if (b_coltop) { emlrtPushRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); check_forloop_overflow_error(); emlrtPopRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); } emlrtPopRtStackR2012b(&wi_emlrtRSI, emlrtRootTLSGlobal); for (j = coltop; j + 1 <= b; j++) { b_V = V->data[j]; V->data[j] = eml_div(b_V, colnorm); } } } emxFree_creal_T(&b_A); emlrtPopRtStackR2012b(&si_emlrtRSI, emlrtRootTLSGlobal); emlrtPopRtStackR2012b(&oi_emlrtRSI, emlrtRootTLSGlobal); emlrtPushRtStackR2012b(&pi_emlrtRSI, emlrtRootTLSGlobal); j = alpha1->size[0]; emxEnsureCapacity((emxArray__common *)alpha1, j, (int32_T)sizeof(creal_T), &wb_emlrtRTEI); loop_ub = alpha1->size[0]; for (j = 0; j < loop_ub; j++) { alpha1_re = alpha1->data[j].re; alpha1_im = alpha1->data[j].im; beta1_re = beta1->data[j].re; beta1_im = beta1->data[j].im; if (beta1_im == 0.0) { if (alpha1_im == 0.0) { alpha1->data[j].re = alpha1_re / beta1_re; alpha1->data[j].im = 0.0; } else if (alpha1_re == 0.0) { alpha1->data[j].re = 0.0; alpha1->data[j].im = alpha1_im / beta1_re; } else { alpha1->data[j].re = alpha1_re / beta1_re; alpha1->data[j].im = alpha1_im / beta1_re; } } else if (beta1_re == 0.0) { if (alpha1_re == 0.0) { alpha1->data[j].re = alpha1_im / beta1_im; alpha1->data[j].im = 0.0; } else if (alpha1_im == 0.0) { alpha1->data[j].re = 0.0; alpha1->data[j].im = -(alpha1_re / beta1_im); } else { alpha1->data[j].re = alpha1_im / beta1_im; alpha1->data[j].im = -(alpha1_re / beta1_im); } } else { brm = muDoubleScalarAbs(beta1_re); colnorm = muDoubleScalarAbs(beta1_im); if (brm > colnorm) { s = beta1_im / beta1_re; colnorm = beta1_re + s * beta1_im; alpha1->data[j].re = (alpha1_re + s * alpha1_im) / colnorm; alpha1->data[j].im = (alpha1_im - s * alpha1_re) / colnorm; } else if (colnorm == brm) { if (beta1_re > 0.0) { s = 0.5; } else { s = -0.5; } if (beta1_im > 0.0) { colnorm = 0.5; } else { colnorm = -0.5; } alpha1->data[j].re = (alpha1_re * s + alpha1_im * colnorm) / brm; alpha1->data[j].im = (alpha1_im * s - alpha1_re * colnorm) / brm; } else { s = beta1_re / beta1_im; colnorm = beta1_im + s * beta1_re; alpha1->data[j].re = (s * alpha1_re + alpha1_im) / colnorm; alpha1->data[j].im = (s * alpha1_im - alpha1_re) / colnorm; } } } emxFree_creal_T(&beta1); loop_ub = alpha1->size[0]; coltop = alpha1->size[0]; j = D->size[0] * D->size[1]; D->size[0] = loop_ub; emxEnsureCapacity((emxArray__common *)D, j, (int32_T)sizeof(creal_T), &wb_emlrtRTEI); j = D->size[0] * D->size[1]; D->size[1] = coltop; emxEnsureCapacity((emxArray__common *)D, j, (int32_T)sizeof(creal_T), &wb_emlrtRTEI); loop_ub *= coltop; for (j = 0; j < loop_ub; j++) { D->data[j].re = 0.0; D->data[j].im = 0.0; } emlrtPushRtStackR2012b(&cl_emlrtRSI, emlrtRootTLSGlobal); if (1 > alpha1->size[0]) { overflow = FALSE; } else { overflow = (alpha1->size[0] > 2147483646); } if (overflow) { emlrtPushRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); check_forloop_overflow_error(); emlrtPopRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal); } emlrtPopRtStackR2012b(&cl_emlrtRSI, emlrtRootTLSGlobal); for (j = 0; j + 1 <= alpha1->size[0]; j++) { D->data[j + D->size[0] * j] = alpha1->data[j]; } emxFree_creal_T(&alpha1); emlrtPopRtStackR2012b(&pi_emlrtRSI, emlrtRootTLSGlobal); if (info < 0.0) { emlrtPushRtStackR2012b(&qi_emlrtRSI, emlrtRootTLSGlobal); c_eml_warning(); emlrtPopRtStackR2012b(&qi_emlrtRSI, emlrtRootTLSGlobal); } else { if (info > 0.0) { emlrtPushRtStackR2012b(&ri_emlrtRSI, emlrtRootTLSGlobal); d_eml_warning(); emlrtPopRtStackR2012b(&ri_emlrtRSI, emlrtRootTLSGlobal); } } emlrtHeapReferenceStackLeaveFcnR2012b(emlrtRootTLSGlobal); } /* End of code generation (eig.c) */