CommonFunctions.m

Go to the documentation of this file.

%%   Eotvos Quantum Transport Utilities - CommonFunctions
%    Copyright (C) 2009-2015 Peter Rakyta, Ph.D.
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see http://www.gnu.org/licenses/.
%
%> @addtogroup basic Basic Functionalities
%> @{
%> @file CommonFunctions.m
%> @brief A class containing common basic functionalities.
%> @} 
%>  @brief A class containing common basic functionalities.
%%
classdef CommonFunctions < handle


properties (Access=protected)
%> True if MKL component is built, false otherwise.
     DeployedMKL
%> The current version of the package.
    version = 'v4.9.146';
%> Name of the package.
    ProgramName = 'Eotvos Quantum Transport Utilities';
%> Short name of the package.
    ProgramShortName = 'EQuUs';
%> Maximal size of full matrixes to be handled
    MaxSize = 1e4;
end


methods


%% xmlread
%> @brief Function to load XML files (based on xerces libraries), providing octave compatibility.
%> @param filename Absolute path to the file to be opened. (In matlab relative path is sufficient)
%> @return The loaded document object model (see http://www.mathworks.com/help/matlab/ref/xmlread.html#outputarg_DOMnode for details.)
function tree = xmlread( obj, filename ) 
    
     if obj.isOctave()
        try
            parser = javaObject('org.apache.xerces.parsers.DOMParser');
        catch err
            error ('xmlread: could not load Xerces parser object, you may have to add required .jar files to your javapath. See the documentation for further details')
        end
        
        parser.parse(filename); 
        tree = parser.getDocument;
     else        
         tree = xmlread(filename); 
     end
    
end


%% xmlwrite
%> @brief Function to export XML files (based on xerces libraries), providing octave compatibility.
%> @param filename Absolute path to the file to be opened. (In matlab relative path is sufficient)
%> @param DOM The loaded document object model (see http://www.mathworks.com/help/matlab/ref/xmlread.html#outputarg_DOMnode for details.)
function xmlwrite( obj, filename, DOM )
    
    if obj.isOctave()
        try
            jfile = javaObject ('java.io.File', filename);
            jos = javaObject ('java.io.FileOutputStream', jfile);
            jfmt =  javaObject ('org.apache.xml.serialize.OutputFormat', DOM);            
        catch err
            error ('xmlwrite: could not load Xerces serializer object, you may have to add required .jar files to your javapath. See the documentation for further details')
        end
        
        jfmt.setIndenting (1);
        serializer = javaObject ('org.apache.xml.serialize.XMLSerializer', ...
                             jos, jfmt);
        serializer.serialize (DOM);
  
        if (nargout > 0 &&  strcmp (class (jos), 'java.io.StringWriter'))
            str = char (jos.toString ());
        end                                                  
                         
    else        
        xmlwrite(filename, DOM);
    end
    
end




%% IsDeployedMKL
%> @brief Checks whether the MKL component is deployed
%> @return Returns with true if MKL component is deployed, false otherwise.
function ret = IsDeployedMKL( obj )

    if isempty(obj.DeployedMKL)
          obj.DeployedMKL = obj.checkMEXfile( 'dgetPartialInv' ) & ...
                      obj.checkMEXfile( 'zgetPartialInv' );
    end

    ret = obj.DeployedMKL;

end



%% partialInv
%> @brief Calculates a partial inverse of a sparse matrix. If MKL component is not developed, the backslash operator is used insted.
%> @param A A sparse matrix to be inverted
%> @param sizeInv The size of partial inverse to be calculated.
%> @return Returns the calculated partial inverse
function invA = partialInv( obj, A, sizeInv )
    
    if ~issparse(A) 
       if ( size(A,1) > obj.MaxSize )
           error('EQuUs:CommonFunctions:partialInv', ['Matrix A has to be sparse or a full matrix smaller than ', num2str(obj.MaxSize), 'x', num2str(obj.MaxSize)]);
       end
    end

     if (~obj.IsDeployedMKL()) && issparse(A)
        evec = sparse( size(A,1)-sizeInv+1:size(A,1), 1:sizeInv, 1, size(A,1), sizeInv);
        invA = full(A \ evec);
        invA = invA( size(A,1)-sizeInv+1:end, :);
    elseif ~issparse(A)
        evec = zeros( size(A,1), sizeInv);
        evec( end-sizeInv+1:end, : ) = eye(sizeInv);
        invA = full(A \ evec);
        invA = invA( size(A,1)-sizeInv+1:end, :);        
     else
     try
          if isreal(A)
               invA = dgetPartialInv(A, sizeInv);
          else
               invA = zgetPartialInv(A, sizeInv);
          end
          catch errCause
               save('matrix.mat', 'A');
               error(errCause)
          end
     end

end


%% diagInv
%> @brief Calculates the diagonal part of the inverse of a sparse matrix. If MKL component is advised to build.
%> @param A A sparse matrix to be inverted
%> @return Returns the diagonal elements of the inverse
function invA = diagInv( obj, A ) 
    
    if ~issparse(A) && ( size(A,1) > obj.MaxSize )
       error('EQuUs:CommonFunctions:diagInv', ['Matrix A has to be sparse or smaller full matrix than ', num2str(obj.MaxSize), 'x', num2str(obj.MaxSize)]);
    end    
    
     if ~obj.IsDeployedMKL() && (~issparse(A))
        invA = diag(inv(A));
    else
        
     try
          if isreal(A)
               invA = dgetDiagInv(A);
          else
               invA = zgetDiagInv(A);
          end
          catch errCause
               save('matrix.mat', 'A');
               throw(errCause)
          end
     end

end



%% eig
%> @brief Calculates the generalized eigenvalue problem based on the zggev and dggev LAPACK functions.
%> @param A A square matrix on the left side.
%> @param B A square matrix on the right side.
%> @return Returns the calculated generalized eigenvalues and the right and left sided eigenvectors.
function [eigenvecs_right, eigenvecs_left, eigvals] = eig( obj, A, B )
%TODO check whether LAPACK is deployed    
    % if both the matrices A and B are reals
    if isreal(A) && isreal(B)
        % check the presence of the compiled MEX file
        fncname = 'dggev';
        [directory, name, ext] = fileparts(which(fncname));
        if (~strcmpi(name, fncname) || ~strcmpi(ext, ['.' mexext]))
           warning('CommonFunction:eig', ...
               'Function dggev is not compiled yet. Please compile first with the MKL component');  
           [eigenvecs_right, eigenvecs_left, eigvals] = eig_MATLAB();
           return
        end
        
        [alphavarR, alphavarI, betavar, eigenvecs_left, eigenvecs_right] = dggev(A,B);
        eigenvecs_left = eigenvecs_left';
        eigvals = (alphavarR + 1i*alphavarI)./betavar;
        return
    end
    
    if isreal(A)
        A = complex(A);
    end
        
    if isreal(B)        
        B = complex(B);
    end
    
    % check the presence of the compiled MEX file 
    fncname = 'zggev';
    [directory, name, ext] = fileparts(which(fncname));
    if (~strcmpi(name, fncname) || ~strcmpi(ext, ['.' mexext]))
       warning('CommonFunction:eig', ...
               'Function zggev is not compiled yet. Please compile first with the MKL component');    
       [eigenvecs_right, eigenvecs_left, eigvals] = eig_MATLAB();
       return
    end
        
    [alphavar, betavar, eigenvecs_left, eigenvecs_right] = zggev(A,B);
    eigenvecs_left = eigenvecs_left';
    eigvals = alphavar./betavar;
    
    % nested functions
    
    %--------------------------------
    % Calculates the right and left eigevectors with MATLAB
    function [eigenvecs_right, eigenvecs_left, eigvals] = eig_MATLAB()
        
        [modusmtx_right,k]  = eig(A,B);
        [modusmtx_left,k_left]  = eig(transpose(A),transpose(B));
        modusmtx_left = transpose(modusmtx_left);
        modusmtx_left_tmp = zeros(size(modusmtx_left));
        
        k = diag(k);
        k_left = diag(k_left);
        
        for idx = 1:length(k)
           [~,jdx] = min( abs(k_left -k(idx)));
           modusmtx_left_tmp(idx,:) = modusmtx_left(jdx,:);           
        end
        
        eigenvecs_right = modusmtx_right;
        eigenvecs_left = modusmtx_left_tmp;
        eigvals = k;
        
    end
    
    % end nested functions
    
end


%% getVersion
%> @brief Gets the current version of the package.
%> @return Returns the current version of the package. 
function ret = getVersion( obj )
     ret = obj.version;
end

%% getProgramName
%> @brief Gets the name of the package.
%> @return Returns the name of the package.
function ret = getProgramName( obj )
     ret = obj.ProgramName;
end

%% getProgramShortName
%> @brief Gets the short name of the package.
%> @return Returns the short name of the package.
function ret = getProgramShortName( obj )
     ret = obj.ProgramShortName;
end



end


methods (Access = public, Static = true )
    
    
%> @brief Simphson integral: y is a vector of function values at equal distant points: y_i = f(x_i), x_i-x_{i-1}=h
%> @param y Function values to be integrated
%> @param h Increment between the x_i points.
    function ret = SimphsonInt(y,h)
        if ~exist( 'h', 'var' )
            h = 1;
        end
        
        if mod(length(y),2) == 0
            warning('length of the vector to integral has to be (2n+1)')
            y(end) = [];
        end
        
        ya = y(1);
        yb = y(end);
        y(1) = [];
        
        ret = h/3*( ya + yb ) + 2*h/3*sum(y(2:2:end-2)) + 4*h/3*sum(y(1:2:end-1));
        
        
    end
    

%% inv_SVD
%> @brief Inverts badly conditioned matrix A with SVD regularization
%> @param A A square matrix to be inverted.
%> @return Returns the calculated inverse.
function invA = inv_SVD( A )
    SVDtolerance = 1e-15;
    
    if issparse(A)
        [U,S,V] = svd(full(A));
    else
          [U,S,V] = svd(A);
    end 
    
    S = diag(S);
     smax = max(abs(S));
     indexes = logical( abs(S) < smax*SVDtolerance );
     S(indexes) = smax*SVDtolerance;
    
    S = 1./S;                
                
    invA = V*diag(S)*U';    

end

%% checkMEXfile  
%> @brief Checks the presence of a given MEX file.
%> @param fncname The name of the function.
%> @return Returns true if the MEX file is present, false otherwise.
function ret = checkMEXfile( fncname )
     [directory, name, ext] = fileparts(which(fncname));
     if (~strcmpi(name, fncname) || ~strcmpi(ext, ['.' mexext]))
          ret = false;
     else
          ret = true;
     end
end


%% isOctave
%> @brief Checks whether the running environment is matlab or octave.
%> @return Returns true if running environment is octave, false otherwise.
function retval = isOctave()
  persistent cacheval;  % speeds up repeated calls

  if isempty (cacheval)
    cacheval = (exist ('OCTAVE_VERSION', 'builtin') > 0);
  end

  retval = cacheval;
end




end %methods


end % Global end