1 #======================================================================

     2 #                    L I B Q Z . P L 

     3 #                    doc: Thu Mar 12 15:23:15 2015

     4 #                    dlm: Thu Mar 12 20:52:43 2015

     5 #                    (c) 2015 A.M. Thurnherr

     6 #                    uE-Info: 36 0 NIL 0 0 72 2 2 4 NIL ofnI

     7 #======================================================================

     8 

     9 # adaptation of EISPACK routines

    10 

    11 # www.netlib.org/eispack

    12 

    13 sub eig($$)

    14 {

    15 	my($aR,$bR) = @_;							# args passed as refs

    16 

    17 	my($N) = scalar(@{aR});

    18 	croak("eig(A,B): A & B must be matching square matrices\n")

    19 		unless (@{bR} == $N) && (@{$aR->[0]} == $N) && (@{$bR->[0]} == $N);

    20 

    21 		ar = new double[n];

    22 		ai = new double[n];

    23 		beta = new double[n];

    24 		var A = (double[,])a.Clone();

    25 		var B = (double[,])b.Clone();

    26 

    27 	my($matZ) = 1;

    28 	my($iErr) = 0;

    29 

    30 	my(@Z);

    31 	QZhes($aR,$bR,\@Z);								# reduce A/B to upper Hessenberg/triangular forms

    32 	QZit($aR,$bR,\@Z,\$iErr);						# reduce Hess A to quasi-triangular form

    33 	QZval($aR,$bR,\@Z,\@alphaR,\@alphaI,\@beta);	# reduce A further

    34 	QZvec($aR,$bR,\@Z,\@alphaR,\@alphaI,\@beta);	# compute eigenvectors & eigenvalues

    35 }

    36 

    37 

    38 	/// <summary>Returns the real parts of the alpha values.</summary>

    39 	public double[] RealAlphas

    40 	{

    41 		get { return ar; }

    42 	}

    43 

    44 	/// <summary>Returns the imaginary parts of the alpha values.</summary>

    45 	public double[] ImaginaryAlphas

    46 	{

    47 		get { return ai; }

    48 	}

    49 

    50 	/// <summary>Returns the beta values.</summary>

    51 	public double[] Betas

    52 	{

    53 		get { return beta; }

    54 	}

    55 

    56 	/// <summary>

    57 	///   Returns true if matrix B is singular.

    58 	/// </summary>

    59 	/// <remarks>

    60 	///   This method checks if any of the generated betas is zero. It

    61 	///   does not says that the problem is singular, but only that one

    62 	///   of the matrices of the pencil (A,B) is singular.

    63 	/// </remarks>

    64 	public bool IsSingular

    65 	{

    66 		get

    67 		{

    68 			for (int i = 0; i < n; i++)

    69 				if (beta[i] == 0)

    70 					return true;

    71 			return false;

    72 		}

    73 	}

    74 

    75 	/// <summary>

    76 	///   Returns true if the eigenvalue problem is degenerate (ill-posed).

    77 	/// </summary>

    78 	public bool IsDegenerate

    79 	{

    80 		get

    81 		{

    82 			for (int i = 0; i < n; i++)

    83 				if (beta[i] == 0 && ar[i] == 0)

    84 					return true;

    85 			return false;

    86 		}

    87 	}

    88 

    89 	/// <summary>Returns the real parts of the eigenvalues.</summary>

    90 	/// <remarks>

    91 	///   The eigenvalues are computed using the ratio alpha[i]/beta[i],

    92 	///   which can lead to valid, but infinite eigenvalues.

    93 	/// </remarks>

    94 	public double[] RealEigenvalues

    95 	{

    96 		get

    97 		{

    98 			// ((alfr+i*alfi)/beta)

    99 			double[] eval = new double[n];

   100 			for (int i = 0; i < n; i++)

   101 				eval[i] = ar[i] / beta[i];

   102 			return eval;

   103 		}

   104 	}

   105 

   106 	/// <summary>Returns the imaginary parts of the eigenvalues.</summary>  

   107 	/// <remarks>

   108 	///   The eigenvalues are computed using the ratio alpha[i]/beta[i],

   109 	///   which can lead to valid, but infinite eigenvalues.

   110 	/// </remarks>

   111 	public double[] ImaginaryEigenvalues

   112 	{

   113 		get

   114 		{

   115 			// ((alfr+i*alfi)/beta)

   116 			double[] eval = new double[n];

   117 			for (int i = 0; i < n; i++)

   118 				eval[i] = ai[i] / beta[i];

   119 			return eval;

   120 		}

   121 	}

   122 

   123 	/// <summary>Returns the eigenvector matrix.</summary>

   124 	public double[,] Eigenvectors

   125 	{

   126 		get

   127 		{

   128 			return Z;

   129 		}

   130 	}

   131 

   132 	/// <summary>Returns the block diagonal eigenvalue matrix.</summary>

   133 	public double[,] DiagonalMatrix

   134 	{

   135 		get

   136 		{

   137 			double[,] x = new double[n, n];

   138 

   139 			for (int i = 0; i < n; i++)

   140 			{

   141 				for (int j = 0; j < n; j++)

   142 					x[i, j] = 0.0;

   143 

   144 				x[i, i] = ar[i] / beta[i];

   145 				if (ai[i] > 0)

   146 					x[i, i + 1] = ai[i] / beta[i];

   147 				else if (ai[i] < 0)

   148 					x[i, i - 1] = ai[i] / beta[i];

   149 			}

   150 

   151 			return x;

   152 		}

   153 	}

   154 

   155 

   156 

   157 	#region EISPACK Routines

   158 	/// <summary>

   159 	///   Adaptation of the original Fortran QZHES routine from EISPACK.

   160 	/// </summary>

   161 	/// <remarks>

   162 	///   This subroutine is the first step of the qz algorithm

   163 	///   for solving generalized matrix eigenvalue problems,

   164 	///   siam j. numer. anal. 10, 241-256(1973) by moler and stewart.

   165 	///

   166 	///   This subroutine accepts a pair of real general matrices and

   167 	///   reduces one of them to upper hessenberg form and the other

   168 	///   to upper triangular form using orthogonal transformations.

   169 	///   it is usually followed by  qzit,	qzval  and, possibly,  qzvec.

   170 	///  

   171 	///   For the full documentation, please check the original function.

   172 	/// </remarks>

   173 	private static int qzhes(int n, double[,] a, double[,] b, bool matz, double[,] z)

   174 	{

   175 		int i, j, k, l;

   176 		double r, s, t;

   177 		int l1;

   178 		double u1, u2, v1, v2;

   179 		int lb, nk1;

   180 		double rho;

   181 

   182 

   183 		if (matz)

   184 		{

   185 			// If we are interested in computing the

   186 			//	eigenvectors, set Z to identity(n,n)

   187 			for (j = 0; j < n; ++j)

   188 			{

   189 				for (i = 0; i < n; ++i)

   190 					z[i, j] = 0.0;

   191 				z[j, j] = 1.0;

   192 			}

   193 		}

   194 

   195 		// Reduce b to upper triangular form

   196 		if (n <= 1) return 0;

   197 		for (l = 0; l < n - 1; ++l)

   198 		{

   199 			l1 = l + 1;

   200 			s = 0.0;

   201 

   202 			for (i = l1; i < n; ++i)

   203 				s += (System.Math.Abs(b[i, l]));

   204 

   205 			if (s == 0.0) continue;

   206 			s += (System.Math.Abs(b[l, l]));

   207 			r = 0.0;

   208 

   209 			for (i = l; i < n; ++i)

   210 			{

   211 				// Computing 2nd power

   212 				b[i, l] /= s;

   213 				r += b[i, l] * b[i, l];

   214 			}

   215 

   216 			r = Special.Sign(System.Math.Sqrt(r), b[l, l]);

   217 			b[l, l] += r;

   218 			rho = r * b[l, l];

   219 

   220 			for (j = l1; j < n; ++j)

   221 			{

   222 				t = 0.0;

   223 				for (i = l; i < n; ++i)

   224 					t += b[i, l] * b[i, j];

   225 				t = -t / rho;

   226 				for (i = l; i < n; ++i)

   227 					b[i, j] += t * b[i, l];

   228 			}

   229 

   230 			for (j = 0; j < n; ++j)

   231 			{

   232 				t = 0.0;

   233 				for (i = l; i < n; ++i)

   234 					t += b[i, l] * a[i, j];

   235 				t = -t / rho;

   236 				for (i = l; i < n; ++i)

   237 					a[i, j] += t * b[i, l];

   238 			}

   239 

   240 			b[l, l] = -s * r;

   241 			for (i = l1; i < n; ++i)

   242 				b[i, l] = 0.0;

   243 		}

   244 

   245 		// Reduce a to upper hessenberg form, while keeping b triangular

   246 		if (n == 2) return 0;

   247 		for (k = 0; k < n - 2; ++k)

   248 		{

   249 			nk1 = n - 2 - k;

   250 

   251 			// for l=n-1 step -1 until k+1 do

   252 			for (lb = 0; lb < nk1; ++lb)

   253 			{

   254 				l = n - lb - 2;

   255 				l1 = l + 1;

   256 

   257 				// Zero a(l+1,k)

   258 				s = (System.Math.Abs(a[l, k])) + (System.Math.Abs(a[l1, k]));

   259 

   260 				if (s == 0.0) continue;

   261 				u1 = a[l, k] / s;

   262 				u2 = a[l1, k] / s;

   263 				r = Special.Sign(Math.Sqrt(u1 * u1 + u2 * u2), u1);

   264 				v1 = -(u1 + r) / r;

   265 				v2 = -u2 / r;

   266 				u2 = v2 / v1;

   267 

   268 				for (j = k; j < n; ++j)

   269 				{

   270 					t = a[l, j] + u2 * a[l1, j];

   271 					a[l, j] += t * v1;

   272 					a[l1, j] += t * v2;

   273 				}

   274 

   275 				a[l1, k] = 0.0;

   276 

   277 				for (j = l; j < n; ++j)

   278 				{

   279 					t = b[l, j] + u2 * b[l1, j];

   280 					b[l, j] += t * v1;

   281 					b[l1, j] += t * v2;

   282 				}

   283 

   284 				// Zero b(l+1,l)

   285 				s = (System.Math.Abs(b[l1, l1])) + (System.Math.Abs(b[l1, l]));

   286 

   287 				if (s == 0.0) continue;

   288 				u1 = b[l1, l1] / s;

   289 				u2 = b[l1, l] / s;

   290 				r = Special.Sign(Math.Sqrt(u1 * u1 + u2 * u2), u1);

   291 				v1 = -(u1 + r) / r;

   292 				v2 = -u2 / r;

   293 				u2 = v2 / v1;

   294 

   295 				for (i = 0; i <= l1; ++i)

   296 				{

   297 					t = b[i, l1] + u2 * b[i, l];

   298 					b[i, l1] += t * v1;

   299 					b[i, l] += t * v2;

   300 				}

   301 

   302 				b[l1, l] = 0.0;

   303 

   304 				for (i = 0; i < n; ++i)

   305 				{

   306 					t = a[i, l1] + u2 * a[i, l];

   307 					a[i, l1] += t * v1;

   308 					a[i, l] += t * v2;

   309 				}

   310 

   311 				if (matz)

   312 				{

   313 					for (i = 0; i < n; ++i)

   314 					{

   315 						t = z[i, l1] + u2 * z[i, l];

   316 						z[i, l1] += t * v1;

   317 						z[i, l] += t * v2;

   318 					}

   319 				}

   320 			}

   321 		}

   322 

   323 		return 0;

   324 	}

   325 

   326 	/// <summary>

   327 	///   Adaptation of the original Fortran QZIT routine from EISPACK.

   328 	/// </summary>

   329 	/// <remarks>

   330 	///   This subroutine is the second step of the qz algorithm

   331 	///   for solving generalized matrix eigenvalue problems,

   332 	///   siam j. numer. anal. 10, 241-256(1973) by moler and stewart,

   333 	///   as modified in technical note nasa tn d-7305(1973) by ward.

   334 	///  

   335 	///   This subroutine accepts a pair of real matrices, one of them

   336 	///   in upper hessenberg form and the other in upper triangular form.

   337 	///   it reduces the hessenberg matrix to quasi-triangular form using

   338 	///   orthogonal transformations while maintaining the triangular form

   339 	///   of the other matrix.	it is usually preceded by  qzhes  and

   340 	///   followed by  qzval  and, possibly,  qzvec.

   341 	///  

   342 	///   For the full documentation, please check the original function.

   343 	/// </remarks>

   344 	private static int qzit(int n, double[,] a, double[,] b, double eps1, bool matz, double[,] z, ref int ierr)

   345 	{

   346 

   347 		int i, j, k, l = 0;

   348 		double r, s, t, a1, a2, a3 = 0;

   349 		int k1, k2, l1, ll;

   350 		double u1, u2, u3;

   351 		double v1, v2, v3;

   352 		double a11, a12, a21, a22, a33, a34, a43, a44;

   353 		double b11, b12, b22, b33, b34, b44;

   354 		int na, en, ld;

   355 		double ep;

   356 		double sh = 0;

   357 		int km1, lm1 = 0;

   358 		double ani, bni;

   359 		int ish, itn, its, enm2, lor1;

   360 		double epsa, epsb, anorm = 0, bnorm = 0;

   361 		int enorn;

   362 		bool notlas;

   363 

   364 

   365 		ierr = 0;

   366 

   367 		#region Compute epsa and epsb

   368 		for (i = 0; i < n; ++i)

   369 		{

   370 			ani = 0.0;

   371 			bni = 0.0;

   372 

   373 			if (i != 0)

   374 				ani = (Math.Abs(a[i, (i - 1)]));

   375 

   376 			for (j = i; j < n; ++j)

   377 			{

   378 				ani += Math.Abs(a[i, j]);

   379 				bni += Math.Abs(b[i, j]);

   380 			}

   381 

   382 			if (ani > anorm) anorm = ani;

   383 			if (bni > bnorm) bnorm = bni;

   384 		}

   385 

   386 		if (anorm == 0.0) anorm = 1.0;

   387 		if (bnorm == 0.0) bnorm = 1.0;

   388 

   389 		ep = eps1;

   390 		if (ep == 0.0)

   391 		{

   392 			// Use roundoff level if eps1 is zero

   393 			ep = Special.Epslon(1.0);

   394 		}

   395 

   396 		epsa = ep * anorm;

   397 		epsb = ep * bnorm;

   398 		#endregion

   399 

   400 

   401 		// Reduce a to quasi-triangular form, while keeping b triangular

   402 		lor1 = 0;

   403 		enorn = n;

   404 		en = n - 1;

   405 		itn = n * 30;

   406 

   407 	// Begin QZ step

   408 	L60:

   409 		if (en <= 1) goto L1001;

   410 		if (!matz) enorn = en + 1;

   411 

   412 		its = 0;

   413 		na = en - 1;

   414 		enm2 = na;

   415 

   416 	L70:

   417 		ish = 2;

   418 		// Check for convergence or reducibility.

   419 		for (ll = 0; ll <= en; ++ll)

   420 		{

   421 			lm1 = en - ll - 1;

   422 			l = lm1 + 1;

   423 

   424 			if (l + 1 == 1)

   425 				goto L95;

   426 

   427 			if ((Math.Abs(a[l, lm1])) <= epsa)

   428 				break;

   429 		}

   430 

   431 	L90:

   432 		a[l, lm1] = 0.0;

   433 		if (l < na) goto L95;

   434 

   435 		// 1-by-1 or 2-by-2 block isolated

   436 		en = lm1;

   437 		goto L60;

   438 

   439 	// Check for small top of b

   440 	L95:

   441 		ld = l;

   442 

   443 	L100:

   444 		l1 = l + 1;

   445 		b11 = b[l, l];

   446 

   447 		if (Math.Abs(b11) > epsb) goto L120;

   448 

   449 		b[l, l] = 0.0;

   450 		s = (Math.Abs(a[l, l]) + Math.Abs(a[l1, l]));

   451 		u1 = a[l, l] / s;

   452 		u2 = a[l1, l] / s;

   453 		r = Special.Sign(Math.Sqrt(u1 * u1 + u2 * u2), u1);

   454 		v1 = -(u1 + r) / r;

   455 		v2 = -u2 / r;

   456 		u2 = v2 / v1;

   457 

   458 		for (j = l; j < enorn; ++j)

   459 		{

   460 			t = a[l, j] + u2 * a[l1, j];

   461 			a[l, j] += t * v1;

   462 			a[l1, j] += t * v2;

   463 

   464 			t = b[l, j] + u2 * b[l1, j];

   465 			b[l, j] += t * v1;

   466 			b[l1, j] += t * v2;

   467 		}

   468 

   469 		if (l != 0)

   470 			a[l, lm1] = -a[l, lm1];

   471 

   472 		lm1 = l;

   473 		l = l1;

   474 		goto L90;

   475 

   476 	L120:

   477 		a11 = a[l, l] / b11;

   478 		a21 = a[l1, l] / b11;

   479 		if (ish == 1) goto L140;

   480 

   481 		// Iteration strategy

   482 		if (itn == 0) goto L1000;

   483 		if (its == 10) goto L155;

   484 

   485 		// Determine type of shift

   486 		b22 = b[l1, l1];

   487 		if (Math.Abs(b22) < epsb) b22 = epsb;

   488 		b33 = b[na, na];

   489 		if (Math.Abs(b33) < epsb) b33 = epsb;

   490 		b44 = b[en, en];

   491 		if (Math.Abs(b44) < epsb) b44 = epsb;

   492 		a33 = a[na, na] / b33;

   493 		a34 = a[na, en] / b44;

   494 		a43 = a[en, na] / b33;

   495 		a44 = a[en, en] / b44;

   496 		b34 = b[na, en] / b44;

   497 		t = (a43 * b34 - a33 - a44) * .5;

   498 		r = t * t + a34 * a43 - a33 * a44;

   499 		if (r < 0.0) goto L150;

   500 

   501 		// Determine single shift zeroth column of a

   502 		ish = 1;

   503 		r = Math.Sqrt(r);

   504 		sh = -t + r;

   505 		s = -t - r;

   506 		if (Math.Abs(s - a44) < Math.Abs(sh - a44))

   507 			sh = s;

   508 

   509 		// Look for two consecutive small sub-diagonal elements of a.

   510 		for (ll = ld; ll + 1 <= enm2; ++ll)

   511 		{

   512 			l = enm2 + ld - ll - 1;

   513 

   514 			if (l == ld)

   515 				goto L140;

   516 

   517 			lm1 = l - 1;

   518 			l1 = l + 1;

   519 			t = a[l + 1, l + 1];

   520 

   521 			if (Math.Abs(b[l, l]) > epsb)

   522 				t -= sh * b[l, l];

   523 

   524 			if (Math.Abs(a[l, lm1]) <= (Math.Abs(t / a[l1, l])) * epsa)

   525 				goto L100;

   526 		}

   527 

   528 	L140:

   529 		a1 = a11 - sh;

   530 		a2 = a21;

   531 		if (l != ld)

   532 			a[l, lm1] = -a[l, lm1];

   533 		goto L160;

   534 

   535 	// Determine double shift zeroth column of a

   536 	L150:

   537 		a12 = a[l, l1] / b22;

   538 		a22 = a[l1, l1] / b22;

   539 		b12 = b[l, l1] / b22;

   540 		a1 = ((a33 - a11) * (a44 - a11) - a34 * a43 + a43 * b34 * a11) / a21 + a12 - a11 * b12;

   541 		a2 = a22 - a11 - a21 * b12 - (a33 - a11) - (a44 - a11) + a43 * b34;

   542 		a3 = a[l1 + 1, l1] / b22;

   543 		goto L160;

   544 

   545 	// Ad hoc shift

   546 	L155:

   547 		a1 = 0.0;

   548 		a2 = 1.0;

   549 		a3 = 1.1605;

   550 

   551 	L160:

   552 		++its;

   553 		--itn;

   554 

   555 		if (!matz) lor1 = ld;

   556 

   557 		// Main loop

   558 		for (k = l; k <= na; ++k)

   559 		{

   560 			notlas = k != na && ish == 2;

   561 			k1 = k + 1;

   562 			k2 = k + 2;

   563 

   564 			km1 = Math.Max(k, l + 1) - 1; // Computing MAX

   565 			ll = Math.Min(en, k1 + ish);  // Computing MIN

   566 

   567 			if (notlas) goto L190;

   568 

   569 			// Zero a(k+1,k-1)

   570 			if (k == l) goto L170;

   571 			a1 = a[k, km1];

   572 			a2 = a[k1, km1];

   573 

   574 		L170:

   575 			s = Math.Abs(a1) + Math.Abs(a2);

   576 			if (s == 0.0) goto L70;

   577 			u1 = a1 / s;

   578 			u2 = a2 / s;

   579 			r = Special.Sign(Math.Sqrt(u1 * u1 + u2 * u2), u1);

   580 			v1 = -(u1 + r) / r;

   581 			v2 = -u2 / r;

   582 			u2 = v2 / v1;

   583 

   584 			for (j = km1; j < enorn; ++j)

   585 			{

   586 				t = a[k, j] + u2 * a[k1, j];

   587 				a[k, j] += t * v1;

   588 				a[k1, j] += t * v2;

   589 

   590 				t = b[k, j] + u2 * b[k1, j];

   591 				b[k, j] += t * v1;

   592 				b[k1, j] += t * v2;

   593 			}

   594 

   595 			if (k != l)

   596 				a[k1, km1] = 0.0;

   597 			goto L240;

   598 

   599 			// Zero a(k+1,k-1) and a(k+2,k-1)

   600 		L190:

   601 			if (k == l) goto L200;

   602 			a1 = a[k, km1];

   603 			a2 = a[k1, km1];

   604 			a3 = a[k2, km1];

   605 

   606 		L200:

   607 			s = Math.Abs(a1) + Math.Abs(a2) + Math.Abs(a3);

   608 			if (s == 0.0) goto L260;

   609 			u1 = a1 / s;

   610 			u2 = a2 / s;

   611 			u3 = a3 / s;

   612 			r = Special.Sign(Math.Sqrt(u1 * u1 + u2 * u2 + u3 * u3), u1);

   613 			v1 = -(u1 + r) / r;

   614 			v2 = -u2 / r;

   615 			v3 = -u3 / r;

   616 			u2 = v2 / v1;

   617 			u3 = v3 / v1;

   618 

   619 			for (j = km1; j < enorn; ++j)

   620 			{

   621 				t = a[k, j] + u2 * a[k1, j] + u3 * a[k2, j];

   622 				a[k, j] += t * v1;

   623 				a[k1, j] += t * v2;

   624 				a[k2, j] += t * v3;

   625 

   626 				t = b[k, j] + u2 * b[k1, j] + u3 * b[k2, j];

   627 				b[k, j] += t * v1;

   628 				b[k1, j] += t * v2;

   629 				b[k2, j] += t * v3;

   630 			}

   631 

   632 			if (k == l) goto L220;

   633 			a[k1, km1] = 0.0;

   634 			a[k2, km1] = 0.0;

   635 

   636 		// Zero b(k+2,k+1) and b(k+2,k)

   637 		L220:

   638 			s = (Math.Abs(b[k2, k2])) + (Math.Abs(b[k2, k1])) + (Math.Abs(b[k2, k]));

   639 			if (s == 0.0) goto L240;

   640 			u1 = b[k2, k2] / s;

   641 			u2 = b[k2, k1] / s;

   642 			u3 = b[k2, k] / s;

   643 			r = Special.Sign(Math.Sqrt(u1 * u1 + u2 * u2 + u3 * u3), u1);

   644 			v1 = -(u1 + r) / r;

   645 			v2 = -u2 / r;

   646 			v3 = -u3 / r;

   647 			u2 = v2 / v1;

   648 			u3 = v3 / v1;

   649 

   650 			for (i = lor1; i < ll + 1; ++i)

   651 			{

   652 				t = a[i, k2] + u2 * a[i, k1] + u3 * a[i, k];

   653 				a[i, k2] += t * v1;

   654 				a[i, k1] += t * v2;

   655 				a[i, k] += t * v3;

   656 

   657 				t = b[i, k2] + u2 * b[i, k1] + u3 * b[i, k];

   658 				b[i, k2] += t * v1;

   659 				b[i, k1] += t * v2;

   660 				b[i, k] += t * v3;

   661 			}

   662 

   663 			b[k2, k] = 0.0;

   664 			b[k2, k1] = 0.0;

   665 

   666 			if (matz)

   667 			{

   668 				for (i = 0; i < n; ++i)

   669 				{

   670 					t = z[i, k2] + u2 * z[i, k1] + u3 * z[i, k];

   671 					z[i, k2] += t * v1;

   672 					z[i, k1] += t * v2;

   673 					z[i, k] += t * v3;

   674 				}

   675 			}

   676 

   677 		// Zero b(k+1,k)

   678 		L240:

   679 			s = (Math.Abs(b[k1, k1])) + (Math.Abs(b[k1, k]));

   680 			if (s == 0.0) goto L260;

   681 			u1 = b[k1, k1] / s;

   682 			u2 = b[k1, k] / s;

   683 			r = Special.Sign(Math.Sqrt(u1 * u1 + u2 * u2), u1);

   684 			v1 = -(u1 + r) / r;

   685 			v2 = -u2 / r;

   686 			u2 = v2 / v1;

   687 

   688 			for (i = lor1; i < ll + 1; ++i)

   689 			{

   690 				t = a[i, k1] + u2 * a[i, k];

   691 				a[i, k1] += t * v1;

   692 				a[i, k] += t * v2;

   693 

   694 				t = b[i, k1] + u2 * b[i, k];

   695 				b[i, k1] += t * v1;

   696 				b[i, k] += t * v2;

   697 			}

   698 

   699 			b[k1, k] = 0.0;

   700 

   701 			if (matz)

   702 			{

   703 				for (i = 0; i < n; ++i)

   704 				{

   705 					t = z[i, k1] + u2 * z[i, k];

   706 					z[i, k1] += t * v1;

   707 					z[i, k] += t * v2;

   708 				}

   709 			}

   710 

   711 		L260:

   712 			;

   713 		}

   714 

   715 		goto L70; // End QZ step

   716 

   717 	// Set error -- all eigenvalues have not converged after 30*n iterations

   718 	L1000:

   719 		ierr = en + 1;

   720 

   721 	// Save epsb for use by qzval and qzvec

   722 	L1001:

   723 		if (n > 1)

   724 			b[n - 1, 0] = epsb;

   725 		return 0;

   726 	}

   727 

   728 	/// <summary>

   729 	///   Adaptation of the original Fortran QZVAL routine from EISPACK.

   730 	/// </summary>

   731 	/// <remarks>

   732 	///   This subroutine is the third step of the qz algorithm

   733 	///   for solving generalized matrix eigenvalue problems,

   734 	///   siam j. numer. anal. 10, 241-256(1973) by moler and stewart.

   735 	///  

   736 	///   This subroutine accepts a pair of real matrices, one of them

   737 	///   in quasi-triangular form and the other in upper triangular form.

   738 	///   it reduces the quasi-triangular matrix further, so that any

   739 	///   remaining 2-by-2 blocks correspond to pairs of complex

   740 	///   eigenvalues, and returns quantities whose ratios give the

   741 	///   generalized eigenvalues.	it is usually preceded by  qzhes

   742 	///   and  qzit  and may be followed by  qzvec.

   743 	///  

   744 	///   For the full documentation, please check the original function.

   745 	/// </remarks>

   746 	private static int qzval(int n, double[,] a, double[,] b, double[] alfr, double[] alfi, double[] beta, bool matz, double[,] z)

   747 	{

   748 		int i, j;

   749 		int na, en, nn;

   750 		double c, d, e = 0;

   751 		double r, s, t;

   752 		double a1, a2, u1, u2, v1, v2;

   753 		double a11, a12, a21, a22;

   754 		double b11, b12, b22;

   755 		double di, ei;

   756 		double an = 0, bn;

   757 		double cq, dr;

   758 		double cz, ti, tr;

   759 		double a1i, a2i, a11i, a12i, a22i, a11r, a12r, a22r;

   760 		double sqi, ssi, sqr, szi, ssr, szr;

   761 

   762 		double epsb = b[n - 1, 0];

   763 		int isw = 1;

   764 

   765 

   766 		// Find eigenvalues of quasi-triangular matrices.

   767 		for (nn = 0; nn < n; ++nn)

   768 		{

   769 			en = n - nn - 1;

   770 			na = en - 1;

   771 

   772 			if (isw == 2) goto L505;

   773 			if (en == 0) goto L410;

   774 			if (a[en, na] != 0.0) goto L420;

   775 

   776 		// 1-by-1 block, one real root

   777 		L410:

   778 			alfr[en] = a[en, en];

   779 			if (b[en, en] < 0.0)

   780 			{

   781 				alfr[en] = -alfr[en];

   782 			}

   783 			beta[en] = (Math.Abs(b[en, en]));

   784 			alfi[en] = 0.0;

   785 			goto L510;

   786 

   787 		// 2-by-2 block

   788 		L420:

   789 			if (Math.Abs(b[na, na]) <= epsb) goto L455;

   790 			if (Math.Abs(b[en, en]) > epsb) goto L430;

   791 			a1 = a[en, en];

   792 			a2 = a[en, na];

   793 			bn = 0.0;

   794 			goto L435;

   795 

   796 		L430:

   797 			an = Math.Abs(a[na, na]) + Math.Abs(a[na, en]) + Math.Abs(a[en, na]) + Math.Abs(a[en, en]);

   798 			bn = Math.Abs(b[na, na]) + Math.Abs(b[na, en]) + Math.Abs(b[en, en]);

   799 			a11 = a[na, na] / an;

   800 			a12 = a[na, en] / an;

   801 			a21 = a[en, na] / an;

   802 			a22 = a[en, en] / an;

   803 			b11 = b[na, na] / bn;

   804 			b12 = b[na, en] / bn;

   805 			b22 = b[en, en] / bn;

   806 			e = a11 / b11;

   807 			ei = a22 / b22;

   808 			s = a21 / (b11 * b22);

   809 			t = (a22 - e * b22) / b22;

   810 

   811 			if (Math.Abs(e) <= Math.Abs(ei))

   812 				goto L431;

   813 

   814 			e = ei;

   815 			t = (a11 - e * b11) / b11;

   816 

   817 		L431:

   818 			c = (t - s * b12) * .5;

   819 			d = c * c + s * (a12 - e * b12);

   820 			if (d < 0.0) goto L480;

   821 

   822 			// Two real roots. Zero both a(en,na) and b(en,na)

   823 			e += c + Special.Sign(Math.Sqrt(d), c);

   824 			a11 -= e * b11;

   825 			a12 -= e * b12;

   826 			a22 -= e * b22;

   827 

   828 			if (Math.Abs(a11) + Math.Abs(a12) < Math.Abs(a21) + Math.Abs(a22))

   829 				goto L432;

   830 

   831 			a1 = a12;

   832 			a2 = a11;

   833 			goto L435;

   834 

   835 		L432:

   836 			a1 = a22;

   837 			a2 = a21;

   838 

   839 		// Choose and apply real z

   840 		L435:

   841 			s = Math.Abs(a1) + Math.Abs(a2);

   842 			u1 = a1 / s;

   843 			u2 = a2 / s;

   844 			r = Special.Sign(Math.Sqrt(u1 * u1 + u2 * u2), u1);

   845 			v1 = -(u1 + r) / r;

   846 			v2 = -u2 / r;

   847 			u2 = v2 / v1;

   848 

   849 			for (i = 0; i <= en; ++i)

   850 			{

   851 				t = a[i, en] + u2 * a[i, na];

   852 				a[i, en] += t * v1;

   853 				a[i, na] += t * v2;

   854 

   855 				t = b[i, en] + u2 * b[i, na];

   856 				b[i, en] += t * v1;

   857 				b[i, na] += t * v2;

   858 			}

   859 

   860 			if (matz)

   861 			{

   862 				for (i = 0; i < n; ++i)

   863 				{

   864 					t = z[i, en] + u2 * z[i, na];

   865 					z[i, en] += t * v1;

   866 					z[i, na] += t * v2;

   867 				}

   868 			}

   869 

   870 			if (bn == 0.0) goto L475;

   871 			if (an < System.Math.Abs(e) * bn) goto L455;

   872 			a1 = b[na, na];

   873 			a2 = b[en, na];

   874 			goto L460;

   875 

   876 		L455:

   877 			a1 = a[na, na];

   878 			a2 = a[en, na];

   879 

   880 		// Choose and apply real q

   881 		L460:

   882 			s = System.Math.Abs(a1) + System.Math.Abs(a2);

   883 			if (s == 0.0) goto L475;

   884 			u1 = a1 / s;

   885 			u2 = a2 / s;

   886 			r = Special.Sign(Math.Sqrt(u1 * u1 + u2 * u2), u1);

   887 			v1 = -(u1 + r) / r;

   888 			v2 = -u2 / r;

   889 			u2 = v2 / v1;

   890 

   891 			for (j = na; j < n; ++j)

   892 			{

   893 				t = a[na, j] + u2 * a[en, j];

   894 				a[na, j] += t * v1;

   895 				a[en, j] += t * v2;

   896 

   897 				t = b[na, j] + u2 * b[en, j];

   898 				b[na, j] += t * v1;

   899 				b[en, j] += t * v2;

   900 			}

   901 

   902 		L475:

   903 			a[en, na] = 0.0;

   904 			b[en, na] = 0.0;

   905 			alfr[na] = a[na, na];

   906 			alfr[en] = a[en, en];

   907 

   908 			if (b[na, na] < 0.0)

   909 				alfr[na] = -alfr[na];

   910 

   911 			if (b[en, en] < 0.0)

   912 				alfr[en] = -alfr[en];

   913 

   914 			beta[na] = (System.Math.Abs(b[na, na]));

   915 			beta[en] = (System.Math.Abs(b[en, en]));

   916 			alfi[en] = 0.0;

   917 			alfi[na] = 0.0;

   918 			goto L505;

   919 

   920 			// Two complex roots

   921 		L480:

   922 			e += c;

   923 			ei = System.Math.Sqrt(-d);

   924 			a11r = a11 - e * b11;

   925 			a11i = ei * b11;

   926 			a12r = a12 - e * b12;

   927 			a12i = ei * b12;

   928 			a22r = a22 - e * b22;

   929 			a22i = ei * b22;

   930 

   931 			if (System.Math.Abs(a11r) + System.Math.Abs(a11i) +

   932 				System.Math.Abs(a12r) + System.Math.Abs(a12i) <

   933 				System.Math.Abs(a21) + System.Math.Abs(a22r)

   934 				+ System.Math.Abs(a22i))

   935 				goto L482;

   936 

   937 			a1 = a12r;

   938 			a1i = a12i;

   939 			a2 = -a11r;

   940 			a2i = -a11i;

   941 			goto L485;

   942 

   943 		L482:

   944 			a1 = a22r;

   945 			a1i = a22i;

   946 			a2 = -a21;

   947 			a2i = 0.0;

   948 

   949 		// Choose complex z

   950 		L485:

   951 			cz = System.Math.Sqrt(a1 * a1 + a1i * a1i);

   952 			if (cz == 0.0) goto L487;

   953 			szr = (a1 * a2 + a1i * a2i) / cz;

   954 			szi = (a1 * a2i - a1i * a2) / cz;

   955 			r = System.Math.Sqrt(cz * cz + szr * szr + szi * szi);

   956 			cz /= r;

   957 			szr /= r;

   958 			szi /= r;

   959 			goto L490;

   960 

   961 		L487:

   962 			szr = 1.0;

   963 			szi = 0.0;

   964 

   965 		L490:

   966 			if (an < (System.Math.Abs(e) + ei) * bn) goto L492;

   967 			a1 = cz * b11 + szr * b12;

   968 			a1i = szi * b12;

   969 			a2 = szr * b22;

   970 			a2i = szi * b22;

   971 			goto L495;

   972 

   973 		L492:

   974 			a1 = cz * a11 + szr * a12;

   975 			a1i = szi * a12;

   976 			a2 = cz * a21 + szr * a22;

   977 			a2i = szi * a22;

   978 

   979 		// Choose complex q

   980 		L495:

   981 			cq = System.Math.Sqrt(a1 * a1 + a1i * a1i);

   982 			if (cq == 0.0) goto L497;

   983 			sqr = (a1 * a2 + a1i * a2i) / cq;

   984 			sqi = (a1 * a2i - a1i * a2) / cq;

   985 			r = System.Math.Sqrt(cq * cq + sqr * sqr + sqi * sqi);

   986 			cq /= r;

   987 			sqr /= r;

   988 			sqi /= r;

   989 			goto L500;

   990 

   991 		L497:

   992 			sqr = 1.0;

   993 			sqi = 0.0;

   994 

   995 		// Compute diagonal elements that would result if transformations were applied

   996 		L500:

   997 			ssr = sqr * szr + sqi * szi;

   998 			ssi = sqr * szi - sqi * szr;

   999 			i = 0;

  1000 			tr = cq * cz * a11 + cq * szr * a12 + sqr * cz * a21 + ssr * a22;

  1001 			ti = cq * szi * a12 - sqi * cz * a21 + ssi * a22;

  1002 			dr = cq * cz * b11 + cq * szr * b12 + ssr * b22;

  1003 			di = cq * szi * b12 + ssi * b22;

  1004 			goto L503;

  1005 

  1006 		L502:

  1007 			i = 1;

  1008 			tr = ssr * a11 - sqr * cz * a12 - cq * szr * a21 + cq * cz * a22;

  1009 			ti = -ssi * a11 - sqi * cz * a12 + cq * szi * a21;

  1010 			dr = ssr * b11 - sqr * cz * b12 + cq * cz * b22;

  1011 			di = -ssi * b11 - sqi * cz * b12;

  1012 

  1013 		L503:

  1014 			t = ti * dr - tr * di;

  1015 			j = na;

  1016 

  1017 			if (t < 0.0)

  1018 				j = en;

  1019 

  1020 			r = Math.Sqrt(dr * dr + di * di);

  1021 			beta[j] = bn * r;

  1022 			alfr[j] = an * (tr * dr + ti * di) / r;

  1023 			alfi[j] = an * t / r;

  1024 			if (i == 0) goto L502;

  1025 

  1026 		L505:

  1027 			isw = 3 - isw;

  1028 

  1029 		L510:

  1030 			;

  1031 		}

  1032 

  1033 		b[n - 1, 0] = epsb;

  1034 

  1035 		return 0;

  1036 	}

  1037 

  1038 	/// <summary>

  1039 	///   Adaptation of the original Fortran QZVEC routine from EISPACK.

  1040 	/// </summary>

  1041 	/// <remarks>

  1042 	///   This subroutine is the optional fourth step of the qz algorithm

  1043 	///   for solving generalized matrix eigenvalue problems,

  1044 	///   siam j. numer. anal. 10, 241-256(1973) by moler and stewart.

  1045 	///  

  1046 	///   This subroutine accepts a pair of real matrices, one of them in

  1047 	///   quasi-triangular form (in which each 2-by-2 block corresponds to

  1048 	///   a pair of complex eigenvalues) and the other in upper triangular

  1049 	///   form.  It computes the eigenvectors of the triangular problem and

  1050 	///   transforms the results back to the original coordinate system.

  1051 	///   it is usually preceded by  qzhes,  qzit, and	qzval.

  1052 	///  

  1053 	///   For the full documentation, please check the original function.

  1054 	/// </remarks>

  1055 	private static int qzvec(int n, double[,] a, double[,] b, double[] alfr, double[] alfi, double[] beta, double[,] z)

  1056 	{

  1057 		int i, j, k, m;

  1058 		int na, ii, en, jj, nn, enm2;

  1059 		double d, q;

  1060 		double r = 0, s = 0, t, w, x = 0, y, t1, t2, w1, x1 = 0, z1 = 0, di;

  1061 		double ra, dr, sa;

  1062 		double ti, rr, tr, zz = 0;

  1063 		double alfm, almi, betm, almr;

  1064 

  1065 		double epsb = b[n - 1, 0];

  1066 		int isw = 1;

  1067 

  1068 

  1069 		// for en=n step -1 until 1 do --

  1070 		for (nn = 0; nn < n; ++nn)

  1071 		{

  1072 			en = n - nn - 1;

  1073 			na = en - 1;

  1074 			if (isw == 2) goto L795;

  1075 			if (alfi[en] != 0.0) goto L710;

  1076 

  1077 			// Real vector

  1078 			m = en;

  1079 			b[en, en] = 1.0;

  1080 			if (na == -1) goto L800;

  1081 			alfm = alfr[m];

  1082 			betm = beta[m];

  1083 

  1084 			// for i=en-1 step -1 until 1 do --

  1085 			for (ii = 0; ii <= na; ++ii)

  1086 			{

  1087 				i = en - ii - 1;

  1088 				w = betm * a[i, i] - alfm * b[i, i];

  1089 				r = 0.0;

  1090 

  1091 				for (j = m; j <= en; ++j)

  1092 					r += (betm * a[i, j] - alfm * b[i, j]) * b[j, en];

  1093 

  1094 				if (i == 0 || isw == 2)

  1095 					goto L630;

  1096 

  1097 				if (betm * a[i, i - 1] == 0.0)

  1098 					goto L630;

  1099 

  1100 				zz = w;

  1101 				s = r;

  1102 				goto L690;

  1103 

  1104 			L630:

  1105 				m = i;

  1106 				if (isw == 2) goto L640;

  1107 

  1108 				// Real 1-by-1 block

  1109 				t = w;

  1110 				if (w == 0.0)

  1111 					t = epsb;

  1112 				b[i, en] = -r / t;

  1113 				goto L700;

  1114 

  1115 			// Real 2-by-2 block

  1116 			L640:

  1117 				x = betm * a[i, i + 1] - alfm * b[i, i + 1];

  1118 				y = betm * a[i + 1, i];

  1119 				q = w * zz - x * y;

  1120 				t = (x * s - zz * r) / q;

  1121 				b[i, en] = t;

  1122 				if (Math.Abs(x) <= Math.Abs(zz)) goto L650;

  1123 				b[i + 1, en] = (-r - w * t) / x;

  1124 				goto L690;

  1125 

  1126 			L650:

  1127 				b[i + 1, en] = (-s - y * t) / zz;

  1128 

  1129 			L690:

  1130 				isw = 3 - isw;

  1131 

  1132 			L700:

  1133 				;

  1134 			}

  1135 			// End real vector