Vandermonde rectangular, symeig

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);

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));
 

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

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

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

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

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

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

Stabilization/adfschur.m

@@ -1,4 +1,4 @@
-function [outputArg1,outputArg2] = adfschur(n,nel,nu)
+function [] = adfschur(n,nel,nu)
 
 [D,x,w]=legD(n);
 

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

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);

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;

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;

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

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))<k))
             ff=find(itopo(:,ee)==e);
             itopo(ff,ee)=ee;
             iflux(ff,ee)=0;

Visualization/setlatex.m

@@ -1,6 +1,7 @@
-function [] = setlatex()
+function [] = setlatex(fs)
 % Sets latex as default text interpreter
-set(groot, 'DefaultAxesFontSize', 14);
+if(nargin<1), fs=14; end
+set(groot, 'DefaultAxesFontSize', fs);
 set(groot, 'DefaultTextInterpreter',          'latex');
 set(groot, 'DefaultLegendInterpreter',        'latex');
 set(groot, 'DefaultAxesTickLabelInterpreter', 'latex');