%% gram_schmidt: function description
function outputs = gram_schmidt(arg)
  c1 = arg([1 : rows(arg)],1);
  c2 = arg([1 : rows(arg)],2);

  res1 = c1 / sqrt(dot(c1, c1));

  res2 = c2 - dot(c2, c1) / dot(c1, c1) * c1;
  res2 = res2 /sqrt(dot(res2,res2));

  outputs =transpose([transpose(res1); transpose(res2)]);
end

function [Z,B]=dwaww (mat,max_iter_count =1000,TOL=0.001)
    if(norm(mat-transpose(mat))!=0)
      error("not a symmetric matrix")
    endif

    iter_count = 0;
    columns_count = columns(mat);
    Z = eye(columns_count,2);

    while(true)
      if(iter_count > max_iter_count)
        error('exceeded the maximal iteration limit')
      endif
      iter_count = iter_count+1;

      W = mat*Z;
      B = transpose(Z) * mat * Z;
      if (abs(B(2, 1)) < TOL)
        return
      endif
      Z = gram_schmidt(W);
      endwhile
end


%% test1: function description
function [] = test1(i)
  % for i=2:arg
    sim_matrix=randi(20,i,i);
    sim_matrix=sim_matrix-tril(sim_matrix,-1)+triu(sim_matrix,1)';
    [z,b]=dwaww(sim_matrix)

    all_eigenvalues= eig(sim_matrix)
    eigenvalues=-sort(-abs(eig(sim_matrix)))([1:2])
    my_eigenvalues= -sort(abs([b(1,1),b(2,2)]))
    if(norm(eigenvalues-my_eigenvalues)>0.1)
      error('eigenvalues are not equal');
    endif
  % endfor
 end