#
# program that find eigenvalue matrix and 
# eigenvector matrix of a given symmetric matrix
#
# That is: for given symmetric matrix A, the program
# finds matrix Lambda and X such that
#    
#      A = X*Lambda*(transpose(X))
#
# where transpose(X)*X = I
#
eigenmat:=proc(
               A::matrix,  # input matrix
          Lambda::evaln,   # output: eigenvalue matrix
               X::evaln    # output: eigenvector matrix
             )
   local n, s, i, j, k, l, ll, m, u, v, lambda;

   ### get matrix size
   n:=linalg[rowdim](A);

   ### quick return if not square                   
   if n<> linalg[coldim](A) then           
      print(`expecting a square matrix`);
      RETURN();
   fi;

   ### get eigenelements as a list
   s:=[linalg[eigenvects](A)];     

   ### number of distinct eigenvalues        
   m:=nops(s);                             

   ### open matrices for output
   X:=matrix(n,n);
   Lambda:=matrix(n,n,[seq(0,i=1..n*n)]);

   if m = n then
      ### 
      ### no multiple eigenvalues, 
      ### normalize all eigenvectors
      ### 
      for j from 1 to n do

         Lambda[j,j]:=op(1,op(j,s));  # get eigenvalue
         v:=op(1,op(3,op(j,s)));      # get eigenvector
         v:=linalg[normalize](v);     # normalize

         for i from 1 to n do         # put into X
            X[i,j]:=v[i]
         od;

      od;

   else
      ### 
      ### there are multiple eigenvalues
      ### need to orthogonalization
      ### 
      ll:=1;
      for j from 1 to m do

         k:=op(2,op(j,s));        # get multiplicity
         lambda:=op(1,op(j,s));   # get eigenvalue
         u:=op(3,op(j,s));        # get eigenvector set
         if k > 1 then
            u:=GramSchmidt(u);    # make eigenvectors orthogonal
         fi;

         for l from 1 to k do
            Lambda[ll, ll]:=lambda;
            v:=linalg[normalize](op(l,u)); 
            for i from 1 to n do
               X[i,ll]:=v[i]
            od;
            ll:=ll+1;
         od;
      od;

   fi;

   print(evalm(Lambda),evalm(X));

end;