% FIGURE 3.2.1
      function [A,X] = DEGSYM(A,N)
%
%  FUNCTION DEGSYM SOLVES THE EIGENVALUE PROBLEM
%
%                A*X = LAMBDA*X
%
%    WHERE A IS A SYMMETRIC MATRIX.
%
%
%  ARGUMENTS
%
%             ON INPUT                          ON OUTPUT
%             --------                          ---------
%
%    A      - THE N BY N SYMMETRIC MATRIX.      A DIAGONAL MATRIX,
%                                               WITH THE EIGENVALUES
%                                               OF A ON THE DIAGONAL.
%
%    N      - THE SIZE OF MATRIX A.
%
%    X      -                                   AN N BY N MATRIX WHICH
%                                               CONTAINS THE EIGEN-
%                                               VECTORS OF A IN ITS
%                                               COLUMNS, IN THE SAME
%                                               ORDER AS THE EIGENVALUES
%                                               APPEAR ON THE DIAGONAL.
%
%-----------------------------------------------------------------------
%                              EPS = MACHINE FLOATING POINT RELATIVE
%                                    PRECISION
% *****************************
      EPS = eps;
% *****************************
%                              ANORM = SUM OF ALL SQUARES
%                              X INITIALIZED TO IDENTITY
      ANORM = 0.0;
      for I=1:N
         for J=1:N
            ANORM = ANORM + A(I,J)^2;
            X(I,J) = 0.0;
         end
         X(I,I) = 1.0;
      end
      ERRLIM = 1000*EPS*ANORM;
%                              EK = SUM OF OFF-DIAGONAL SQUARES
      EK = 0.0;
      for I=1:N
      for J=1:N
         if (I ~= J)
            EK = EK + A(I,J)^2;
         end
      end
      end
      if (EK <= ERRLIM)
         return
      end
      THRESH = 0.5*EK/N/(N-1);
      while (1 > 0)
         for I=1:N-1
            for J=I+1:N
%                              IF A(J,I)^2 LESS THAN HALF THE
%                              AVERAGE FOR OFF-DIAGONALS, SKIP IT.
               if (A(J,I)^2 <= THRESH)
                  continue
               end
%                              KNOCKING OUT A(J,I) WILL DECREASE OFF-
%                              DIAGONAL SUM OF SQUARES BY 2*A(J,I)^2.
               EK = EK - 2*A(J,I)^2;
%                              CALCULATE NEW THRESHOLD.
               THRESH = 0.5*EK/N/(N-1);
%                              CALCULATE C,S
               BETA = (A(I,I)-A(J,J))/(2.*A(J,I));
               FRACT = 0.5*BETA/sqrt(1.0+BETA^2);
               S = sqrt(max(0.5-FRACT,0.0));
               C = sqrt(max(0.5+FRACT,0.0));
%                              PREMULTIPLY A BY Qij^T
               for K=1:N
                  PIK =  C*A(I,K)+S*A(J,K);
                  PJK = -S*A(I,K)+C*A(J,K);
                  A(I,K) = PIK;
                  A(J,K) = PJK;
               end
%                              POSTMULTIPLY A AND X BY Qij
               for K=1:N
                  BKI =  C*A(K,I)+S*A(K,J);
                  BKJ = -S*A(K,I)+C*A(K,J);
                  A(K,I) = BKI;
                  A(K,J) = BKJ;
                  XKI =  C*X(K,I)+S*X(K,J);
                  XKJ = -S*X(K,I)+C*X(K,J);
                  X(K,I) = XKI;
                  X(K,J) = XKJ;
               end
%                              CHECK FOR CONVERGENCE
               if (EK <= ERRLIM)
                  return
               end
            end
         end
%                              RETURN TO BEGINNING OF CYCLE
      end