
%  Generalized_Eigen_sparse.m  ver 2.4  by Tom Irvine

function[fn,omega,ModeShapes,MST]=Generalized_Eigen_sparse(K,M,dof,ijk)
%
    nsz=size(K);

    n=nsz(1);

    fn=zeros(n,1);
    ModeShapes=zeros(n,n);


    sparsityThreshold = 0.5;
    nnz_M = nnz(M);
    nnz_K = nnz(K);
    total_M = numel(M);
    total_K = numel(K);

    fprintf('M sparsity: %.2f%%\n', 100 * (1 - nnz_M / total_M));
    fprintf('K sparsity: %.2f%%\n', 100 * (1 - nnz_K / total_K));

    isSparseM = (nnz_M / total_M) < sparsityThreshold;
    isSparseK = (nnz_K / total_K) < sparsityThreshold;

    if isSparseM
        fprintf('Converting M to sparse format.\n');
        M = sparse(M);
    else
        fprintf('M is treated as dense.\n');
    end

    if isSparseK
        fprintf('Converting K to sparse format.\n');
        K = sparse(K);
    else
        fprintf('K is treated as dense.\n');
    end    


    [V, D] = eigs(K, M, dof, 'smallestabs'); 


    lambda = diag(D);  % Extract eigenvalues

    % Sort eigenvalues and corresponding eigenvectors in ascending order
    [lambdaSorted, idx] = sort(lambda, 'ascend');
    V = V(:, idx);

    % Mass normalization: Vn^T * M * Vn = I
    for i = 1:size(V,2)
        normFactor = sqrt(V(:,i)' * M * V(:,i));
        if normFactor ~= 0
            if sum(V(:,i)) < 0
                V(:,i) = -V(:,i);  % Consistent sign convention
            end
            V(:,i) = V(:,i) / normFactor;
        end
    end

    % Return sorted and normalized results
    lambda = lambdaSorted;
    Vn = V;


    aModeShapes=Vn;

    ndof=length(lambda);

    aomega=zeros(ndof,1);
    
    for i=1:ndof
        aomega(i) = sqrt(lambda(i,1));
    end
    
    afn=aomega/(2*pi);


   fn(1:dof,1)=afn;
   omega=2*pi*fn;

   ModeShapes(:,1:dof)=aModeShapes;
   MST=ModeShapes';

    
    %
    if(ijk==1 || ijk==3)
	    disp(' Natural Frequencies ');
        disp(' No.      f(Hz)');
    
            
        for i=1:dof
            if(ijk==1 || ijk==3)
                
                if(imag(fn)==0)
                    out1=sprintf('%d.  %10.5g ',i,fn(i));
                else
                    out1=sprintf('%d.  %10.5g +%10.5gi ',i,real(fn(i)),imag(fn(i)));                
                end
                
                disp(out1);
            end
        end    
    
    end
    %

    
    if(ijk==1)
       disp(' ');
       disp('  Modes Shapes (column format)');
       ModeShapes
    end