function [Sple,theta_true]=GenerateData_ggm(p,s,n)
	%p number of variable; s sparsity rate
	epsilon=1;
	%%%%First we generate the "true" precision matrix theta_true
    A = sprandn(p,p,0.5*(s/p^2));
    B=A+A'; 
    ind=B>0; B(ind)=B(ind)+3; %magnify the non-zeros
    ind=B<0; B(ind)=B(ind)-3;
    Bfun=@(x)(B*x);
    %find the spectral range of B
    [min_eig, max_eig] = EigK(Bfun, p);
    min_eig = min_eig(length(min_eig));
    theta_true = full(B + (epsilon-min_eig)*speye(p));
	
	%Sple=GenerateGaussian3(theta_true,n);
	L_th=chol(theta_true,'lower'); Sple= (L_th')\randn(p,n); %using Cholesky
    Sple = Sple';
end



function [min_eig_hist, max_eig_hist, iter_hist, convg] = EigK(A, n, opts)
%function [min_eig_hist, max_eig_hist, iter_hist, convg] = EigK(A, n, opts)
%
% EigK computes the smallest and the largest eigenvelue of an n*n
% Hermitian operator A. The input A can either be a function handle or
% a matrix. The optional input 'opts' specifies the parameters of the
% Lanczos algorithm. By default,
%
%   opts.maxiter  1000      : Maximum number of Lanczos steps
%   opts.tol      sqrt(eps) : Tolerance for convergence
%   opts.t_step   10        : Freqeuency for checking convergence
%
% On return,
%
%   min_eig_hist : History of the smallest eigenvalue. The last
%                  element gives the final computed value.
%   max_eig_hist : History of the largest eigenvalue. The last element
%                  gives the final computed value.
%   iter_hist    : At which iterations are the values in
%                  min_eig_hist and max_eig_hist recorded?
%   convg        : Whether Lanczos has converged
%
% This routine is a quick hack for easy eigenvalue problems, where
% efficiency is substantially improved. For robust calculation, use
% eigs.
%

if ~exist('opts','var')      opts = struct;             end
if ~isfield(opts,'maxiter')  opts.maxiter = 1000;       end
if ~isfield(opts,'tol')      opts.tol     = sqrt(eps);  end
if ~isfield(opts,'t_step')   opts.t_step  = 10;         end

maxiter = opts.maxiter;
tol     = opts.tol;
t_step  = opts.t_step;

beta    = zeros(maxiter+1,1);
alpha   = zeros(maxiter,1);
Q       = zeros(n,maxiter+1);
T       = zeros(maxiter+1);

Q(:,1)  = randn(n,1);
Q(:,1)  = Q(:,1)/norm(Q(:,1));

if isnumeric(A)
  A = @(x)(A*x);
end

i = 1;
min_eig_hist = zeros(maxiter,1);
max_eig_hist = zeros(maxiter,1);
iter_hist = zeros(maxiter,1);
inc = 0;
while i <= maxiter
  
  % lanczos step
  if i > 1
    w = A(Q(:,i)) - beta(i)*Q(:,i-1);
  else
    w = A(Q(:,i));
  end
  alpha(i) = w'*Q(:,i);
  w = w - alpha(i)*Q(:,i);
  
  % full re-orthogonalization
  w = w - Q(:,1:i-1)*(w'*Q(:,1:i-1))';

  % continue lanczos step
  beta(i+1) = norm(w);
  if beta(i+1) < 10*eps
    T(i,i) = alpha(i);
    eig_val = eig(T(1:i,1:i));
    inc = inc + 1;
    min_eig_hist(inc) = min(eig_val);
    max_eig_hist(inc) = max(eig_val);
    iter_hist(inc) = i;
    min_eig_hist(inc+1:end) = [];
    max_eig_hist(inc+1:end) = [];
    iter_hist(inc+1:end) = [];
    convg = 1;
    return;
  end
  Q(:,i+1) = w/beta(i+1);
  
  T(i,i) = alpha(i);
  T(i,i+1) = beta(i+1);
  T(i+1,i) = beta(i+1);
  
  % convergence test
  if (mod(i,t_step)==0) || (i==maxiter)
    eig_val = eig(T(1:i,1:i));
    inc = inc + 1;
    min_eig_hist(inc) = min(eig_val);
    max_eig_hist(inc) = max(eig_val);
    iter_hist(inc) = i;
    if inc > 1 & ...
          abs(min_eig_hist(inc)-min_eig_hist(inc-1))/abs(min_eig_hist(inc-1)) < tol & ...
          abs(max_eig_hist(inc)-max_eig_hist(inc-1))/abs(max_eig_hist(inc-1)) < tol
      min_eig_hist(inc+1:end) = [];
      max_eig_hist(inc+1:end) = [];
      iter_hist(inc+1:end) = [];
      convg = 1;
      return;
    end
  end
  
  % increase i
  i = i + 1;
  
end

% if ever reaching here, Lanczos has not converged
fprintf(1, 'EigK: Eigenvalue computation has not converged. You are advised to plot semilogy(iter_hist,min_eig_hist) and semilogy(iter_hist,max_eig_hist) and see whether you have set a too small opts.maxiter or a too small opts.tol.\n\n');
min_eig_hist(inc+1:end) = [];
max_eig_hist(inc+1:end) = [];
iter_hist(inc+1:end) = [];
convg = 0;
end

