I assume that you are interested in the eigenvalues (and possibly the eigen vectors) of a SYMMETRIC matrix. You may write all the equations that define the eigenvectors, but most likely the NLP solver will not be able to solve the model because it is highly non linear and non-convex. However, you may compute the eigenvalues one by one by maximizing a norm. The following little GAMS program should show you the principle:

set i dimension of the symmetric matrix / i1 * i6 /,
    io(i) subset with already defined eigenvectors;
alias (i,ip,ipp);
parameter c(i,i) matrix for which eigenvalues are computed
          ev(i,i) eigenvectors
          la(i) eigenvalues;
* Generate random positive semidefinite symmetric matrix, c.
* The rank will be less than min( card(i), card(o) ):
*
set o / o1 * o6 /;
parameter d(i,o) data observations;
d(i,o) = normal(0,1);
c(i,ip) = sum(o, d(i,o) * d(ip,o) );
display c;
*
* Make sure io and ev are defined before we enter the loop
*
io(i)    = no;
ev(i,io) = 0;
variable x(i) eigenvector
         lambda eigenvalue;
equation edef     definition of max eigenvalue
         norm     eigenvector has norm 1
         ortho(i) eigen vector is orthonormal to already defined vectors;
edef ..      sum( (i,ip), x(i) * c(i,ip) * x(ip) ) =e= lambda;
norm ..      sum( i, sqr( x(i) ) ) =e= 1;
ortho(io)..  sum( i, ev(i,io) * x(i) ) =e= 0;
model eigen / all /;
option limrow = 0, limcol = 0, solprint = off;
* Now generate the eigen vectors one by one starting from the largest
* eigenvalue. Each new eigen vector must be orthogonal to the previous
* ones and maximize the weighted norm.
* If maximizing is replaced by minimizing we get them in the reverse
* order.
loop( ipp,
* Initialize with random vector and normalize it so the only nonlinear
* constraint is satiesfied:
         x.l(i) = normal(0,1);
         x.l(i) = x.l(i) / sum( ip, sqr( x.l(ip) ) );
solve eigen using nlp maximizing lambda;
* Store the current eigen vector and eigen value and increment the
* set of vectors that future vectors must be orthogonal to:
         ev(i,ipp) = x.l(i);
         io(ipp) = yes;
         la(ipp) = lambda.l;
);
display ev, la;
*
* Test the accuracy of the solution:
* If dtest below is too large you may reduce the optimality tolerance.
* In CONOPT this is done with the line "set rtredg 1.d-10" in a conopt
* options file, conopt.opt.
*
parameter test(i,i) computed tr(ev) * diag(la) * ev -- should be equal to c
          dtest(i,i) error in test;
test(i,ip)  = sum( ipp, ev(i,ipp) * la(ipp) * ev(ip,ipp) );
dtest(i,ip) = test(i,ip) - c(i,ip);
display test, dtest;

Another approach submitted by Arne Drud: Essentially, the eigenvalues are found one at a time from the largest one using maximization. Each time an eigenvector has been found we require that the next eigenvectors are orthogonal. A little trick in the generation of initial values removes components parallel to already found eigenvectors and ensures that the initial point to each solve is feasible. All constraints except one are linear and the objective is quadratic, so the model is easy to solve.

$Title First Covariance Matrix and its Eigen-values and -vectors.
* Contributed by Arne Drud (adrud@arki.dk)
SETS
 j          variables                            /j1*j25/
 n          observations                         /n1*n40/

ALIAS (j,jj,jjj);

PARAMETER
 jn1(j,n)    observations of the variables
 meanj(j)    mean value of the observations
 covj(j,jj)  First covariance matrix;

* forming a covariance matrix
jn1(j,n)   = NORMAL(0,1);
meanj(j)   = SUM(n, jn1(j,n)) /CARD(n);
covj(j,jj) = SUM(n, (jn1(j,n) - meanj(j))*(jn1(jj,n) - meanj(jj)))/CARD(n);
DISPLAY covj;

scalar OK;

Parameter evecj(j,j)  Eigenvectors found so far;
Parameter evalj(j)    Eigenvalues found so far;

$batinclude eigen.inc j covj evecj evalj ok call1

display evecj, evalj;

Parameter testj(j,j) New covariance matrix
          difj(j,j)  Difference between original and new covariance;

testj(j,jj) = sum( jjj, evecj(j,jjj) * evalj(jjj) * evecj(jj,jjj) );
difj(j,jj)  = covj(j,jj) - testj(j,jj);
display difj;

$Title Second Covariance Matrix and its Eigen-values and -vectors.
SETS
 i          variables                            /i1*i15/

ALIAS (i,ii,iii);

PARAMETER
 in1(i,n)    observations of variables
 meani(i)    mean value of the observations
 covi(i,ii)  Second covariance matrix;

* forming a covariance matrix
in1(i,n)   = NORMAL(0,1);
meani(i)   = SUM(n, in1(i,n)) /CARD(n);
covi(i,ii) = SUM(n, (in1(i,n) - meani(i))*(in1(ii,n) - meani(ii)))/CARD(n);
DISPLAY covi;

Parameter eveci(i,i)  Eigenvectors found so far;
Parameter evali(i)    Eigenvalues found so far;

$batinclude eigen.inc i covi eveci evali ok call2

display eveci, evali;

Parameter testi(i,i) New covariance matrix
          difi(i,i)  Difference between original and new covariance;

testi(i,ii) = sum( iii, eveci(i,iii) * evali(iii) * eveci(ii,iii) );
difi(i,ii)  = covi(i,ii) - testi(i,ii);
display difi; 

$offlisting
$ontext
  Batinclude file eigen.inc

  The program computes the eigenvalues and eigen vector of a square
  symmetric matrix.

  Parameters (all required):
  1: The set over which the matrix is defined (input)
  2: The square matrix (input)
  3: The eigenvectors (output) declared as a square matrix over the set
     in argument 1. The second set indicates the different eigenvalues.
  4: The eigenvalues (output) declared as a one dimensional parameter
     over the set in argument 1.
  5: Return Code. A scalar into which a return code will be written.
     A +1 indicates that not all eigenvectors and eigenvalues were found.
  6: A name that can be used to make all identifiers unique for multiple
     calls
$offtext

alias (%1, %6__1, %6__2 );
%5 = 1;
*
*  Define a model that finds the largest eigen value and eigen vector:
*
set       %6__3(%1)  Subset already having an associated eigenvector-eigenvalue pair;
Scalar    %6__4      Norm of initial random vector;

%6__3(%1) = no;
%3(%1,%1) = 0;

variable %6__5(%1)   Eigenvector
         %6__6       Eigenvalue;

equation %6__7       Definition of Objective;
Equation %6__8(%1)   The new Eigenvector must be orthogonal to earlier Eigenvectors;
Equation %6__9       The norm of the eigen vector must be 1;

%6__7        .. %6__6 =e= sum((%1,%6__1), %6__5(%1)*%2(%1,%6__1)*%6__5(%6__1) );

%6__8(%6__3) .. sum( %1, %3(%1,%6__3) * %6__5(%1) ) =e= 0;

%6__9        .. sum(%1, sqr(%6__5(%1)) ) =E= 1;

Model %6__10 / %6__7, %6__8, %6__9 /;
%6__10.limrow=0;
%6__10.limcol=0;
%6__10.solprint=2;
%6__10.solvelink=1;

parameter %6__11(%1);

loop( %6__2$%5,
*
*  We initialize the eigenvector with random numbers. To prevent using
*  a lot of time in phase 1 we make sure the orthogonality constraints
*  are satisfied by subtracting a multiple of the already found
*  eigenvectors. We do not check that the remaining vector still is
*  nonzero.
*
   %6__5.l(%1)     = uniform(-1,1);
   %6__11(%6__3)   = sum(%1, %6__5.l(%1)*%3(%1,%6__3));
   %6__5.l(%1)     = %6__5.l(%1) - sum(%6__3,%6__11(%6__3)*%3(%1,%6__3));
   %6__4           = sqrt( sum(%1, sqr(%6__5.l(%1)) ) );
   %6__5.l(%1)     = %6__5.l(%1) / %6__4;
   solve %6__10 using nlp maximizing %6__6;
   %5$(%6__10.solvestat ne 1) = 0;
   %5$(%6__10.modelstat ne 1 and %6__10.modelstat ne 2) = 0;
   %6__3(%6__2)    = Yes;
   %3(%1,%6__2)    = %6__5.l(%1);
   %4(%6__2)       = %6__6.l;
);
*
*  Clear the content of all local variables to preserve space
*
option kill=%6__3, kill=%6__4, kill=%6__5, kill=%6__6, kill=%6__11;
*
*  Clear as much of the symbol table as possible.
*  It seems that aliasses cannot be cleared.
*
$kill %6__3
$kill %6__4
$kill %6__5
$kill %6__6
$kill %6__7
$kill %6__8
$kill %6__9
$kill %6__11