clear;clc;close all;
rng(0);


n = 4000;
p = 3;

fprintf('n:%d, p:%d\n', n, p);

h = 1 / (n + 1);
A = 1 / h^2 * spdiags(ones(n, 1) * [-1 2 -1], [-1 0 1], n, n);    
M = spdiags([rand(n-1, 1); 0] + 0.1, [0], n, n);
% M = speye(n);
LM = chol(M, 'lower');
b = randn(n, 1);
B = [];


Manifold = SymPosSemDefinMatFixRank(n, p);

preconditioner = 2;
problem = GenLyapProb(A,M,LM,B,b,n,p, preconditioner);

problem.M = Manifold;

% checkhessian(problem)
x0 = problem.M.rand();
gf0 = problem.grad(x0);
ngf0 = problem.M.norm(x0, gf0);

options.Delta_bar = 1e10;
options.Delta0 = 10;
options.tolgradnorm = ngf0 * 1e-10;
[x, cost, info, options] = trustregions(problem, x0, options);
sum([info.hessvecevals])

