# csdms-contrib/slepian_juliet

Switch branches/tags
Nothing to show
Fetching contributors…
Cannot retrieve contributors at this time
65 lines (60 sloc) 2.07 KB
 function tracecheck(L,A,mv,tol,meth) % TRACECHECK(L,A,mv,tol,meth) % % Checks the trace (sum of the eigenvalues), or the sum of the squares of % the eigenvalues, of a matrix quadratic form L'*A*L, against an % analytically derived expression. No output for pass, error for fail. % % INPUT: % % L The cell array with the various first matrices % A The cell array with the various second matrices % mv This which should be equal to twice the negative trace or % which should be equal to the sum of the squares of the eigenvalues % tol The negative of the tolerance exponent % meth 1 Checks the sum of the eigenvalues % 2 Checks the sum of the squares of the eigenvalues % % SEE ALSO: % % MAOSL % % Last modified by fjsimons-at-alum.mit.edu, 06/20/2018 for j=1:length(A) % Checks a random wavenumber randi=ceil(rand*length(A{j})); % These are matrices out of MAOSL mrand=mv{j}(min(randi,length(mv{j}))); Arand=[A{j}(randi,1) A{j}(randi,2) ; A{j}(randi,2) A{j}(randi,3)]; % This is a lower-triangular matrix out of a Cholesky decomposition Lrand=[L(randi,1) 0 ; L(randi,2) L(randi,3)]; switch meth case 1 % Eq. (122) in doi: 10.1093/gji/ggt056 chekit=trace(Lrand'*Arand*Lrand); chekot=-2*mrand; case 2 chekit=sum([eig(Lrand'*Arand*Lrand)].^2); chekot=mrand; otherwise error('Specify valid method') end if ~isnan(mrand) && ~isnan(chekit) % This may by chance be at zero wavenumber in which case we skip the % test. Check RELATIVE accuracy, do not return an output to DIFER! % Maybe should only check the leading terms for the tiny ones (D) % since that's where the finite precision is going to be pretty bad relac=(chekit-chekot)/chekot; difer(relac,tol,[],NaN) else if j==1 disp(sprintf('The %ist entity not being checked',j)) elseif j==2 disp(sprintf('The %ind entity not being checked',j)) elseif j==3 disp(sprintf('The %ird entity not being checked',j)) else disp(sprintf('The %ith entity not being checked',j)) end end end