Home > solver > DSA_BD_normOfOperator.m

DSA_BD_normOfOperator

PURPOSE ^

DSA_BD_normOfOperator: compute the first singular value of At

SYNOPSIS ^

function [L] = DSA_BD_normOfOperator(blk,At,tol,m)

DESCRIPTION ^

 DSA_BD_normOfOperator: compute the first singular value of At
%
% DSA_BD:
% C. Ortiz
% Last Modified 3/31/2011
%***********************************************************

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 % DSA_BD_normOfOperator: compute the first singular value of At
0002 %%
0003 %% DSA_BD:
0004 %% C. Ortiz
0005 %% Last Modified 3/31/2011
0006 %%***********************************************************
0007 
0008 function [L] = DSA_BD_normOfOperator(blk,At,tol,m)
0009 
0010 
0011 if nargin < 4
0012     tol = 1e-6;
0013 end
0014 if nargin < 5
0015     m = size(At{1},2);
0016 end
0017 
0018 fprintf('\n Calculating Norm of A');
0019 MAXITER = 2000;
0020 TOLERANCE = min(1e-6,tol);
0021 try
0022     AAt = sparse(m,m);
0023     for p = 1:size(blk,1)
0024         AAt = AAt + At{p,1}'*At{p,1};
0025     end
0026     pertdiag = 1e-13*ones(m,1);
0027     AAt = AAt + spdiags(pertdiag,0,m,m);
0028     options.disp = 0;
0029     options.maxit = MAXITER;
0030     options.tol = TOLERANCE;
0031     L = sqrt(svds(AAt,1,'L',options));
0032 
0033     fprintf('\nNorm of A = %d\n',L);
0034     return;
0035 catch e
0036     fprintf('\nsvds(A,1,L,options) is not used.' );
0037 end
0038 
0039 q = rand(size(At{1},2),1);%q = rand(size(A,1),1);
0040 q = q/norm(q);
0041 Atq = Atyfun1a(blk,At,q);%Atyfun3(At,q,n);
0042 L1 = 0;
0043 L2 = 1;
0044 q1 = q;
0045 
0046 for i = 1:MAXITER
0047     z =DSA_BD_AXfun(blk,At,Atq);% AXfun(blk,At,[],Atq);%AXfun3(A,Atq);
0048     q = z/norm(z);
0049     Atq = DSA_BD_Atyfun(blk,At,q);%Atyfun(blk,At,[],[],q);%Atyfun3(At,q,n);
0050     L1 = q'*DSA_BD_AXfun(blk,At,Atq);%AXfun(blk,At,[],Atq);%AXfun3(A,Atq);
0051     error1 = abs(L1 - L2);
0052     error2 = norm(q - q1);
0053     if i > 1 && error1 < TOLERANCE && error2 < TOLERANCE
0054         break;
0055     end
0056     q1 = q;
0057     L2 = L1;
0058 end
0059 L =sqrt(L1);
0060 if (i == MAXITER)
0061     fprintf('WARNING: Power iteration did not converged. Errors: %f %f',error1,error2);
0062 end
0063 fprintf('\nPower iterations = %d',i);
0064 fprintf('\nNorm of A = %d\n',L);
0065 end

Generated on Mon 19-Sep-2011 21:41:33 by m2html © 2005