clear all error_tol=1.0e-6;%1.0e-7; % termination tolerance (relative error) tol =1.0e-6; %1.0e-6% Tolerance for PCG %%%%%% Weight for Obj. c1=1e0; c2=1e-0; %%%%%%%%%%% Generate the prescribed eigendata %%%%%%%%%%%%%%% %%%% Generate k random eigenvalues with k/2 complex-value ones %%%%%%% n=1000; k=30; Lambda=zeros(k,k); alpha=-abs(randn(k,1)); s=floor(1/4*k); % Number of complex eigenvalues beta=randn(s,1); Lambda=diag(alpha); i=1; while i<=s Lambda(2*(i-1)+1:2*i,2*(i-1)+1:2*i)=[alpha(i) beta(i);-beta(i) alpha(i)]; i=i+1; end %%%% Generate the eigenvector matrix randomly %%% X=randn(n,k); [Q,R]=qr(X);% QR factorization of the eigenvectors X R=R(1:k,:); % R factor i=1; while i<=k if norm(R(i,i))<1.0e-6 R(i,i) =1.0e-6; end i=i+1; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%% Strictly feasible solution given by Bai, Chu and Sun %%%%%% %%%%% See Example 5.1 in our paper %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Produce a strictly feasible solution with the given eigendata %%%% Rinv=inv(R); Mtrue=[Rinv'*Rinv zeros(k,n-k);zeros(n-k,k) eye(n-k)]; Mtrue=c1^0.5*Mtrue; Ctrue=-[Rinv'*(Lambda+Lambda')*Rinv zeros(k,n-k);zeros(n-k,n)]; Ctrue=c2^0.5*Ctrue; Ktrue=Lambda*Rinv; Ktrue=[Ktrue'*Ktrue zeros(k,n-k);zeros(n-k,k) eye(n-k)]; W1=2*rand(n)-ones(n,n);W1=(W1'+W1)/2; W2=2*rand(n)-ones(n,n);W2=(W2'+W2)/2; W3=2*rand(n)-ones(n,n);W3=(W3'+W3)/2; %%%% Generate a perturbed solution (Ma,Ca,Ka) as the analytic solution %%%% tau=0.1; Ma=Mtrue+tau*W1; % generate the given symmetric Ma, Ca,and Ka Ca=Ctrue+tau*W2; Ka=Ktrue+tau*W3; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % %%%%%%%%%%%%%%%%%%%%%%% % Ma: the analytic mass matrix (n by n) % Ca: the analytic damping matrix (n by n) % Ka: the analytic stiffness matrix (n by n) % R: R factor of the given eigenvector matrix X (n by k) % Lambda: the given eigenvalue matrix (k by k) % c1 and c2: weighted parameters % error_tol: stopping crierion for the method: 1.0e-7 % tol: Tolerance for PCG % M,C,K: final solution % val_obj: final objective function value %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% tic; [M,C,K,y,z]= IQEP_Newton(Ma,Ca,Ka,R,Lambda,n,k,c1,c2);toc tic; [M1,C1,K1,y1,z1]= IQEP_Newton_old(Ma,Ca,Ka,R,Lambda,n,k,c1,c2);toc; %tic; [M,C,K,y,z]= IQEP_Newton_old(Ma,Ca,Ka,R,Lambda,n,k,c1,c2);toc %tic; [M,C,K,y,z]= IQEP_Newton(Ma,Ca,Ka,R,Lambda,n,k,c1,c2);toc