<!--Define JavaScript functions. -->

var N = 2; //Global variable, dimension of matrix.

function cdivA(ar, ai, br, bi, A, in1, in2, in3){
// Division routine for dividing one complex number into another:
// This routine does (ar + ai)/(br + bi) and returns the results in the specified
// elements of the A matrix.

 var s, ars, ais, brs, bis;

 s = Math.abs(br) + Math.abs(bi);
 ars = ar/s;
 ais = ai/s;
 brs = br/s;
 bis = bi/s;
 s = brs*brs + bis*bis;
 A[in1][in2] = (ars*brs + ais*bis)/s;
 A[in1][in3] = (-(ars*bis) + ais*brs)/s;
 return;
} // End cdivA

function hqr2(A, B, igh, wi, wr){
/* Computes the eigenvalues and eigenvectors of a real upper-Hessenberg Matrix using the QR method. */

var norm, p, q, r, s, t, tst1, tst2, w, x, y, zz;
var k = 0, l, m, en = igh, enm2, na;

norm = Math.abs(A[0][0]) + Math.abs(A[0][1]) + Math.abs(A[1][0])+ Math.abs(A[1][1]);

if (igh == 0){  // Store eigenvalue already isolated and compute matrix norm.
  wi[1] = 0.0;
  wr[1] = A[1][1];
 } //End if (igh == 0)

// Search for next eigenvalues
while (en >= 0){
 na = en - 1;

  /*Look for single small sub-diagonal element for l = en step -1 until 0 */

  l = en;
  if (l != 0){
   s = Math.abs(A[l - 1][l - 1]) + Math.abs(A[l][l]);
   if (s == 0.0)
    s = norm;
   tst1 = s;
   tst2 = tst1 + Math.abs(A[l][l - 1]);
   if (tst2 != tst1)
    l--;
  } //End if l != 0)

 x = A[en][en];

 if (l == en){  //One root found
  wr[en] = A[en][en] = x;
  wi[en] = 0.0;
  en--;
  continue;
 } //End if (l == en)

 y = A[na][na];
 w = A[en][na]*A[na][en];
 
 if (l == na){  //Two roots found
  p = (-x + y)/2;
  q = p*p + w;
  zz = Math.sqrt(Math.abs(q));
  A[en][en] = x;
  A[na][na] = y;
  if (q >= 0.0){//Real Pair
   zz = ((p < 0.0) ? (-zz + p) : (p + zz));
   wr[en] = wr[na] = x + zz;
   if (zz != 0.0)
    wr[en] = -(w/zz) + x;
   wi[en] = wi[na] = 0.0;
   x = A[en][na];
   s = Math.abs(x) + Math.abs(zz);
   p = x/s;
   q = zz/s;
   r = Math.sqrt(p*p + q*q);
   p /= r;
   q /= r;
   for (var j = na; j < N; j++){ //Row modification
    zz = A[na][j];
    A[na][j] = q*zz + p*A[en][j];
    A[en][j] = -(p*zz) + q*A[en][j];
   }//End for j
   for (var j = 0; j <= en; j++){ // Column modification
    zz = A[j][na];
    A[j][na] = q*zz + p*A[j][en];
    A[j][en] = -(p*zz) + q*A[j][en];
   }//End for j
   for (var j = 0; j <= igh; j++){//Accumulate transformations
    zz = B[j][na];
    B[j][na] = q*zz + p*B[j][en];
    B[j][en] = -(p*zz) + q*B[j][en];
   }//End for j
  } //End if (q >= 0.0)
  else {//else q < 0.0, Complex Pair
   wr[en] = wr[na] = x + p;
   wi[na] = zz;
   wi[en] = -zz;
  } //End else
  en--;
  en--;
 }//End if (l == na)
 
}//End while (en >= 0)

if (norm == 0.0)
 return;

//Step from (N - 1) to 0 in steps of -1

for (var i = 0; i < N; i++){
 en = N - 1 - i;
 p = wr[en];
 q = wi[en];
 na = en - 1;

 if (q > 0.0)
  continue;

 if (q == 0.0){//Real vector
  m = en;
  A[en][en] = 1.0;

  for (var j = na; j >= 0; j--){
   w = -p + A[j][j];
   r = 0.0;
   for (var ii = m; ii <= en; ii++)
    r += A[j][ii]*A[ii][en];

   if (wi[j] < 0.0){
    zz = w;
    s = r;
   }//End wi[j] < 0.0

   else {//wi[j] >= 0.0
    m = j;
    if (wi[j] == 0.0){
     t = w;
     if (t == 0.0){
      t = tst1 = norm;
      do {
       t *= 0.01;
       tst2 = norm + t;
      } while (tst2 > tst1);
     } //End if t == 0.0
     A[j][en] = -(r/t);
    }//End if wi[j] == 0.0

    else { //wi[j] > 0.0; Solve real equations
     x = A[j][j + 1];
     y = A[j + 1][j];
     q = (-p + wr[j])*(-p + wr[j]) + wi[j]*wi[j];
     A[j][en] = t = (-(zz*r) + x*s)/q;
     A[j + 1][en] = ((Math.abs(x) > Math.abs(zz)) ? -(r + w*t)/x : -(s + y*t)/zz);
    }//End  else wi[j] > 0.0

    // Overflow control
    t = Math.abs(A[j][en]);
    if (t == 0.0)
     continue; //go up to top of for j loop
    tst1 = t;
    tst2 = tst1 + 1.0/tst1;
    if (tst2 > tst1)
     continue; //go up to top of for j loop
    for (var ii = j; ii <= en; ii++)
     A[ii][en] /= t;

   }//End else wi[j] >= 0.0

  }//End for j
  
 }      //End q == 0.0

 else {//else q < 0.0, complex vector
 //Last vector component chosen imaginary so that eigenvector matrix is triangular
 m = na;

 if (Math.abs(A[en][na]) <= Math.abs(A[na][en]))
  cdivA(0.0, -A[na][en], -p + A[na][na], q, A, na, na, en);
 else {
  A[na][na] = q/A[en][na];
  A[na][en] = -(-p + A[en][en])/A[en][na];
 } //End else (Math.abs(A[en][na] > Math.abs(A[na][en])

 A[en][na] = 0.0;
 A[en][en] = 1.0;

 }//End else q < 0.0
 
}//End for i

//End back substitution. Vectors of isolated roots.

if (igh == 0)
 B[1][1] = A[1][1];

// Multiply by transformation matrix to give vectors of original full matrix.

//Step from (N - 1) to 0 in steps of -1.

for (var i = (N - 1); i >= 0; i--){
 m = ((i < igh) ? i : igh);
 
 for (var ii = 0; ii <= igh; ii++){
  zz = 0.0;
  for (var jj = 0; jj <= m; jj++)
   zz += B[ii][jj]*A[jj][i];
  B[ii][i] = zz;
 }//End for ii
}//End for i loop

return;
} //End of function hqr2

function norVecC(Z, wi){// Normalizes the eigenvectors

// This subroutine is based on the LINPACK routine SNRM2, written 25 October 1982, modified
// on 14 October 1993 by Sven Hammarling of Nag Ltd.
// I have further modified it for use in this Javascript routine, for use with a column
// of an array rather than a column vector.
//
// Z - eigenvector Matrix
// wi - eigenvalue vector

 var scale, ssq, absxi, dummy, norm;

 for (var j = 0; j < N; j++){ //Go through the columns of the vector array 
  scale = 0.0;
  ssq = 1.0;

  for (var i = 0; i < N; i++){
   if (Z[i][j] != 0){
    absxi = Math.abs(Z[i][j]);
    dummy = scale/absxi;
    if (scale < absxi){
     ssq = 1.0 + ssq*dummy*dummy;
     scale = absxi;
    }//End if (scale < absxi)
    else
     ssq += 1.0/dummy/dummy;
   }//End if (Z[i][j] != 0)
  } //End for i

  if (wi[j] != 0){// If complex eigenvalue, take into account imaginary part of eigenvector
   for (var i = 0; i < N; i++){
    if (Z[i][j + 1] != 0){
     absxi = Math.abs(Z[i][j + 1]);
     dummy = scale/absxi;
     if (scale < absxi){
      ssq = 1.0 + ssq*dummy*dummy;
      scale = absxi;
     }//End if (scale < absxi)
     else
      ssq += 1.0/dummy/dummy;
     }//End if (Z[i][j + 1] != 0)
   } //End for i
  }//End if (wi[j] != 0)

  norm = scale*Math.sqrt(ssq); //This is the norm of the (possibly complex) vector

  for (var i = 0; i < N; i++)
   Z[i][j] /= norm;

  if (wi[j] != 0){// If complex eigenvalue, also scale imaginary part of eigenvector
   j++;
   for (var i = 0; i < N; i++)
   Z[i][j] /= norm;
  }//End if (wi[j] != 0)

 }// End for j
 
 return;
} // End norVecC

function Eig2Solve(dataForm){

var A = new Array(N);
var B = new Array(N);
for (var i = 0; i < N; i++){
 A[i] = new Array(N);
 B[i] = new Array(N);
}

A[0][0] = parseFloat(dataForm.a11.value);
A[0][1] = parseFloat(dataForm.a12.value);

A[1][0] = parseFloat(dataForm.a21.value);
A[1][1] = parseFloat(dataForm.a22.value);

var wi = new Array(N);
var wr = new Array(N);

/* Balance the matrix to improve accuracy of eigenvalues. Introduces no rounding errors, since it scales A by powers of the radix.
*/

var scale = new Array(N);    //Contains information about transformations.
var trace = new Array(N);    //Records row and column interchanges
var radix = 2;               //Base of machine floating point representation.
var c, f, g, r, s, b2 = radix*radix;
var igh, l = N - 1, noconv;

if (A[1][0] == 0.0){
 scale[1] = 1;
 scale[0] = 0;
 l--;
}//End if (A[1][0] == 0.0)

else{
 if (A[0][1] == 0.0){
  scale[0] = scale[1] = 0;
  l--;

  //Swap rows and columns

  f = A[0][0];
  A[0][0] = A[1][1];
  A[1][1] = f;

  f = A[1][0];
  A[1][0] = A[0][1];
  A[0][1] = f;
  
 }//End if (A[0][1] == 0.0)

}//End else (A[1][0] != 0.0)

if (l > 0){

 //Balance the matrix in both rows.

 scale[1] = scale[0] = 1.0;

 //Iterative loop for norm reduction
 do {
  noconv = 0;
  for (var i = 0; i <= l; i++) {
   c = r = 0.0;
   for (var j = 0; j <= l; j++){
    if (j == i) continue;
    c += Math.abs(A[j][i]);
    r += Math.abs(A[i][j]);
   } // End for j
   if ((c == 0.0) || (r == 0.0)) continue;   //guard against zero c or r due to underflow
   g = r/radix;
   f = 1.0;
   s = c + r;
   while (c < g) {
    f *= radix;
    c *= b2;
   } // End while (c < g)
   g = r*radix;
   while (c >= g) {
    f /= radix;
    c /= b2;
   } // End while (c >= g)

   //Now balance
   if ((c + r)/f < 0.95*s) {
    g = 1.0/f;
    scale[i] *= f;
    noconv = 1;
    for (var j = 0; j < N; j++)
     A[i][j] *= g;
    for (var j = 0; j <= l; j++)
     A[j][i] *= f;
   } //End if ((c + r)/f < 0.95*s)
  } // End for i
 } while (noconv);  // End of do-while loop.

}//End if (l > 0)

igh = l;

/* Accumulate the stabilized elementary similarity transformations used in the reduction of A to upper-Hessenberg form. Introduces no rounding errors since it only transfers the multipliers used in the reduction process into the eigenvector matrix. */

//Initialize B to the identity Matrix
B[1][1] = B[0][0] = 1.0;
B[1][0] = B[0][1] = 0.0;

hqr2(A, B, igh, wi, wr);

dataForm.lam1Re.value = wr[0];
dataForm.lam1Im.value = wi[0];
dataForm.lam2Re.value = wr[1];
dataForm.lam2Im.value = wi[1];

if (igh != 0){
 for (var i = 0; i <= igh; i++){
  s = scale[i];
  for (var j = 0; j < N; j++)
   B[i][j] *= s;
 }//End for i
}//End if (igh != 0)

for (var i = (igh + 1); i < N; i++){
 k = scale[i];
 if (k != i){
  for (var j = 0; j < N; j++){
   s = B[i][j];
   B[i][j] = B[k][j];
   B[k][j] = s;
  }//End for j
 }//End if k != i
}//End for i

norVecC(B, wi);  //Normalize the eigenvectors

dataForm.v11Re.value = B[0][0];
dataForm.v12Re.value = B[1][0];

if (wi[0] == 0.0){ // Real eigenvalues; therefore, real eigenvectors
 dataForm.v22Im.value = dataForm.v21Im.value = dataForm.v12Im.value = dataForm.v11Im.value = 0;
 dataForm.v21Re.value = B[0][1];
 dataForm.v22Re.value = B[1][1];
} // End if (wi[0] == 0.0)
else { // Complex eigenvalues; therefore, complex eigenvectors
 dataForm.v21Re.value = B[0][0];
 dataForm.v22Re.value = B[1][0];
 dataForm.v21Im.value = dataForm.v11Im.value = B[0][1];
 dataForm.v22Im.value = dataForm.v12Im.value = B[1][1];
} // End else (wi[0] != 0.0)

return;
}  //End of Eig2Solve

// end of JavaScript-->