% FIGURE 3.4.4
      function A = HESSM(A,N)
      if (N <= 2)
         return
      end
%                              USE Mij TRANSFORMATIONS TO REDUCE A
%                              TO UPPER HESSENBERG FORM
      for I=2:N-1
%                              SEARCH FROM A(I,I-1) ON DOWN FOR
%                              LARGEST POTENTIAL PIVOT, A(L,I-1)
         BIG = abs(A(I,I-1));
         L = I;
         for J=I+1:N
            if (abs(A(J,I-1)) > BIG)
               BIG = abs(A(J,I-1));
               L = J;
            end
         end
%                              IF ALL SUBDIAGONAL ELEMENTS IN COLUMN
%                              I-1 ALREADY ZERO, GO ON TO NEXT COLUMN
         if (BIG == 0.0)
            continue
         end
%                              PREMULTIPLY BY Pil
%                              (SWITCH ROWS I AND L)
         for K=I-1:N
            TEMP = A(L,K);
            A(L,K) = A(I,K);
            A(I,K) = TEMP;
         end
%                              POSTMULTIPLY BY Pil^(-1) = Pil
%                              (SWITCH COLUMNS I AND L)
         for K=1:N
            TEMP = A(K,L);
            A(K,L) = A(K,I);
            A(K,I) = TEMP;
         end
         for J=I+1:N
            R = A(J,I-1)/A(I,I-1);
            if (R == 0.0)
               continue
            end
%                              PREMULTIPLY BY Mij^(-1)
%                              (SUBTRACT R TIMES ROW I FROM ROW J)
            for K=I-1:N
               A(J,K) = A(J,K) - R*A(I,K);
            end
%                              POSTMULTIPLY BY Mij
%                              (ADD R TIMES COLUMN J TO COLUMN I)
            for K=1:N
               A(K,I) = A(K,I) + R*A(K,J);
            end
         end
      end