function  [nu,nnull,MatReta_d,Rmu,MatRVector] = sub_PCA_function(n_d,N_d,MatRx_d)

    %----------------------------------------------------------------------------------------------------------------------------------------------------
    % sub_PCA_function : PCA REDUCTION OF the random vector X = (X_1,...,X_{n_d}) for which N_d samples are available in MatRx_d(n_d,N_d)
    %
    % INPUT
    %      n_d                    = dimension of X
    %      N_d                    = number of realizations in the database
    %      MatRx_d(n_d,N_d)       = N_d realizations of X = (X_1, ..., X_{n_d}) 
    %
    % OUTPUT
    %      nu                     = order of the PCA reduction
    %      nnull = n_d - nu       = dimension of the null space  
    %      MatReta_d(nu,N_d)      = N_d samples of H = (H_1,...,H_nu)  
    %      Rmu(nu,1)              = vector of eigenvalues in descending order
    %      MatRVector(n_d,nu)     = matrix of the eigenvectors
    %----------------------------------------------------------------------------------------------------------------------------------------------------

      error_PCA = 1e-6;                                      % relative error

      iexec_plot_PCA = 1;                                    % iexec_plot_PCA = 1: plot the relative error 
      iexec_debugg   = 0;                                    % iexec_debugg = 1;   debugg

      RXmean        = mean(MatRx_d,2);                       % RXmean(n_d,1)
      MatRXmean     = repmat(RXmean,1,N_d);                  % MatRXmean(n_d,N_d)
      MatRtemp      = (MatRx_d - MatRXmean).^2;              % traceMatRXcov = trace(MatRXcov);  
      Rtemp         = sum(MatRtemp,2)/(N_d-1);
      clear MatRtemp
      traceMatRXcov = sum(Rtemp); 
      clear Rtemp

      if n_d <= N_d                                          %--- solving with an eigenvalue problem
         MatRXcov    = cov(MatRx_d');                        %    MatRXcov(n_d,n_d)
         MatRXcov    = 0.5*(MatRXcov+MatRXcov');     
         [MatRVectorTemp,MatRmuTemp] = eig(MatRXcov);        %    MatRVectorTemp(n_d,n_d),MatRmuTemp(n_d,n_d)
         RmuTemp     = diag(MatRmuTemp);
         [Rmu,Index] = sort(RmuTemp,'descend');              %    Rmu(n_d,1)
         MatRVector  = MatRVectorTemp(:,Index);              %    MatRVector(n_d,n_d)                       
                                                             %    computing nu and the related matrices
         RerrPCA        = 1-cumsum(Rmu,1)/traceMatRXcov;     %    RerrPCA(n_d,1)
         Ind = find(RerrPCA <= error_PCA);
                                                             %    Adapting the dimension of Rmu(nu,1) and MatRVector(n_d,nu)
         if isempty(Ind) == 0                                %    nu < n_d
            nu = Ind(1);
            if nu+1 <= n_d
               Rmu(nu+1:n_d) = [];                           %    Rmu(nu,1)
               MatRVector(:,nu+1:n_d)   = [];                %    MatRVector(n_d,nu)
            end
         else                                                %    nu = n_d
            nu = n_d;
         end  
         nnull = n_d - nu;
      end

      if n_d > N_d                                                        %--- solving with a "thin SVD" without assembling the covariance matrix of X
         [MatRVectorTemp,MatRSigma,~] = svd(MatRx_d-MatRXmean ,'econ');   % MatRVectorTemp(n_d,N_d), MatRSigma(N_d,N_d) 
         [RSigma,Index] = sort(diag(MatRSigma),'descend');                % Rsigma(N_d,1)
         Rmu        = (RSigma.^2)/(N_d-1);                                % Rmu(N_d,1)          
         MatRVector = MatRVectorTemp(:,Index);                            % MatRVector(n_d,N_d)  
         RerrPCA    = 1-cumsum(Rmu,1)/traceMatRXcov;                      % RerrPCA(N_d,1)
         Ind = find(RerrPCA <= error_PCA);
         if isempty(Ind) == 0                                             % nu < N_d
            nu = Ind(1);
            if nu+1 <= N_d
               Rmu(nu+1:N_d) = [];                                        % Rmu(nu,1)
               MatRVector(:,nu+1:N_d)   = [];                             % MatRVector(n_d,nu)
            end
         else
            nu = N_d;
         end
         nnull = n_d - nu;
      end

      if iexec_plot_PCA == 1                                           %--- plot the eigenvalues and the relative error 

         h = figure;                                                                                                                      
         plot((1:1:nu)',Rmu(1:nu,1),'-')                                                 
         title('Graph of the eigenvalues for the PCA of scaled X for generating additional realizations','FontSize',16,'Interpreter','tex','FontWeight','normal')                                         
         xlabel('j','FontSize',16,'Interpreter','tex')                                                                
         ylabel('\mu(j)','FontSize',16,'Interpreter','tex')
         saveas(h,'figureB1_Rmu.fig'); 
     %    saveas(h,'figureB1_Rmu.eps','epsc');
		 hold off
		 close(h);
        

         h = figure;                                                                                                                      
         plot((1:1:nu)',RerrPCA(1:nu,1),'-')                                                 
         title('Graph of the error function err_{PCA} for the PCA of scaled X for generating additional realizations','FontSize',16,'Interpreter','tex','FontWeight','normal')                                         
         xlabel('j','FontSize',16,'Interpreter','tex')                                                                
         ylabel('err_{PCA}(j)','FontSize',16,'Interpreter','tex')                                                           
         saveas(h,'figureB1_RerrPCA.fig'); 
     %    saveas(h,'figureB1_RerrPCA.eps','epsc'); 
         close(h);		 
      end
                                                                       %--- Construction of  the samples of RH = (H_1,...,H_nu)
                                                                       %    MatReta_d(nu,N_d)
      Rcoef     = 1./sqrt(Rmu);
      MatRcoef  = diag(Rcoef);
      MatReta_d = MatRcoef*MatRVector'*(MatRx_d - MatRXmean);

      if iexec_debugg == 1                                             %--- checking the projection
         RHmean       = mean(MatReta_d,2);
         MatReta_dcov = cov(MatReta_d'); 
         MatRP_nu     = MatRXmean + MatRVector*(diag(sqrt(Rmu)))*MatReta_d;  %  MatRVector(n_d,nu), Rmu(nu,1), MatReta_d(nu,N_d)
         MatRXcov     = cov(MatRx_d'); 
         norm(MatRXcov-cov(MatRP_nu'),'fro')
      end  
      return
end