classdef bslv_test < meq_test
  % test of bslv (block linear problem solving algorithm)
  %
  % [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
  
  properties
    a1,b1,as,bs,A,B,b0;
    tola = 1e-2;
    tolb = 1e-2;
    nit = 50;
  end
  
  properties(MethodSetupParameter)
    nr = {2}; % For now these are just random integers
    ne = {2};
    ni = {3};
    nj = {4};
    na = {2};   
  end
  
  methods(TestMethodSetup)
    function bslvTestSetup(testCase,nr,ne,ni,nj,na)
   
      Wre = rand(nr,nj);
      Wrj = rand(nr,na);
      Wee = rand(ne,nj);
      Wij = rand(ni,na);
      
      Yr = rand(nr,1);
      Ye = rand(ne,1);
      Yi = rand(ni,1);
      
      Ajj = Wrj'*Wrj + Wij'*Wij;
      Aee = Wre'*Wre + Wee'*Wee;
      
      testCase.a1 = inv(Ajj)*(Wrj'*Yr + Wij'*Yi);
      testCase.b1 = (inv(Aee)*Wre')*Yr + (inv(Aee)*Wee')*Ye;
      
      testCase.A = inv(Ajj)*Wrj'*Wre;
      testCase.B = (inv(Aee)*Wre')*Wrj;
      
      % True solution
      s = [Wre Wrj; Wee zeros(ne,na); zeros(ni,nj) Wij] \ [Yr;Ye;Yi];
      testCase.bs = s(1:nj);
      testCase.as = s(nj+1:end);
      
      % Generate initial guess close to the real solution
      testCase.b0 = s(1:nj)+0.1*rand(nj,1).*s(1:nj);
      
    end
  end
  
  methods(Test, TestTags={'Unit'})
    function bslvTest(testCase)
      
      [a,b] = bslvmex(testCase.a1,testCase.b1,testCase.b0,testCase.A,testCase.B,testCase.nit);
      
      b2 = testCase.b0;
      for k = 1:testCase.nit
        a2 = testCase.a1 - testCase.A*b2;
        b2 = testCase.b1 - testCase.B*a2;
      end

      testCase.verifyLessThan(norm(a-testCase.as),testCase.tola,'test result failure')
      testCase.verifyLessThan(norm(a2-a),testCase.tola,'test result failure')
      testCase.verifyLessThan(norm(b-testCase.bs),testCase.tola,'test result failure')
      testCase.verifyLessThan(norm(b2-b),testCase.tola,'test result failure')
      
    end
  end
  
end

