public void value2x2 ( double a11 , double a12 , double a21 , double a22 ) { // apply a rotators such that th a11 and a22 elements are the same double c , s ; if ( a12 + a21 == 0 ) { // is this pointless since c = s = 1.0 / Math . sqrt ( 2 ) ; } else { double aa = ( a11 - a22 ) ; double bb = ( a12 + a21 ) ; double t_hat = aa / bb ; double t = t_hat / ( 1.0 + Math . sqrt ( 1.0 + t_hat * t_hat ) ) ; c = 1.0 / Math . sqrt ( 1.0 + t * t ) ; s = c * t ; } double c2 = c * c ; double s2 = s * s ; double cs = c * s ; double b11 = c2 * a11 + s2 * a22 - cs * ( a12 + a21 ) ; double b12 = c2 * a12 - s2 * a21 + cs * ( a11 - a22 ) ; double b21 = c2 * a21 - s2 * a12 + cs * ( a11 - a22 ) ; // double b22 = c2*a22 + s2*a11 + cs*(a12+a21); // apply second rotator to make A upper triangular if real eigenvalues if ( b21 * b12 >= 0 ) { if ( b12 == 0 ) { c = 0 ; s = 1 ; } else { s = Math . sqrt ( b21 / ( b12 + b21 ) ) ; c = Math . sqrt ( b12 / ( b12 + b21 ) ) ; } // c2 = b12;//c*c; // s2 = b21;//s*s; cs = c * s ; a11 = b11 - cs * ( b12 + b21 ) ; // a12 = c2*b12 - s2*b21; // a21 = c2*b21 - s2*b12; a22 = b11 + cs * ( b12 + b21 ) ; value0 . real = a11 ; value1 . real = a22 ; value0 . imaginary = value1 . imaginary = 0 ; } else { value0 . real = value1 . real = b11 ; value0 . imaginary = Math . sqrt ( - b21 * b12 ) ; value1 . imaginary = - value0 . imaginary ; } }
if |a11 - a22| >> |a12 + a21| there might be a better way . see pg371