// function computes all eigenvalues and eigenvectors of a real symmetric matrix A[0..n-1][0..n-1]
// stored in the top left of a 10x10 array A[10][10]
// A[][] is changed on output.
// eigval[0..n-1] returns the eigenvalues of A[][].
// eigvec[0..n-1][0..n-1] returns the normalized eigenvectors of A[][]
// the eigenvectors are not sorted by value
void eigencompute(float A[][10], float eigval[], float eigvec[][10], int8 n)
{
    // maximum number of iterations to achieve convergence: in practice 6 is typical
#define NITERATIONS 15

    // various trig functions of the jacobi rotation angle phi
    float cot2phi, tanhalfphi, tanphi, sinphi, cosphi;
    // scratch variable to prevent over-writing during rotations
    float ftmp;
    // residue from remaining non-zero above diagonal terms
    float residue;
    // matrix row and column indices
    int8 ir, ic;
    // general loop counter
    int8 j;
    // timeout ctr for number of passes of the algorithm
    int8 ctr;

    // initialize eigenvectors matrix and eigenvalues array
    for (ir = 0; ir < n; ir++)
    {
        // loop over all columns
        for (ic = 0; ic < n; ic++)
        {
            // set on diagonal and off-diagonal elements to zero
            eigvec[ir][ic] = 0.0F;
        }

        // correct the diagonal elements to 1.0
        eigvec[ir][ir] = 1.0F;

        // initialize the array of eigenvalues to the diagonal elements of m
        eigval[ir] = A[ir][ir];
    }

    // initialize the counter and loop until converged or NITERATIONS reached
    ctr = 0;
    do
    {
        // compute the absolute value of the above diagonal elements as exit criterion
        residue = 0.0F;
        // loop over rows excluding last row
        for (ir = 0; ir < n - 1; ir++)
        {
            // loop over above diagonal columns
            for (ic = ir + 1; ic < n; ic++)
            {
                // accumulate the residual off diagonal terms which are being driven to zero
                residue += fabs(A[ir][ic]);
            }
        }

        // check if we still have work to do
        if (residue > 0.0F)
        {
            // loop over all rows with the exception of the last row (since only rotating above diagonal elements)
            for (ir = 0; ir < n - 1; ir++)
            {
                // loop over columns ic (where ic is always greater than ir since above diagonal)
                for (ic = ir + 1; ic < n; ic++)
                {
                    // only continue with this element if the element is non-zero
                    if (fabs(A[ir][ic]) > 0.0F)
                    {
                        // calculate cot(2*phi) where phi is the Jacobi rotation angle
                        cot2phi = 0.5F * (eigval[ic] - eigval[ir]) / (A[ir][ic]);

                        // calculate tan(phi) correcting sign to ensure the smaller solution is used
                        tanphi = 1.0F / (fabs(cot2phi) + sqrtf(1.0F + cot2phi * cot2phi));
                        if (cot2phi < 0.0F)
                        {
                            tanphi = -tanphi;
                        }

                        // calculate the sine and cosine of the Jacobi rotation angle phi
                        cosphi = 1.0F / sqrtf(1.0F + tanphi * tanphi);
                        sinphi = tanphi * cosphi;

                        // calculate tan(phi/2)
                        tanhalfphi = sinphi / (1.0F + cosphi);

                        // set tmp = tan(phi) times current matrix element used in update of leading diagonal elements
                        ftmp = tanphi * A[ir][ic];

                        // apply the jacobi rotation to diagonal elements [ir][ir] and [ic][ic] stored in the eigenvalue array
                        // eigval[ir] = eigval[ir] - tan(phi) *  A[ir][ic]
                        eigval[ir] -= ftmp;
                        // eigval[ic] = eigval[ic] + tan(phi) * A[ir][ic]
                        eigval[ic] += ftmp;

                        // by definition, applying the jacobi rotation on element ir, ic results in 0.0
                        A[ir][ic] = 0.0F;

                        // apply the jacobi rotation to all elements of the eigenvector matrix
                        for (j = 0; j < n; j++)
                        {
                            // store eigvec[j][ir]
                            ftmp = eigvec[j][ir];
                            // eigvec[j][ir] = eigvec[j][ir] - sin(phi) * (eigvec[j][ic] + tan(phi/2) * eigvec[j][ir])
                            eigvec[j][ir] = ftmp - sinphi * (eigvec[j][ic] + tanhalfphi * ftmp);
                            // eigvec[j][ic] = eigvec[j][ic] + sin(phi) * (eigvec[j][ir] - tan(phi/2) * eigvec[j][ic])
                            eigvec[j][ic] = eigvec[j][ic] + sinphi * (ftmp - tanhalfphi * eigvec[j][ic]);
                        }

                        // apply the jacobi rotation only to those elements of matrix m that can change
                        for (j = 0; j <= ir - 1; j++)
                        {
                            // store A[j][ir]
                            ftmp = A[j][ir];
                            // A[j][ir] = A[j][ir] - sin(phi) * (A[j][ic] + tan(phi/2) * A[j][ir])
                            A[j][ir] = ftmp - sinphi * (A[j][ic] + tanhalfphi * ftmp);
                            // A[j][ic] = A[j][ic] + sin(phi) * (A[j][ir] - tan(phi/2) * A[j][ic])
                            A[j][ic] = A[j][ic] + sinphi * (ftmp - tanhalfphi * A[j][ic]);
                        }
                        for (j = ir + 1; j <= ic - 1; j++)
                        {
                            // store A[ir][j]
                            ftmp = A[ir][j];
                            // A[ir][j] = A[ir][j] - sin(phi) * (A[j][ic] + tan(phi/2) * A[ir][j])
                            A[ir][j] = ftmp - sinphi * (A[j][ic] + tanhalfphi * ftmp);
                            // A[j][ic] = A[j][ic] + sin(phi) * (A[ir][j] - tan(phi/2) * A[j][ic])
                            A[j][ic] = A[j][ic] + sinphi * (ftmp - tanhalfphi * A[j][ic]);
                        }
                        for (j = ic + 1; j < n; j++)
                        {
                            // store A[ir][j]
                            ftmp = A[ir][j];
                            // A[ir][j] = A[ir][j] - sin(phi) * (A[ic][j] + tan(phi/2) * A[ir][j])
                            A[ir][j] = ftmp - sinphi * (A[ic][j] + tanhalfphi * ftmp);
                            // A[ic][j] = A[ic][j] + sin(phi) * (A[ir][j] - tan(phi/2) * A[ic][j])
                            A[ic][j] = A[ic][j] + sinphi * (ftmp - tanhalfphi * A[ic][j]);
                        }
                    }   // end of test for matrix element already zero
                }   // end of loop over columns
            }   // end of loop over rows
        }  // end of test for non-zero residue
    } while ((residue > 0.0F) && (ctr++ < NITERATIONS)); // end of main loop

    return;
}