function mgt(D,dprim,KJ,LJ1,ntry1,nvru,lab,IN,inel,ngen,map,nmberg,dfi,multab,G,steer)
%
% This function determines a unique eigencolumn which transforms
% irreducibly. The column is obtained as an eigencolumn of a matrix,
% constructed as a linear combination of matrix group representatives.
%
% Version 1.0  050908
%
fid = fopen('Irrep_res.m', 'a');
nvr = LJ1*dprim;
for N1 = 1:nvr
  fi(1:G,1) = dfi(N1,1:G).';
  IV = 1;
  I2 = 1;
  I1 = 1;
  while I1 <= G
    if ntry1(I2) ~= I1
%
% Operate with group element I1 on fi(1:G,1).
%
      K1 = inel(I1);
      vec(1:G) = fi(multab(K1,1:G),1).';
      J1 = 1;
      while J1 <= IV
        rin = vec(1:G)*conj(fi(1:G,J1));
        rinab = abs(rin);
        if abs(rinab-1) >= 0.000001
          if rinab >= 0.000001
            for J2 = 1:G
              vec(J2) = vec(J2) - rin*fi(J2,J1);
            end
            rnorm = sum((abs(vec(1:G).')).^2);
            if rnorm < 0.000001
              break;
            end
            rnorm2 = 1/sqrt(rnorm);
            vec(1:G) = rnorm2*vec(1:G);
          end
        else
          break;
        end
        J1 = J1 + 1;
      end
      if J1 > IV
        IV = IV + 1;
        fi(1:G,IV) = vec(1:G).';
        if IV >= nvr
          break;
        end
      end
    else
      if I2 < nvru
        I2 = I2 + 1;
      end
    end
    I1 = I1 + 1;
  end
  if I1 <= G & IV >= nvr
    dfin(1:nvr,1:G) = fi(1:G,1:nvr)';
    subm(1:nvr,1:nvr,1) = eye(nvr);
%
% A nvr*nvr dimensional representation subm(1:nvr,1:nvr,1:G) is produced,
% containing the KJth irreducible representation dprim times.
%
    for K1 = 2:nmberg
      K2 = ngen(K1);
      genfi(1:G,1:nvr) = fi(multab(inel(K2),1:G),1:nvr);
      subm(1:nvr,1:nvr,K2) = dfin(1:nvr,1:G)*genfi(1:G,1:nvr);
    end
    for J1 = 1:G
      ML2 = map(J1,2);
      if ML2 ~= 0
        ML1 = map(J1,1);
        subm(1:nvr,1:nvr,J1) = subm(1:nvr,1:nvr,ML1)*subm(1:nvr,1:nvr,ML2);
      end
    end
    K1 = 1;
    K2 = 1;
    for J1 = 1:G
      if ntry1(K1) == J1
        if K1 < nvru  
          K1 = K1 + 1;
        end  
      else
        ntry(K2) = J1;
        K2 = K2 + 1;
      end
    end
    nvrt = K2 - 1;
    L1 = D*dprim;
    L2 = nvr - L1;
    for I = 1:nvru
      if ntry1(I) == IN
        break;
      end
    end
    [zr,wr] = eig(subm(1:nvr,1:nvr,ntry1(I)));
    J = 1;
    for I = 1:nvr
      if abs(lab-wr(I,I)) < 0.000001
        ortr(1:nvr) = zr(1:nvr,I);
        zr(1:nvr,I) = zr(1:nvr,J);
        zr(1:nvr,J) = ortr(1:nvr).';
        J = J + 1;
      end
    end
    fi(1:nvr,1:nvr) = zr(1:nvr,1:nvr);
%
% The representation is transformed into a new one where the eigenvalue lab
% appears in the first D*dprim diagonal places in subm(I1,I2,IN).
%
    for I = 1:G
      genfi(1:nvr,1:nvr) = fi(1:nvr,1:nvr)'*subm(1:nvr,1:nvr,I)*fi(1:nvr,1:nvr);
      subm(1:nvr,1:nvr,I) = genfi(1:nvr,1:nvr);
    end
    grupel(1:nvr,1:nvr,1:G) = subm(1:nvr,1:nvr,1:G);
    ntry2 = ones(1,G);
    ntry3 = zeros(1,G);
%
% The indices of the entries to be put equal to zero are stored in ntry2(I) and
% ntry3(J), respectively.
%
    for I = 1:G
      J = I;
      if I > 2*L1*L2
        break;
      end
      if J <= L1*L2
        while J > L2
          J = J - L2;
          ntry2(I) = ntry2(I) + 1;
        end
        ntry3(I) = J + L1;
      else
        J = J - L1*L2;
        ntry2(I) = ntry2(I) + L1;
        while J > L1
          J = J - L1;
          ntry2(I) = ntry2(I) + 1;
        end
        ntry3(I) = J;
      end
    end
%
% A block diagonal matrix with a D*dprim-dimensional block to split the
% degeneracy is constructed by linear combination of matrix
% representatives.
%
    J1 = 1;
    I1 = 1;
    while I1 <= 2*L1*L2
      J2 = J1;
      while J2 <= nvrt
        if abs(subm(ntry2(I1),ntry3(I1),ntry(J2))) >= 0.000001
          break;
        end
        J2 = J2 + 1;
      end
      if J2 > nvrt
        I1 = I1 + 1;
      else
        if J2 < nvrt
          J3 = J2 + 1;
          while J3 <= nvrt
            if abs(subm(ntry2(I1),ntry3(I1),ntry(J3))) >= 0.000001
              break;
            end
            J3 = J3 + 1;
          end
          if J3 <= nvrt
            rin = subm(ntry2(I1),ntry3(I1),ntry(J2));
            for J4 = (J2+1):nvrt
              LAB = subm(ntry2(I1),ntry3(I1),ntry(J4));
              subm(1:nvr,1:nvr,ntry(J4)) = LAB*subm(1:nvr,1:nvr,ntry(J2)) -...
                                       rin*subm(1:nvr,1:nvr,ntry(J4));
            end
          end
          I1 = I1 + 1;
          J1 = J2;
        else
          break;
        end
      end
    end
    if I1 <= 2*L1*L2
      continue
    end
    I3 = J1;
    while I3 <= G
      I1 = 1;
      while I1 <= L1
        I2 = 1;
        while I2 <= L1
          if abs(subm(I1,I2,ntry(I3))) > 0.000001
            break;
          end
          I2 = I2 + 1;
        end
        if I2 <= L1
          break;
        end
        I1 = I1 + 1;
      end
      if I1 <= L1
        break;
      end
      I1 = L1 + 1;
      while I1 <= nvr
        I2 = L1 + 1;
        while I2 <= nvr
          if abs(subm(I1,I2,ntry(I3))) > 0.000001
            break;
          end
          I2 = I2 + 1;
        end
        if I2 <= nvr
          break;
        end
        I1 = I1 + 1;
      end
      if I1 <= nvr
        break;
      end
      I3 = I3 + 1;
    end
    J1 = I3;
    [zr,wr] = eig(subm(1:L1,1:L1,ntry(J1)));
    J3 = 1;
    while J3 <= L1
      J2 = 1;  
      while J2 <= J3
        if (abs(wr(J3,J3)-wr(J2,J2)) >= 0.000001)
          break;
        end
        J2 = J2 + 1;
      end
      if J2 <= J3
        break;
      end
      J3 = J3 + 1;
    end
    if J3 > L1
      fprintf(fid, 'Error1, J3 > L1 in mgt\n\n');
      continue
    end
    for I1 = 1:nvr
      for I2 = 1:nvr
        if (I2 > L1) | (I1 > L1)
          fi(I2,I1) = 0;
          if I1 == I2
            fi(I2,I2) = 1;
          end
        else
          fi(I2,I1) = zr(I2,I1);
        end
      end
    end
    for J7 = 1:nvru
      if map(ntry1(J7),2) ~= 0
        grupel(1:nvr,1:nvr,ntry1(J7)) = subm(1:nvr,1:nvr,ntry(J1));
        break;
      end
    end
%
% All matrix representatives are transformed by means of the eigencolumns
% of subm(I1,I2,ntry(J1)).
%
    for I = 1:G
      genfi(1:nvr,1:nvr) = fi(1:nvr,1:nvr)'*grupel(1:nvr,1:nvr,I)*...
                           fi(1:nvr,1:nvr);
      subm(1:nvr,1:nvr,I) = genfi(1:nvr,1:nvr);
    end
    [zr,wr] = eig(subm(1:L1,1:L1,ntry1(J7)));
    J3 = 1;
    while J3 <= L1
      for J2 = 1:J3
        if (abs(wr(J3,J3)-wr(J2,J2)) >= 0.000001)
          break;
        end
      end
      if (abs(wr(J3,J3)-wr(J2,J2)) >= 0.000001)
        break;
      end
      J3 = J3 + 1;
    end
    if J3 > L1
      fprintf(fid, 'Error2, J3 > L1 in mgt\n\n');
      continue
    end
    I = L1 + 1;
    zr((L1+1):nvr,J2) = zeros((nvr-L1),1);
    fi(1:nvr,1) = zr(1:nvr,J2);
%
% The irreducible matrix representatives are calculated and stored in
% grupel(I1,I2,1:G).
%
    IV = 1;
    K3 = 1;
    K11 = 1;
    while K3 <= G
      if ntry1(K3) == 0
        vec(1:nvr) = subm(1:nvr,1:nvr,K3)*fi(1:nvr,1);
        while K11 <= IV
          rin = conj(fi(1:nvr,K11))*vec(1:nvr);
          rinab = abs(rin);
          if abs(rinab-1) >= 0.000001
            if rinab >= 0.000001
              for K4 = 1:nvr
                vec(K4) = vec(K4) - rin*fi(K4,K11);
              end
              rnorm = sum(abs(vec(1:nvr).').^2);
              if rnorm < 0.000001
                break;
              end
              rnorm2 = 1/sqrt(rnorm);
              vec(1:nvr) = vec(1:nvr)*rnorm2;
            end
          else
            break;
          end
          K11 = K11 + 1;
        end
        if K11 > IV
          IV = IV + 1;
          fi(1:nvr,IV) = vec(1:nvr).';
        end
        if IV >= LJ1
          break;
        end
      end
      K3 = K3 + 1;
    end
    if K3 > G
      continue
    end
    dfin(1:LJ1,1:nvr) = fi(1:nvr,1:LJ1)';
    grupel(1:LJ1,1:LJ1,1) = eye(LJ1);
    for K1 = 2:nmberg
      grupel(1:LJ1,1:LJ1,ngen(K1)) = dfin(1:LJ1,1:nvr)*...
                                     subm(1:nvr,1:nvr,ngen(K1))*fi(1:nvr,1:LJ1);
    end
    for LG = 1:G
      ML2 = map(LG,2);
      if ML2 ~= 0
        ML1 = map(LG,1);
        grupel(1:LJ1,1:LJ1,LG) = grupel(1:LJ1,1:LJ1,ML1)*...
                                 grupel(1:LJ1,1:LJ1,ML2);
      end
    end
    if steer(9)
      fprintf(fid, '\n\n');
      for I3 = 1:G
        for I4 = 1:LJ1
          for I5 = 1:LJ1
            grupelord(1,I5,I4,I3) = real(grupel(I5,I4,I3));
            grupelord(2,I5,I4,I3) = imag(grupel(I5,I4,I3));
          end
        end
      end
      for I3 = 1:G
        fprintf(fid, 'Representation number %3d, element number %3d\n', KJ, I3);
        for I4 = 1:LJ1
          fprintf(fid, '(%7.4f %7.4f) (%7.4f %7.4f) (%7.4f %7.4f) (%7.4f %7.4f)\n',...
                  grupelord(1:2,I4,1:LJ1,I3));
          fprintf(fid, '\n');
        end
      end
      return;
    end
  end
end
fprintf(fid, 'No irreducible representation for representation number %3d\n', KJ);
fclose(fid);
return