You are an expert in writing code in many different languages. Your goal is to perform code completion for the following code, keeping in mind the rest of the code and the file meta data.
    Metadata details:
    filename: spm_ness.m
    path of file: ./repos/spm
    
    the code of the file until where you have to start completion: function [H,R,J,G] = spm_ness(J,G)
% Evaluation of hessian and solenoidal operators at NESS
% FORMAT [H,R]     = spm_ness(J,G)
% FORMAT [H,R,J,G] = spm_ness(J,G)   %%%% complex
% J  - Jacobian (dfdx)
% G  - diffusion tensor (amplitude of random fluctuations)
%
% H  - Hessian matrix (i.e., precision of a Gaussian density)
% R  - Skew symmetric solenoidal operator (-Q')
%
% if called with four output arguments, complex forms are returned
%__________________________________________________________________________
% This routine evaluates the Hessian (i.e., precision) of a nonequilibrium
% steady-state density (using a local linear approximation, under Gaussian
% assumptions). This is evaluated  under linear constraints on the
% solenoidal term of a Helmholtz decomposition. In short, given the flow
% (encoded by the systems Jacobian) and amplitude of random fluctuations,
% one can evaluate the steady-state density under nonequilibrium dynamics
% implied by solenoidal flow.
%
% There are additional notes using symbolic maths and numerical examples in
% the main body of the script.
%
% flow constraints (Jacobian J)(R = -Q')
%--------------------------------------------------------------------------
% where flow   f = (R + G)*d log(p(x))/dx and
% log(p(x))      = -(1/2)*x'*H*x =>
% d log(p(x))/dx = -H*x =>
% df/dx = J      = -(R + G)*H =>
% H              = -(R + G)\J =>
% J*R + R*J'     = J*G - G*J'
%__________________________________________________________________________

% Karl Friston
% Copyright (C) 2008-2022 Wellcome Centre for Human Neuroimaging


% solve for solenoidal (R) operator J*R + R*J' = J*G - G*J'
%==========================================================================
if nargout < 3
    n  = size(J,1);
    L  = speye(n*n)/512;
    I  = speye(n,n);
    X  = kron(I,J) + kron(conj(J),I);
    Y  = spm_vec(J*G - G*J');
    R  = reshape((X'*X + L)\(X'*Y),n,n);
    
    % precision (inverse covariance) of steady-state density
    %----------------------------------------------------------------------
    H  = -(R + G)\J;
    
else
    
    % complex form
    %======================================================================
    
    % eigenbasis
    %----------------------------------------------------------------------
    E  = eig(J,'nobalance','vector');
    J  = pinv(E)*J*E;
    G  = pinv(E)*G*E;
    
    % solenoidal and Hessian
    %----------------------------------------------------------------------

    
function [H,R,J,G] = spm_ness(J,G)
% Evaluation of hessian and solenoidal operators at NESS
% FORMAT [H,R]     = spm_ness(J,G)
% FORMAT [H,R,J,G] = spm_ness(J,G)   %%%% complex
% J  - Jacobian (dfdx)
% G  - diffusion tensor (amplitude of random fluctuations)
%
% H  - Hessian matrix (i.e., precision of a Gaussian density)
% R  - Skew symmetric solenoidal operator (-Q')
%
% if called with four output arguments, complex forms are returned
%__________________________________________________________________________
% This routine evaluates the Hessian (i.e., precision) of a nonequilibrium
% steady-state density (using a local linear approximation, under Gaussian
% assumptions). This is evaluated  under linear constraints on the
% solenoidal term of a Helmholtz decomposition. In short, given the flow
% (encoded by the systems Jacobian) and amplitude of random fluctuations,
% one can evaluate the steady-state density under nonequilibrium dynamics
% implied by solenoidal flow.
%
% There are additional notes using symbolic maths and numerical examples in
% the main body of the script.
%
% flow constraints (Jacobian J)(R = -Q')
%--------------------------------------------------------------------------
% where flow   f = (R + G)*d log(p(x))/dx and
% log(p(x))      = -(1/2)*x'*H*x =>
% d log(p(x))/dx = -H*x =>
% df/dx = J      = -(R + G)*H =>
% H              = -(R + G)\J =>
% J*R + R*J'     = J*G - G*J'
%__________________________________________________________________________

% Karl Friston
% Copyright (C) 2008-2022 Wellcome Centre for Human Neuroimaging


% solve for solenoidal (R) operator J*R + R*J' = J*G - G*J'
%==========================================================================
if nargout < 3
    n  = size(J,1);
    L  = speye(n*n)/512;
    I  = speye(n,n);
    X  = kron(I,J) + kron(conj(J),I);
    Y  = spm_vec(J*G - G*J');
    R  = reshape((X'*X + L)\(X'*Y),n,n);
    
    % precision (inverse covariance) of steady-state density
    %----------------------------------------------------------------------
    H  = -(R + G)\J;
    
else
    
    % complex form
    %======================================================================
    
    % eigenbasis
    %----------------------------------------------------------------------
    E  = eig(J,'nobalance','vector');
    J  = pinv(E)*J*E;
    G  = pinv(E)*G*E;
    
    % solenoidal and Hessian
    %----------------------------------------------------------------------
    R  = 1j*real(G)*imag(J)/real(J);
    H  = - real(J)/G;
    
end