%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2020-12-16 15:40:22
 % @LastEditTime: 2020-12-16 15:40:30
 
 % I = eye(3);
 % PP = u*(F + F'-2*I) + l * trace(F-I) * I;
%--------------------------
function P = PK1_fast_linear(u, l, F)
    m = size(F,2);
    
    FT  = F([1:3:9,2:3:9,3:3:9], :); % [9,m]
    trF = F(1,:) + F(5,:) + F(9,:) - 3; % 1*m

    I  = eye(3);
    I1 = repmat(I(:), 1, m); % 9*m

    P = u(:)' .* (F + FT - 2*I1) + l(:)' .* trF .* I1;
end