function Result = inverse4x4(M)
  % ========================
  % https://guiard.org/pages/octave/inverseX.html
  % =========================
  %
  % Parallel inverte multiple matrices of size [4 4].
  % size(A) = [4 4 x]
  % Version 1.00
  %

  A = M([1:2],[1:2],:);

  B = M([1:2],[3:4],:);

  C = M([3:4],[1:2],:);

  D = M([3:4],[3:4],:);


  function Mult = mult2x2(A,B)
    % Ermöglicht vektorisierte Matrizenmultiplikation von 2x2 Matrizen.
    % size(A) = size(B) = [2 2 x]

    Mult = zeros(size(A));

    Mult(1,1,:) = A(1,1,:) .* B(1,1,:) + A(1,2,:) .* B(2,1,:);
    Mult(1,2,:) = A(1,1,:) .* B(1,2,:) + A(1,2,:) .* B(2,2,:);
    Mult(2,1,:) = A(2,1,:) .* B(1,1,:) + A(2,2,:) .* B(2,1,:);
    Mult(2,2,:) = A(2,1,:) .* B(1,2,:) + A(2,2,:) .* B(2,2,:);

  endfunction



  InvA = inverse2x2(A);
  InvB = inverse2x2(B);
  InvC = inverse2x2(C);
  InvD = inverse2x2(D);


  % Explizites Verfahren zur Inversion von 4x4 Matrizen nach Wikipedia: https://en.wikipedia.org/wiki/Invertible_matrix#Inversion_of_4_%C3%97_4_matrices

  %rA = InvA + InvA * B * inv(D-C*InvA*B)*C*InvA
  rA = InvA + mult2x2(mult2x2(mult2x2(mult2x2(InvA , B) , inverse2x2(D-mult2x2(mult2x2(C,InvA),B))) ,C),InvA )   ;

  %rB =-InvA * B * inv(D-C*InvA*B)
  rB =-mult2x2(mult2x2(InvA,B), inverse2x2(D-mult2x2(mult2x2(C,InvA),B)));

  %rC =-inv(D-C*InvA*B)*C*InvA
  rC =-mult2x2(mult2x2(inverse2x2(D-mult2x2(mult2x2(C,InvA),B)),C),InvA);

  %rD = inv(D-C*InvA*B)
  rD = inverse2x2(D-mult2x2(mult2x2(C,InvA),B));


  Result = zeros(size(M));

  Result([1:2],[1:2],:) = rA;

  Result([1:2],[3:4],:) = rB;

  Result([3:4],[1:2],:) = rC;

  Result([3:4],[3:4],:) = rD;

endfunction
