function [A,gam,fval] = gam3ind(c1,c2,h,m)

global nstate lenseq
nstate = m;
%lenseq = n;

T = 1;

adim = m^2*T;
gdim = m^2;
xdim = adim + gdim;


% box constraints
%LB = zeros(adim,1)+1e-10;
%UB = ones(adim,1);
LB = zeros(xdim,1)+1e-10;
LB(adim+1:end) = LB(adim+1:end) - 1e9;  % no box constraints on phi
UB = ones(xdim,1);
UB(adim+1:end) = UB(adim+1:end) + 1e9;  % no box constraints on phi


% normalization constraints
Aeq = zeros(m*T,xdim);
for i=1:m*T
  for j=1:m
    Aeq(i,m*(i-1)+j) = 1;
  end
end
Beq = ones(m*T,1);


[Aina,Bina] = constrict(m,h,T);

rind=0;
Aing = zeros(0,gdim);
% gamma sum constraint
%for i=1:m
%  for s=[-1 1]
%    rind = rind+1;
%    for j=1:m
%      cind = coord2ind([i j],[m m]);
%      Aing(rind,cind) = s;
%    end
%    Bing(rind,1) = c1;
%  end
%end

mm = repmat([m],[1 m]);
minds = list_all_ind(mm);
for s=[-1 1]
  for jind = 1:size(minds,1)
    rind = rind+1;
    iseq = minds(jind,:);
    rho = ~all(iseq==iseq(1));
    for j=1:m
      i = iseq(j);
      cind = coord2ind([i j],[m m]);
      Aing(rind,cind) = s;
      Bing(rind,1) = c1+c2*rho;
    end
  end
end

Ain = assembleblock({Aina,Aing});
Bin = [Bina;Bing];

A0 = pnormdim(ones(m,m,T),1);
A0 = A0+randn(size(A0))*1e-2;
A0 = pnormdim(A0,1);

G0 = randn(m,m);
%G0(:,1) = G0(:,1) + 1e6;
%G0(:,2) = G0(:,2) - 1e6;

x0 = [A0(:); G0(:)];

[x,fval] = fmincon(@myobjm,x0,Ain,Bin,Aeq,Beq,LB,UB);
fval = sqrt(abs(fval));

y = myobjm(x);

global A FN

A = reshape(x(1:adim),[m m T]);
gam = reshape(x(adim+1:end),[m m]);
return


function y = myobjm(x)
global nstate lenseq

m = nstate;
%n = lenseq;
T = 1;

adim = m^2*T;
gdim = m^2;
xdim = adim + gdim;

A = reshape(x(1:adim),[m m T]);
gam = reshape(x(adim+1:end),[m m]);

s = sum(sum( A.*gam));

y = -s^2;
return