/*****************************************************************************
*                                                                            *
*   EIG.SPL      Copyright (C) 2001-2014 Development Corporation             *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    Returns eigenvalues and eigenvectors of an array            *
*                                                                            *
*   Revisions:   24 Jul 2001  RRR  Creation                                  *
*                 8 Jul 2011  RRR  use built-in                              *
*                31 Aug 2012  RRR  generalized eigenvalues                   *
*                22 Jun 2014  RRR  explicit "hess" flag                      *
*                                                                            *
*****************************************************************************/

#if @HELP_EIG

    Purpose: Computes the Eigenvalues and Eigenvectors of a square matrix.

    Syntax:  EIG(a, "opt")

             (vec, val) = EIG(a)

               a   - A square matrix.

             "opt" - Optional. A string, the calculation method.

                       "nobalance": prevents the preliminary matrix
                                    balancing step

                       "hess":      use the Hessenberg form

    Alternate Syntax:  

             EIG(a, b)

             (vec, val) = EIG(a, b)

               a   - A square matrix.

               b   - A square matrix.

    Returns: A single column real or complex series that contains the
             eigenvalues.

             (vec, val) = EIG(a)

             returns both the eigenvector array and the eigenvalues as the
             diagonal of an array.

             (vec, val) = EIG(a, b)

             returns the generalized eigenvectors and eigenvalues for
             arrays a and b.

    Example:
              a = {{1, 3, 4},
                   {5, 6, 7},
                   {8, 9, 12}}

              val = eig(a)

              val == {19.964160,
                      -1.473921,
                       0.509760}

    Example:
             a = {{1, 3, 4},
                  {5, 6, 7},
                  {8, 9, 12}}

             (x, lambda) = eig(a)


                  x == {{ -0.253874, -0.896277,  0.046508},
                        { -0.504564,  0.270278, -0.801862},
                        { -0.825205,  0.351621,  0.595697}}

             lambda == {{ 19.964160,  0.000000,  0.000000},
                        {  0.000000, -1.473921,  0.000000},
                        {  0.000000,  0.000000,  0.509760}}


              a *^ x == x *^ lambda ==

                       {{ -5.068385,  1.321041,  0.023708},
                        {-10.073188, -0.398369, -0.408757},
                        {-16.474527, -0.518261,  0.303662}}

    Remarks:
              For (x, lambda) = eig(a), the nth entry of the diagonal
              of the eigenvalue array lambda corresponds to the eigenvector
              in column n of x.

              For (v, d) = eig(a, b), 

                a *^ v == b *^ v *^ d

              and eig(a, b) == eig(inv(b) *^ a).

    See Also:
             *^
             Balance
             Eigval
             Nbeigval
             Nbeigvec
#endif


/* eigenvalues and eigenvectors */
eig(a, b, opt)
{
        local V, D, H;

        if (not(isarray(a))) error("eig - input array required");

        if (argc < 3)
        {
                if (argc < 2)
                {
                        b = {};
                }

                opt = "";
        }

        if (isstring(b))
        {
                /* eig(A, opt) */
                opt = b;
                b   = {};
        }

        if (not(isarray(b))) error("eig - input array required");

        if (not(isstring(opt))) opt = "";
                        
        opt = tolower(opt);

        if (outargc > 1)
        {
                if (opt == "nobalance")
                {
                        /* eigenvalue without balancing */
                        (V, D) = _eig(a, "nobalance");
                }
                else if (opt == "hess")
                {
                        /* eigenvalue using Hessenburg form */
                        (Q, H) = hess(a);
                        (V, D) = _eig(H);
                        V = Q *^ V;
                }
                else
                {
                        if (isempty(b))
                        {
                                /* standard */
                                (V, D) = _eig(a);
                        }
                        else
                        {
                                /* generalized */
                                (V, D) = _eig(a, b);
                        }
                }
                
                return(V, D);
        }
        else
        {
                if (isempty(b))
                {
                        /* standard */
                        return(_eig(a));
                }
                else
                {
                        /* generalized */
                        return(_eig(a, b));
                }
        }
}