From aac3867c315079391bd3d965280e8f6779ef075d Mon Sep 17 00:00:00 2001 From: pbrubeck Date: Mon, 30 Nov 2020 18:10:21 +0000 Subject: [PATCH] Vandermonde rectangular, symeig --- DD/SchurNKP/curvedquad.m | 1 + DD/SchurNKP/meshtopo.m | 2 +- Linear Algebra/normc.m | 2 +- Linear Algebra/symeig.m | 36 +++++++++++++++++++++++++++ Linear Algebra/trideigs.m | 19 ++++++++++---- Root Finding/NewtonRaphsonOpt.m | 2 +- Spectral/VandermondeLeg.m | 10 +++++--- Spectral/legD.m | 44 ++++++++++++++++++++------------- Stabilization/adfschur.m | 2 +- Stabilization/box_topo.m | 25 +++++++++++++------ Stabilization/bubfilt.m | 2 +- Stabilization/color_greedy.m | 4 ++- Stabilization/crs2.m | 1 - Stabilization/get_graph.m | 7 ++---- Stabilization/toposort_loops.m | 7 +++--- Visualization/setlatex.m | 5 ++-- 16 files changed, 118 insertions(+), 51 deletions(-) create mode 100644 Linear Algebra/symeig.m diff --git a/DD/SchurNKP/curvedquad.m b/DD/SchurNKP/curvedquad.m index 77bcc2b..f622a98 100644 --- a/DD/SchurNKP/curvedquad.m +++ b/DD/SchurNKP/curvedquad.m @@ -6,6 +6,7 @@ c=(b+a)/2; d=(b-a)/2; h=sqrt(r.^2-abs(d).^2); + w=atan2(sign(r).*abs(d), h); s=1i*sign(r.*d); f=cell(4,1); diff --git a/DD/SchurNKP/meshtopo.m b/DD/SchurNKP/meshtopo.m index 032a3f5..e8821a5 100644 --- a/DD/SchurNKP/meshtopo.m +++ b/DD/SchurNKP/meshtopo.m @@ -115,7 +115,7 @@ % Add attached face to the vertex verts_faces{vert} = add_face_to_vertex_faces(verts_faces{vert}, [corner, face_idx]); - % Add the neighboring vertices to the vertix + % Add the neighboring vertices to the vertex verts_vertices{vert}= add_vertex_to_vertex_vertices(verts_vertices{vert}, face_verts(bitxor(corner-1,1)+1)); verts_vertices{vert}= add_vertex_to_vertex_vertices(verts_vertices{vert}, face_verts(bitxor(corner-1,2)+1)); diff --git a/Linear Algebra/normc.m b/Linear Algebra/normc.m index 9d1454d..92f3806 100644 --- a/Linear Algebra/normc.m +++ b/Linear Algebra/normc.m @@ -5,7 +5,7 @@ a = sqrt(G(:)'*(real(A).^2+imag(A).^2)); B=bsxfun(@rdivide, A, a); else - a = sqrt(diag(A'*G*A)).'; + a = sqrt(sum(A.*(G*A)),2).'; B=bsxfun(@rdivide, A, a); end else diff --git a/Linear Algebra/symeig.m b/Linear Algebra/symeig.m new file mode 100644 index 0000000..28b9306 --- /dev/null +++ b/Linear Algebra/symeig.m @@ -0,0 +1,36 @@ +function [varargout] = symeig(A,B) +if(nargout<2) + JOBZ='N'; +else + JOBZ='V'; +end +UPLO='L'; +N=size(A,1); +LWORK=max(20,N)*N; +L=zeros(N, 1); +WORK=zeros(LWORK, 1); +INFO=0; +if(nargin==1) + out=lapack('dsyev',JOBZ,UPLO,N,A,N,L,WORK,LWORK,INFO); + INFO=out{9}; + if(JOBZ=='V') + varargout{1}=out{4}; % V + varargout{2}=out{6}; % L + else + varargout{1}=out{6}; % L + end +else + ITYPE=1; + out=lapack('dsygv',ITYPE,JOBZ,UPLO,N,A,N,B,N,L,WORK,LWORK,INFO); + INFO=out{12}; + if(JOBZ=='V') + varargout{1}=out{5}; % V + varargout{2}=out{9}; % L + else + varargout{1}=out{9}; % L + end +end +if(INFO) + warning('dsygv::INFO = %d',INFO); +end +end \ No newline at end of file diff --git a/Linear Algebra/trideigs.m b/Linear Algebra/trideigs.m index a47839b..b5b89d0 100644 --- a/Linear Algebra/trideigs.m +++ b/Linear Algebra/trideigs.m @@ -1,6 +1,5 @@ function [W, Z] = trideigs(D, E) % Calculates the eigenvalues and eigenvectors of a symmetric tridiagonal matrix. - N=numel(D); IL=1; IU=N; M=N; LDZ=N; @@ -13,8 +12,18 @@ Z=zeros(LDZ, M); WORK=zeros(1, LWORK); -out=lapack('dstevr','V','A',N,D,E,VL,VU,IL,IU,ABSTOL,M,W,Z,LDZ,ISUPPZ, ... - WORK,LWORK,IWORK,LIWORK,INFO); -W=out{12}; -Z=out{13}; +if(nargout==2) + out=lapack('dstevr','V','A',N,D,E,VL,VU,IL,IU,ABSTOL,M,W,Z,LDZ,ISUPPZ, ... + WORK,LWORK,IWORK,LIWORK,INFO); + W=out{12}; + Z=out{13}; +else +% out=lapack('dstevr','N','A',N,D,E,VL,VU,IL,IU,ABSTOL,M,W,Z,LDZ,ISUPPZ, ... +% WORK,LWORK,IWORK,LIWORK,INFO); +% W=out{12}; +% Z=out{13}; + out=lapack('dsteqr','N',N,D,E,Z,LDZ,WORK,INFO); + W=out{3}; + Z=out{5}; +end end \ No newline at end of file diff --git a/Root Finding/NewtonRaphsonOpt.m b/Root Finding/NewtonRaphsonOpt.m index 200cda8..e404014 100755 --- a/Root Finding/NewtonRaphsonOpt.m +++ b/Root Finding/NewtonRaphsonOpt.m @@ -9,7 +9,7 @@ vars=num2cell(x); y=g(vars{:}); dx=H(vars{:})\y; - err=norm(dx'*y); + err=norm(y); x=x-dx; i=i+1; end diff --git a/Spectral/VandermondeLeg.m b/Spectral/VandermondeLeg.m index 763ea1e..5568760 100644 --- a/Spectral/VandermondeLeg.m +++ b/Spectral/VandermondeLeg.m @@ -1,6 +1,8 @@ -function [ V ] = VandermondeLeg(x) -V=zeros(length(x)); -for j=1:length(x) - V(:,j)=sqrt((2*j-1)/2)*LegendreP([zeros(j-1,1);1], x); +function [ V ] = VandermondeLeg(x,m) +if(nargin<2), m=length(x); end +V=zeros(length(x),m); +for j=1:m + V(:,j) = LegendreP([zeros(j-1,1);1], x); + V(:,j) = sqrt((2*j-1)/2) * V(:,j); end end \ No newline at end of file diff --git a/Spectral/legD.m b/Spectral/legD.m index e7ed8fb..21eecda 100644 --- a/Spectral/legD.m +++ b/Spectral/legD.m @@ -1,23 +1,33 @@ -function [D, x, w] = legD(N) +function [D, x, w] = legD(N,itype) % Legendre spectral differentiation matrix -x=zeros(1,N); x([1,N])=[-1,1]; -w=zeros(1,N); w([1,N])=2/(N*(N-1)); -if(N==2) - D=[-1,1;-1,1]/2; - x=x(:); w=w(:); - return; -end +% CG (GLL nodes) if itype==0 (default) +% DG (GL nodes) if otherwise +if(nargin<2), itype=0; end -k=1:N-3; -E=sqrt((k.*(k+2))./((2*k+1).*(2*k+3))); -[x(2:N-1),V]=trideigs(zeros(1,N-2), E); -w(2:N-1)=4/3*V(1,:).^2./(1-x(2:N-1).^2); +if(itype==0) + x=zeros(1,N); x([1,N])=[-1,1]; + w=zeros(1,N); w([1,N])=2/(N*(N-1)); + if(N==2) + D=[-1,1;-1,1]/2; + x=x(:); w=w(:); + return; + end + k=1:N-3; + E=sqrt((k.*(k+2))./((2*k+1).*(2*k+3))); + [x(2:N-1),V]=trideigs(zeros(1,N-2), E); + w(2:N-1)=4/3*V(1,:).^2./(1-x(2:N-1).^2); + %p=[1; ((-1).^(1:N-2)./abs(V(1,:)))'.*sqrt((1-x(2:N-1).^2)*3/(2*N*(N-1))); -(-1)^N]; +else + k=1:N-1; + E=k./sqrt(4*k.*k-1); + [x,V]=trideigs(zeros(1,N), E); + w=2*V(1,:).^2; +end x=x(:); w=w(:); - -X=repmat(x, [1, N]); -p=[1; ((-1).^(1:N-2)./abs(V(1,:)))'.*sqrt((1-x(2:N-1).^2)*3/(2*N*(N-1))); -(-1)^N]; -dX=X-X'+eye(N); -D=(p*(1./p)')./dX; +D=repmat(x, [1, N]); +D=D-D'+eye(N); +p=prod(D,2); +D=(p*(1./p)')./D; D=D-diag(sum(D,2)); end \ No newline at end of file diff --git a/Stabilization/adfschur.m b/Stabilization/adfschur.m index e1c93e6..74ead3a 100644 --- a/Stabilization/adfschur.m +++ b/Stabilization/adfschur.m @@ -1,4 +1,4 @@ -function [outputArg1,outputArg2] = adfschur(n,nel,nu) +function [] = adfschur(n,nel,nu) [D,x,w]=legD(n); diff --git a/Stabilization/box_topo.m b/Stabilization/box_topo.m index 0f38f46..23047e9 100644 --- a/Stabilization/box_topo.m +++ b/Stabilization/box_topo.m @@ -1,11 +1,20 @@ -function [itopo] = box_topo(nex,ney) -ndim=2; +function [itopo] = box_topo(nex,ney,nez) +ndim=nargin; +if(ndim==2) + nez=1; +end nfaces=2*ndim; -nel=nex*ney; -eid=reshape(1:nel,nex,ney); +nel=nex*ney*nez; itopo=zeros(nfaces,nel); -tt=eid; tt(2:end ,:)=eid(1:end-1,:); itopo(1,:)=tt(:); -tt=eid; tt(1:end-1,:)=eid(2:end ,:); itopo(2,:)=tt(:); -tt=eid; tt(:,2:end )=eid(:,1:end-1); itopo(3,:)=tt(:); -tt=eid; tt(:,1:end-1)=eid(:,2:end ); itopo(4,:)=tt(:); +eid=reshape(1:nel,nex,ney,nez); +tt=eid; tt(2:end ,:,:)=eid(1:end-1,:,:); itopo(1,:)=tt(:); +tt=eid; tt(1:end-1,:,:)=eid(2:end ,:,:); itopo(2,:)=tt(:); +if(ndim>=2) +tt=eid; tt(:,2:end ,:)=eid(:,1:end-1,:); itopo(3,:)=tt(:); +tt=eid; tt(:,1:end-1,:)=eid(:,2:end ,:); itopo(4,:)=tt(:); +end +if(ndim>=3) +tt=eid; tt(:,:,2:end )=eid(:,:,1:end-1); itopo(5,:)=tt(:); +tt=eid; tt(:,:,1:end-1)=eid(:,:,2:end ); itopo(6,:)=tt(:); +end end \ No newline at end of file diff --git a/Stabilization/bubfilt.m b/Stabilization/bubfilt.m index 12a612b..1a79dda 100644 --- a/Stabilization/bubfilt.m +++ b/Stabilization/bubfilt.m @@ -11,7 +11,7 @@ sigma = ones(N,1); p = 1/16; p=0; cut = max(1,ceil(N*p)); -alpha = (5E-1/2)*max(1/2,1-p*cut); +alpha = 0.25; for k=1:cut j = N+1-k; sigma(j) = 1-alpha*(((cut+1-k)/cut)^2); diff --git a/Stabilization/color_greedy.m b/Stabilization/color_greedy.m index 6732a61..be77e9e 100644 --- a/Stabilization/color_greedy.m +++ b/Stabilization/color_greedy.m @@ -7,7 +7,9 @@ icolor=zeros(nel,1); for e=elems k=1; - while any(k==icolor(itopo(:,e))) + id=itopo(:,e); + id=id(id>0); + while any(k==icolor(id)) k=k+1; end icolor(e)=k; diff --git a/Stabilization/crs2.m b/Stabilization/crs2.m index 237b600..32667ba 100644 --- a/Stabilization/crs2.m +++ b/Stabilization/crs2.m @@ -9,7 +9,6 @@ tau=(h./v).*(1-1./peh); tau(peh<=1)=0; - nel=numel(hx); N=2; n=N*N; diff --git a/Stabilization/get_graph.m b/Stabilization/get_graph.m index d3e9cab..35337d0 100644 --- a/Stabilization/get_graph.m +++ b/Stabilization/get_graph.m @@ -23,11 +23,9 @@ heat(e)=ux*ux+uy*uy; end - %heat=-heat; - %heat(:)=0; iflux=sign(iflux).*(abs(iflux)>tol); iflux(itopo==repmat(1:nel,nfaces,1))=0; - + igrad=zeros(nfaces,nel); for e=1:nel for f=1:nfaces @@ -49,5 +47,4 @@ end end end -end - +end \ No newline at end of file diff --git a/Stabilization/toposort_loops.m b/Stabilization/toposort_loops.m index d50b4c6..7ee27e6 100644 --- a/Stabilization/toposort_loops.m +++ b/Stabilization/toposort_loops.m @@ -1,6 +1,5 @@ function [isweep,icolor] = toposort_loops(itopo,iflux) -ndim=2; -nfaces=2*ndim; +nfaces=size(itopo,1); nel=size(itopo,2); isweep=zeros(nel,1); @@ -27,6 +26,7 @@ icolor(e)=k; else % apply greedy coloring of root (only needed for DG) + % FIXME greedy coloring of all neighbors iroot=find(icolor==1); icolor=color_greedy(itopo,iroot); iroot=find(icolor==1); @@ -60,7 +60,8 @@ e=isweep(p); for f=1:nfaces ee=itopo(f,e); - if(iflux(f,e)>=0 && ee~=e) + if(iflux(f,e)>=0 && ee~=e && ... + all(icolor(itopo(:,ee))