mp3: done 1-3

mp2/Project_2_Maggioni_Claudio/ex2.m

@@ -1,4 +1,18 @@
-n = 5;
+A = createA(5);
+
+spy(A);
+matlab2tikz('../ex2_2_spy.tex');
+
+A = createA(100);
+
+figure;
+subplot(1,2,1)
+spy(A);
+subplot(1,2,2)
+spy(chol(A));
+matlab2tikz('../ex2_3_spy.tex');
+
+function A = createA(n)
 A = sparse(n, n);
 A(1,:) = 1;
 A(end,:) = 1;
@@ -7,15 +21,7 @@ A(:,end) = 1;
 
 % make sure diagonal is 0
 A(1,1) = 0;
-A(n,n) = 0;
+A(end,end) = 0;
 
 A = A + diag(n:(2*n-1));
-spy(A);
+end
-matlab2tikz('../ex2_2_spy.tex');
-
-figure;
-subplot(1,2,1)
-spy(A);
-subplot(1,2,2)
-spy(chol(A));
-matlab2tikz('../ex2_3_spy.tex');

mp2/project_2_Maggioni_Claudio.tex

@@ -122,7 +122,7 @@ to solve $Ax = b$
 for a given righthand-side vector would be problematic.}
 
 Here is the plot of \texttt{spy(A)} (on the left) and \texttt{chol(spy(A))} (on
-the right).
+the right) for $n = 100$.
 
 \centering{\input{ex2_3_spy.tex}}
 

mp3/Makefile

@@ -0,0 +1,15 @@
+filename=template
+
+pdf:
+	pdflatex ${filename}
+	#bibtex ${filename}
+	pdflatex ${filename}
+	pdflatex ${filename}
+	make clean
+
+read:
+	evince ${filename}.pdf &
+
+
+clean:
+	rm -f *.out *.log *.bbl *.blg *.aux ${filename}.log  ${filename}.ps ${filename}.aux ${filename}.out ${filename}.dvi ${filename}.bbl ${filename}.blg

mp3/Project_3_Maggioni_Claudio/datasets/Mesh_generation/grid3d.m

@@ -0,0 +1,31 @@
+function [A,xyz] = grid3d(k)
+% GRID3D : Generate 3-dimensional 7-point finite difference mesh.
+%
+% [A,xyz] = GRID3D(k) returns a k^3-by-k^3 symmetric positive definite 
+%        matrix A with the structure of the k-by-k-by-k 7-point grid,
+%        and an array xyz of coordinates for the grid points.
+
+
+a  = blockdiags ([-1 6 -1], -1:1, k, k);
+I  = speye (k, k);
+aa = blockdiags ([-I a -I], -1:1, k, k);
+II = speye(k^2,k^2);
+A  = blockdiags ([-II aa -II], -1:1, k, k);
+
+A  = diag(diag(A)) - A;
+
+
+xyz = zeros(k^3,3);
+x   = ones(k,1) * (1:k);
+y   = x';
+
+j   = 1;
+for i = 1:k
+    xyz(j:j+k^2-1 , 1) = x(:);
+    xyz(j:j+k^2-1 , 2) = y(:);
+    xyz(j:j+k^2-1 , 3) = i * ones(k^2,1);
+    j = j + k^2;
+end
+
+
+end

mp3/Project_3_Maggioni_Claudio/datasets/Mesh_generation/grid3dt.m

@@ -0,0 +1,38 @@
+function [A,xyz] = grid3dt(k)
+% GRID3DT : Generate 3-dimensional tetrahedral finite element mesh.
+%
+% [A,xyz] = GRID3DT(k) returns a k^3-by-k^3 symmetric positive definite 
+%        matrix A of the k-by-k-by-k grid, with cells divided into tetrahedra, 
+%        and an array xyz of coordinates for the grid points.
+
+% 1-d mesh
+a = blockdiags ([-1 14 -1], -1:1, k, k);
+
+% glue to laminate 1-d meshes into 2-d meshes
+b = blockdiags ([-1 -1], [0 1], k, k);
+
+% 2-d mesh
+aa = blockdiags ([b a b'], -1:1, k, k);
+
+% glue to laminate 2-d meshes into 3-d meshes
+bb = blockdiags ([b b'], [0 1], k, k);
+
+% 3-d mesh
+A = blockdiags ([bb aa bb'], -1:1, k, k);
+A = diag(diag(A)) - A;
+
+
+% xyz coordinates of nodes
+xyz = zeros(k^3,3);
+x = ones(k,1) * (1:k);
+y = x';
+j = 1;
+for i = 1:k
+    slab = j:j+k^2-1;
+    xyz(slab,1) = x(:);
+    xyz(slab,2) = y(:);
+    xyz(slab,3) = i * ones(k^2,1);
+    j = j + k^2;
+end
+
+end

mp3/Project_3_Maggioni_Claudio/datasets/Mesh_generation/grid5.m

@@ -0,0 +1,21 @@
+function [A,xy] = grid5(k)
+% GRID5 : Generate 5-point finite difference mesh.
+%
+% [A,xy] = GRID5(k) returns a k^2-by-k^2 symmetric positive definite 
+%        matrix A with the structure of the k-by-k 5-point grid,
+%        and an array xy of coordinates for the grid points.
+
+a = blockdiags ([-1 4 -1], -1:1, k, k);
+I = speye (k, k);
+A = blockdiags ([-I a -I], -1:1, k, k);
+A  = diag(diag(A)) - A;
+
+
+xy = zeros(k^2,2);
+x  = ones(k,1) * (1:k);
+y  = x';
+xy(:,1) = x(:);
+xy(:,2) = y(:);
+
+
+end

mp3/Project_3_Maggioni_Claudio/datasets/Mesh_generation/grid5rec.m

@@ -0,0 +1,20 @@
+function [A,xy] = grid5rec(k1, k2)
+% GRID5lrec : Generate 5-point finite difference on a rectangular mesh.
+%
+% [A,xy] = GRID5REC(k) returns a k1*k2-by-k1*k2 symmetric positive definite 
+%        matrix A with the structure of the k-by-k 5-point grid,
+%        and an array xy of coordinates for the grid points.
+
+a = blockdiags ([-1 4 -1], -1:1, k1, k1);
+I = speye (k1, k1);
+A = blockdiags ([-I a -I], -1:1, k2, k2);
+A = diag(diag(A)) - A;
+
+
+[x,y] = meshgrid(1:k2, 1:k1);
+xy    = zeros(k1*k2,2);
+xy(:,1) = x(:);
+xy(:,2) = y(:);
+
+
+end

mp3/Project_3_Maggioni_Claudio/datasets/Mesh_generation/grid5recRotate.m

@@ -0,0 +1,25 @@
+function [A,xy] = grid5recRotate(k1, k2, angle)
+% GRID5lrec : Generate 5-point finite difference on a rectangular mesh.
+%
+% [A,xy] = GRID5REC(k) returns a k1*k2-by-k1*k2 symmetric positive definite 
+%        matrix A with the structure of the k-by-k 5-point grid,
+%        and an array xy of coordinates for the grid points.
+
+a = blockdiags ([-1 4 -1], -1:1, k1, k1);
+I = speye (k1, k1);
+A = blockdiags ([-I a -I], -1:1, k2, k2);
+A = diag(diag(A)) - A;
+
+
+[x,y] = meshgrid(1:k2, 1:k1);
+
+rad = angle*pi/180;
+
+newX = x * cos(rad) - y * sin(rad);
+newY = y * cos(rad) + x * sin(rad);
+
+xy    = zeros(k1*k2,2);
+xy(:,1) = newX(:);
+xy(:,2) = newY(:);
+
+end

mp3/Project_3_Maggioni_Claudio/datasets/Mesh_generation/grid7.m

@@ -0,0 +1,21 @@
+function [A,xy] = grid7(k)
+% GRID7 : Generate 7-point finite difference mesh.
+%
+% [A,xy] = GRID7(k) returns a k^2-by-k^2 symmetric positive definite 
+%        matrix A with the structure of the k-by-k 7-point grid,
+%        and an array xy of coordinates for the grid points.
+
+a = blockdiags ([-1 6 -1], -1:1, k, k);
+b = blockdiags ([-1 -1], [0 1], k, k);
+A = blockdiags ([b a b'], -1:1, k, k);
+A = diag(diag(A)) - A;
+
+
+xy = zeros(k^2,2);
+x  = ones(k,1) * (1:k);
+y  = x';
+xy(:,1) = x(:);
+xy(:,2) = y(:);
+
+
+end

mp3/Project_3_Maggioni_Claudio/datasets/Mesh_generation/grid9.m

@@ -0,0 +1,23 @@
+function [A,xy] = grid9(k)
+% GRID9 : Generate 9-point finite difference mesh.
+%
+% [A,xy] = GRID9(k) returns a k^2-by-k^2 symmetric positive definite 
+%        matrix A with the structure of the k-by-k 9-point grid,
+%        and an array xy of coordinates for the grid points.
+
+
+a = blockdiags ([-4 20 -4], -1:1, k, k);
+b = blockdiags ([-1 -4 -1], -1:1, k, k);
+A = blockdiags ([b a b], -1:1, k, k);
+A = diag(diag(A)) - A;
+
+xy = zeros(k^2,2);
+x  = ones(k,1) * (1:k);
+y  = x';
+
+
+xy(:,1) = x(:);
+xy(:,2) = y(:);
+
+
+end

mp3/Project_3_Maggioni_Claudio/datasets/Mesh_generation/gridt.m

@@ -0,0 +1,25 @@
+function [A,xy] = gridt(k)
+% GRIDT : Generate triangular mesh.
+%
+% [A,xy] = GRIDT(k) returns a k*(k+1)/2-square symmetric positive 
+% definite matrix A with the structure of an equilateral triangular 
+% mesh with side k, and an array xy of coordinates of the mesh points.
+
+
+% Start with a square 7-point grid.
+
+[A,xy] = grid7(k);
+
+% Mask off one triangle of it.
+
+xy = xy-1;
+f  = find(xy(:,1)+xy(:,2)<k);
+A  = A(f,f);
+xy = xy(f,:);
+
+% Make the other triangle equilateral.
+
+T = [ 1 0 ; 1/2 2/sqrt(5)];
+xy = xy*T;
+
+end

mp3/Project_3_Maggioni_Claudio/external/metismex.m

@@ -0,0 +1,30 @@
+function metismex
+% METISMEX Establish an interface between METIS and Matlab
+% 
+% [map,edgecut] = metismex('PartGraphRecursive',A,nparts,options);
+% [map,edgecut] = metismex('PartGraphKway',A,nparts,options);
+% [perm,iperm] = metismex('EdgeND',A,options);
+% [perm,iperm] = metismex('NodeND',A,options);
+%
+% options is now a structure with the following fields:
+%   'seed' : an integer for the random seed used in metis
+%   'ctype' : 'rm' or 'shem' [default]
+%   'iptype' : 'grow' or 'random' (only applys for recursive bisection)
+%   'objtype' : 'cut' or 'vol' [default] (only applys for part kway)
+%   'rtype' : '1sided' [default] or '2sided' 
+%   'ufactor' : an integer for balance
+%   'pfactor' : an integer for minimum degree of vertex to be ordered last
+%   'ccorder' : a flag (value ignored) to order connected components
+%               separately
+%   'nseps' : an integer for the number of separators tried at each level
+%               (default 1)
+%   'niter' : an integer for the number of refinement iterations  
+%               (default 10)
+%   'ncuts' : number of initial partitions to test (default 1)
+%   'dbglvl' : the debug level (default 0)
+% 
+% The output and options commands are optional.
+%
+% Note that error checking is not done: make sure A is structurally
+% symmetric or it will crash.
+

mp3/Project_3_Maggioni_Claudio/external/metispart.m

@@ -0,0 +1,31 @@
+function [p1,p2] = metispart(A,xy);
+% METISPART : Partition a graph using Metis default method.
+%
+% p = metispart(A) returns a list of the vertices on one side of
+%     a partition obtained by Metis 4.0 applied to the graph of A.
+%     
+% Optional arguments:
+%   metispart(A,xy) draws a picture of the partitioned graph,
+%                   using the rows of xy as vertex coordinates.
+%   [p1,p2] = metispart(...)   also returns the list of vertices
+%                              on the other side of the partition.
+%
+% See also METISMEX (which accepts all the Metis options), 
+% GEOPART, GSPART, SPECPART, METISDICE, METISND.
+%
+% John Gilbert  3 Jul 01
+% Copyright (c) 1990-2001 by Xerox Corporation.  All rights reserved.
+% HELP COPYRIGHT for complete copyright and licensing notice.
+
+if nargin < 2
+    xy = 0;
+end;
+picture = max(size(xy)) > 1;
+
+map = metismex('PartGraphRecursive',A,2);
+[p1,p2] = other(map);
+
+if picture
+    gplotpart(A,xy,p1);
+    title('Metis Partition')
+end;

mp3/Project_3_Maggioni_Claudio/external/other.m

@@ -0,0 +1,31 @@
+function [out1,out2] = other(in1,in2);
+% OTHER : Find the other part of a partition, or
+%         convert a partition to the other representation.
+%
+% part2 = other(part1,n)     : part1 is a list of subscripts in 1:n;
+%                              part2 becomes the complementary list.
+% part2 = other(part1,A)     : Same, except n is taken as max(size(A)).
+% [part2,p] = other(part1,n) : Also returns 0/1 partition vector of length n.
+% [part2,p] = other(part1,A) : Same.
+% [part1,part2] = other(p)   : Converts 0/1 vector to lists of subscripts.
+%
+% John Gilbert, 1993.
+% Copyright (c) 1990-1996 by Xerox Corporation.  All rights reserved.
+% HELP COPYRIGHT for complete copyright and licensing notice.
+
+% Modified 3 Jul 01 by JRG to fix a bug in the "map" case
+
+if nargin == 1
+
+    out1 = find(in1 == 0);
+    out2 = find(in1 ~= 0);
+
+else
+
+    if max(size(in2)) > 1
+        in2 = max(size(in2));
+    end;
+    out2 = full(spones(sparse(1,in1,1,1,in2)));
+    out1 = find(out2 == 0);
+
+end;

mp3/Project_3_Maggioni_Claudio/src/Bench_bisection.m

@@ -0,0 +1,196 @@
+function Bench_bisection()
+% Compare various graph bisection algorithms
+%
+% D.P & O.S for Numerical Computing at USI
+
+% add the necessaty paths
+addpaths_GP;
+
+warning('off','all');
+picture = 1;
+
+if nargin < 1
+    whichdemos = [1 2 3];
+end;
+
+format compact;
+
+disp('          *********************************************')
+disp('          ***      Graph bisection benchmark        ***');
+disp('          *********************************************')
+disp(' ');
+disp(' The file "Toy_meshes.mat" contains sample meshes with coordinates.');
+disp(' ');
+
+% load meshes
+load Toy_meshes;
+whos;
+
+
+for nmesh = 1:7
+    close all; clf reset;
+    
+    if (nmesh==1)
+        disp(' ');
+        disp(' Function "grid5rec" produces a rectangular grid:');
+        disp(' ');
+        disp('[A,xy] = grid5rec(8,80);');
+        disp(' ');
+        [W,coords] = grid5rec(8, 80);
+    end
+    if (nmesh==2)
+        disp(' ');
+        disp(' Function "grid5rec" produces a rectangular grid:');
+        disp(' ');
+        disp('[A,xy] = grid5rec(80,8);');
+        disp(' ');
+        [W,coords] = grid5rec(80, 8);
+    end
+    if (nmesh==3)
+        disp(' ');
+        disp(' Function "grid5recRotate" produces a rotated rectangular grid:');
+        disp(' ');
+        disp('[A,xy] = grid5recRotate(80,8, -45);');
+        disp(' ');
+        [W,coords] = grid5recRotate(80, 8, -45);
+    end
+    if (nmesh==4)
+        disp(' ');
+        disp(' Function "gridt" produces a triangular grid:');
+        disp(' ');
+        disp(' (See also grid5, grid7, grid9, grid3d, grid3dt.)');
+        disp(' ');
+        disp('[A,xy] = gridt(40);');
+        disp(' ');
+        [W,coords] = gridt(20);
+    end
+    if (nmesh==5)
+        disp(' ');
+        disp(' Function "gridt" produces a triangular grid:');
+        disp(' ');
+        disp(' (See also grid5, grid7, grid9, grid3d, grid3dt.)');
+        disp(' ');
+        disp('[A,xy] = gridt(30);');
+        disp(' ');
+        [W,coords] = grid9(30);
+    end
+    if (nmesh==6)
+        W       = Smallmesh;
+        coords = Smallmesh_coords;
+    end
+    if (nmesh==7)
+        W       = Tapir;
+        coords  = Tapir_coords;
+    end
+    if (nmesh==8)
+        W       = Eppstein;
+        coords  = Eppstein_coords;
+    end
+    
+    
+    disp(' ');
+    disp('          *********************************************')
+    disp('          ***        Various Bisection Methods      *** ');
+    disp('          *********************************************')
+    disp(' ');
+    disp(' ');
+    
+    
+    if (nmesh==1)
+        disp('An initial rectangular  grid5rec(8,80) mesh');
+    end
+    if (nmesh==2)
+        disp('An initial rectangular  grid5rec(80,8) mesh');
+    end
+    if (nmesh==3)
+        disp('An initial rectangular  grid5rec(80,8) mesh, rotated by 45 degree');
+    end
+    if (nmesh==4)
+        disp('  gridt(20) mesh');
+    end
+    if (nmesh==5)
+        disp('  gridt9(20) mesh');
+    end
+    if (nmesh==6)
+        disp(' Small mesh ');
+    end
+    if (nmesh==7)
+        disp(' "Tapir" is a test of a no-obtuse-angles mesh generation algorithm');
+        disp(' due to Bern, Mitchell, and Ruppert.  ');
+    end
+    if (nmesh==8)
+        disp('  Eppstein mesh');
+    end
+    
+    figure(1)
+    disp('gplotg(Tmesh,Tmeshxy);');
+    disp('  ');
+    gplotg(W,coords);
+    nvtx  = size(W,1);
+    nedge = (nnz(W)-nvtx)/2;
+    xlabel([int2str(nvtx) ' vertices, ' int2str(nedge) ' edges'],'visible','on');
+    
+    disp(' Hit space to continue ...');
+    pause;
+    
+    disp(' 1. Coordinate bisection of a mesh.  ');
+    disp(' p = coordpart(A,xy) returns a list of the vertices  ');
+    disp(' on one side of a partition obtained by bisection    ');
+    disp(' perpendicular to a coordinate axis.  We try every   ');
+    disp(' coordinate axis and return the best cut. Input W is ');
+    disp(' the adjacency matrix of the mesh; each row of xy is ');
+    disp(' the coordinates of a point in d-space.              ');
+    
+    figure(2)
+    [p1,p2] = bisection_coordinate(W,coords,picture);
+    [cut_coord] = cutsize(W,p1);
+    disp('Space to continue ...');
+    pause;
+    
+    figure(3)
+    disp(' ');
+    disp(' 2. A multilevel method from the "Metis 5.0.2" package.');
+    disp(' This will only work if you have Metis and its Matlab interface.');
+    disp('  ');
+    [p1,p2] = bisection_metis(W,coords,picture);
+    [cut_metis] = cutsize(W,p1);
+    disp('  ');
+    disp(' Hit space to continue ...');
+    pause;
+        
+            
+    disp(' ');
+    disp(' 3. Spectral partitioning, which uses the second eigenvector of');
+    disp(' the Laplacian matrix of the graph, also known as the "Fiedler vector".');
+    disp('  ');
+    figure(6)
+    [p1,p2] = bisection_spectral(W,coords,picture);    
+    [cut_spectral] = cutsize(W,p1);
+    disp('  ');    
+    disp(' Hit space to continue ...');
+    pause;
+    
+    figure(7)
+    disp(' ');
+    disp(' 4. Inertial partitioning, which uses the coordinates to find');
+    disp(' a separating line in the plane.');
+    disp('  ');
+    [p1,p2] = bisection_inertial(W,coords,picture);    
+    [cut_inertial] = cutsize(W,p1);
+    disp('  ');
+    disp(' Hit space to continue ...');
+    pause;
+
+    close all;   
+    format;
+    
+    disp(' ');
+    disp('          *****************************************************************')
+    disp('          ***        Bisection Benchmark                               *** ');
+    disp('          ***        D.P. & O.S. for Numerical Computing, USI Lugano   *** ');
+    disp('          *****************************************************************')
+    disp(' ');
+    disp(' ');
+    
+    
+end

mp3/Project_3_Maggioni_Claudio/src/Bench_metis.m

@@ -0,0 +1,20 @@
+function [cut_recursive,cut_kway] = Bench_metis(picture)
+% Compare recursive bisection and direct k-way partitioning,
+% as implemented in the Metis 5.0.2 library.
+
+%  Add necessary paths
+addpaths_GP;
+
+% Graphs in question
+% load 'airfoil1.mat' ;
+% load 'crack.mat';
+
+% Steps
+% 1. Initialize the cases
+% 2. Call metismex to
+%     a) Recursively partition the graphs in 16 and 32 subsets.
+%     b) Perform direct k-way partitioning of the graphs in 16 and 32 subsets.
+% 3. Visualize the results for 32 partitions
+
+
+end

mp3/Project_3_Maggioni_Claudio/src/Bench_rec_bisection.m

@@ -0,0 +1,80 @@
+% Benchmark for recursively partitioning meshes, based on various
+% bisection approaches
+%
+% D.P & O.S for Numerical Computing in USI
+
+
+
+% add necessary paths
+addpaths_GP;
+nlevels_a = 3;
+nlevels_b = 4;
+
+fprintf('       *********************************************\n')
+fprintf('       ***  Recursive graph bisection benchmark  ***\n');
+fprintf('       *********************************************\n')
+
+% load cases
+cases = {
+    'airfoil1.mat';
+    '3elt.mat';
+    'barth4.mat';
+    'mesh3e1.mat';
+    'crack.mat';
+    };
+
+nc = length(cases);
+maxlen = 0;
+for c = 1:nc
+    if length(cases{c}) > maxlen
+        maxlen = length(cases{c});
+    end
+end
+
+for c = 1:nc
+    fprintf('.');
+    sparse_matrices(c) = load(cases{c});
+end
+
+
+fprintf('\n\n Report Cases         Nodes     Edges\n');
+fprintf(repmat('-', 1, 40));
+fprintf('\n');
+for c = 1:nc
+    spacers  = repmat('.', 1, maxlen+3-length(cases{c}));
+    [params] = Initialize_case(sparse_matrices(c));
+    fprintf('%s %s %10d %10d\n', cases{c}, spacers,params.numberOfVertices,params.numberOfEdges);
+end
+
+%% Create results table
+fprintf('\n%7s %16s %20s %16s %16s\n','Bisection','Spectral','Metis 5.0.2','Coordinate','Inertial');
+fprintf('%10s %10d %6d %10d %6d %10d %6d %10d %6d\n','Partitions',8,16,8,16,8,16,8,16);
+fprintf(repmat('-', 1, 100));
+fprintf('\n');
+
+
+for c = 1:nc
+    spacers = repmat('.', 1, maxlen+3-length(cases{c}));
+    fprintf('%s %s', cases{c}, spacers);
+    sparse_matrix = load(cases{c});
+    
+
+    % Recursively bisect the loaded graphs in 8 and 16 subgraphs.
+    % Steps
+    % 1. Initialize the problem
+    [params] = Initialize_case(sparse_matrices(c));
+    W      = params.Adj;
+    coords = params.coords;    
+    % 2. Recursive routines
+    % i. Spectral    
+    % ii. Metis
+    % iii. Coordinate    
+    % iv. Inertial
+    % 3. Calculate number of cut edges
+    % 4. Visualize the partitioning result
+    
+    
+    fprintf('%6d %6d %10d %6d %10d %6d %10d %6d\n',0,0,...
+    0,0,0,0,0,0);
+    
+end

mp3/Project_3_Maggioni_Claudio/src/Contents.m

@@ -0,0 +1,45 @@
+% Graph Partitioning Toolbox. 
+% 
+% Datasets
+%   2d_meshes        -  Graphs from aeronautics applications
+%   Mesh_generation  -  Simple 2d and 3d mesh generation routines
+% 
+% External
+%   metismex - Interface between Metis 5.0.2 and Matlab
+% 
+% Bisection methods.
+%   bisection_spectral   - Spectral bisection
+%   bisection_coordinate - Coordinate bisection.
+%   bisection_inertial   - Inertial bisection.
+%   bisection_metis      - Multilevel method from Metis.
+%   partition            - Partition points by a plane.
+% 
+% Multiway partitions.
+%   rec_bisection  - Use any bisection method to get a multiway partition.
+% 
+% 
+% Benchmark scripts
+%   Bench_bisection        - Compare various graph bisection algorithms
+%   Bench_eigen_plot       - Visualize information from the eigenspectrum 
+%                            of the graph Laplacian
+%   Bench_metis_comparison - Compare recursive bisection and direct k-way partitioning,
+%                            as implemented in the Metis 5.0.2 library.
+%   Bench_rec_bisection    - Benchmark for recursively partitioning meshes, based on various
+%                            bisection approaches
+% 
+% 
+% Visualization and graphics.
+%   gplotpart      - Draw a 2-way partition.
+%   gplotmap       - Draw a multiway partition.
+%   highlight      - Draw a mesh with some vertices highlighted.
+%   gplotg         - Draw a 2D or 3D mesh (replaces Matlab's gplot).
+%   etreeplotg     - Draw an elimination tree (replaces Matlab's etreeplot).
+%   spypart        - Matrix spy plot with partition boundaries.
+%  
+% Miscellaneous.
+%   cutsize                 - Find or count edges cut by a partition.
+%   other                   - Other side of a partition.
+%   adjacency_to_incidence  - create the incidence matrix of the graph
+%   Contents                - contents of the partitioning toolbox 
+%   blockdiags              - Create matrix with specified block diagonals.
+

mp3/Project_3_Maggioni_Claudio/src/Initialize_case.m

@@ -0,0 +1,18 @@
+function [params] =  Initialize_case(SM_case)
+
+% Adjacency
+params.Adj    = SM_case.Problem.A;
+% Incidence
+params.Inc    = adjacency_to_incidence(params.Adj);
+% Edges
+params.numberOfEdges = size(params.Inc,1);
+% Vertices
+params.numberOfVertices = size(params.Inc,2);
+% coords (if exist)
+if isfield(SM_case, 'Problem') && isfield(SM_case.Problem, 'aux') && ...
+   isfield(SM_case.Problem.aux, 'coord')
+    params.coords = SM_case.Problem.aux.coord;
+end
+    
+    
+end

mp3/Project_3_Maggioni_Claudio/src/Visualization/gplotg.m

@@ -0,0 +1,49 @@
+function handle = gplotg(A,xy,lc)
+% GPLOTG : Plot a "graph theoretic" graph.
+%
+%	handle = gplotg(A,xy) plots the graph specified by A and xy.
+%	A graph, G, is a set of nodes numbered from 1 to n,
+%	and a set of connections, or edges, between them.
+%	In order to plot G, two matrices are needed.
+%	The adjacency matrix, A, has a(i,j) nonzero if and
+%	only if node i is connected to node j.  The coordinates
+%	array, xy, is an n-by-2 or n-by-3  matrix with the position 
+%       for node i in the i-th row, xy(i,:) = [x(i) y(i)],
+%                               or  xy(i,:) = [x(i) y(i) z(i)].
+%	
+%	gplotg(A,xy,lc) uses line type and color instead of the 
+%	default, 'r-'.   For example, lc = 'g:'.  See PLOT.
+
+
+
+if nargin < 3
+    lc = 'r-';
+end
+
+[i, j] = find(A);
+[~, p] = sort(max(i,j));
+i = i(p);
+j = j(p);
+
+% Create a long, NaN-seperated list of line segments,
+% rather than individual segments.
+
+X = [ xy(i,1) xy(j,1) NaN*ones(size(i))]';
+Y = [ xy(i,2) xy(j,2) NaN*ones(size(i))]';
+X = X(:);
+Y = Y(:);
+
+if size(xy,2) == 2
+    h = plot (X, Y, lc);
+else
+    set(gca,'drawmode','fast');
+    Z = [ xy(i,3) xy(j,3) NaN*ones(size(i))]';
+    Z = Z(:);
+    h = plot3 (X, Y, Z, lc);
+end
+
+axis equal;
+axis off;
+if nargout >= 1
+    handle = h;
+end

mp3/Project_3_Maggioni_Claudio/src/Visualization/gplotmap.m

@@ -0,0 +1,90 @@
+function handle = gplotmap(A,xy,map,pcolor,ecolor)
+% GPLOTMAP : Plot a partitioned graph in 2 or 3 dimensions.
+%
+%	gplotmap(A,xy,map) plots the n-vertex graph specified 
+%       by the n by n adjacency (or Laplacian) matrix A 
+%       and the n by 2 or 3 matrix of coordinates xy.
+%       Argument map is an n-vector of part numbers, which
+%       assigns each vertex to a part.
+%
+%       By default, edges that join different parts are omitted, and
+%       the picture shows each part in a different color.  The call
+%
+%       gplotmap(A,xy,map,pcolor,ecolor)
+%
+%       uses color "pcolor" for the parts and "ecolor" for the
+%       edges between the parts.  If "pcolor" has multiple rows,
+%       each part gets the color of one row (in cyclic order).
+
+
+[n,n] = size(A);
+if nargin < 3, 
+    map = zeros(1,n);
+end;
+if nargin < 4 
+    pcolor = [];
+end;
+if nargin < 5,
+    ecolor = [];
+end;
+
+parts = setfilter(map);
+nparts = length(parts);
+if length(pcolor) == 0,
+    pcolor = hsv(nparts);
+    pcolor = pcolor(randperm(nparts),:);
+end;
+
+clf reset
+colordef(gcf,'black')
+if size(xy,2) == 2
+    axis([min(xy(:,1)) max(xy(:,1)) min(xy(:,2)) max(xy(:,2))]);
+else
+    axis([min(xy(:,1)) max(xy(:,1)) min(xy(:,2)) max(xy(:,2)) ...
+        min(xy(:,3)) max(xy(:,3))]); 
+end;
+
+axis equal;
+axis off;
+
+hold on
+
+% Count and plot the separating edges.
+[i,j] = find(A);
+f = find(map(i) > map(j));
+
+%comm volume
+volume = 0;
+for r=1:size(A,1)
+    idx = find(i == r);
+    volume = volume + length(unique(map(j(idx)))) - 1;
+end
+
+if length(f) 
+    xlabel( [int2str(length(f)) ' cut edges on ' int2str(nparts) ' partitions, '...
+        'communication volume ' int2str(volume) ],'visible','on');
+    if length(ecolor)
+        cut = sparse(i(f),j(f),1,n,n);
+        set(gplotg(cut,xy,'-'),'color',ecolor);
+    end;
+else
+    xlabel('0 cut edges','visible','on');
+end;
+
+% Plot each piece.
+ncolor = 1;
+for partnumber = parts;
+    c = pcolor(ncolor,:);
+    ncolor = rem(ncolor, size(pcolor,1)) + 1;
+    part = find(map == partnumber);
+    set(gplotg(A(part,part),xy(part,:),'-'),'color',c);
+    if n < 500,
+        if size(xy,2) == 2
+            set(plot(xy(part,1),xy(part,2),'o'),'color',c);
+        else
+            set(plot3(xy(part,1),xy(part,2),xy(part,3),'o'),'color',c);
+        end;
+    end;
+end;
+
+hold off

mp3/Project_3_Maggioni_Claudio/src/Visualization/gplotpart.m

@@ -0,0 +1,76 @@
+function handle = gplotpart(A,xy,part1,color1,color2,color3)
+% GPLOTPART : Plot a partitioned graph in 2 or 3 dimensions.
+%
+%	gplotpart(A,xy,part1) plots the n-vertex graph specified 
+%       by the n by n adjacency (or Laplacian) matrix A 
+%       and the n by 2 or 3 matrix of coordinates xy.
+%       Argument part1 is a vector of vertex names (integers 1:n);
+%       the subgraphs induced by part1 and its complement are plotted
+%       in different colors, with the edges joining them in a third color.
+%       Three more optional arguments give the three colors.
+
+
+if nargin < 3
+    part1 = [];
+end
+if nargin < 4
+    color1 = 'yellow';
+end
+if nargin < 5
+    color2 = 'cyan';
+end
+if nargin < 6
+    color3 = 'red';
+end
+
+[n,~] = size(A);
+part1 = part1(:);
+part2 = 1:n;
+part2(part1)=zeros(size(part1));
+part2 = find(part2);
+part2 = part2(:);
+cut = spaugment(A(part1,part2),1);
+cutxy = xy([part1; part2],:);
+
+clf reset
+%colordef(gcf,'black')
+if size(xy,2) == 2
+    axis([min(xy(:,1)) max(xy(:,1)) min(xy(:,2)) max(xy(:,2))]);
+else
+    axis([min(xy(:,1)) max(xy(:,1)) min(xy(:,2)) max(xy(:,2)) ...
+        min(xy(:,3)) max(xy(:,3))]); 
+end
+
+axis equal;
+axis off;
+
+hold on
+if ~isempty(part1) && ~isempty(part2)
+    set(gplotg(cut,cutxy,'-'),'color',color3);
+    xlabel([int2str(cutsize(A,part1)) ' cut edges'],'visible','on');
+else
+    xlabel('0 cut edges','visible','on');
+end
+
+if ~isempty(part1) 
+    set(gplotg(A(part1,part1),xy(part1,:),'-'),'color',color1);
+    if n < 500
+        if size(xy,2) == 2
+            set(plot(xy(part1,1),xy(part1,2),'o'),'color',color1);
+        else
+            set(plot3(xy(part1,1),xy(part1,2),xy(part1,3),'o'),'color',color1);
+        end
+    end
+end
+
+if ~isempty(part2)
+    set(gplotg(A(part2,part2),xy(part2,:),'-'),'color',color2);
+    if n < 500
+        if size(xy,2) == 2
+            set(plot(xy(part2,1),xy(part2,2),'o'),'color',color2);
+        else
+            set(plot3(xy(part2,1),xy(part2,2),xy(part2,3),'o'),'color',color2);
+        end
+    end
+end
+hold off

mp3/Project_3_Maggioni_Claudio/src/Visualization/highlight.m

@@ -0,0 +1,60 @@
+function highlight(A,xy,sep,highcolor,meshcolor,dotsize)
+% HIGHLIGHT : Plot a mesh with subgraph highlighted.
+%
+%  highlight(A,xy,sep) plots a picture of the mesh A with coordinates xy,
+%  highlighting the subgraph induced by the vertices in sep.
+%  Optional fourth and fifth arguments are color/linetype of highlight and mesh.
+%  Optional sixth argument is size of highlighting dot.
+
+
+
+if nargin < 4
+    highcolor = 'w-';
+end
+if nargin < 5
+    meshcolor = 'r-';
+end
+if nargin < 6
+    dotsize = 15;
+end
+
+[n,n] = size(A);
+[i,j] = find(A);
+
+% Plot main graph with vertices before edges.
+
+[ignore, p] = sort(max(i,j));
+i = i(p);
+j = j(p);
+
+% Create a long, NaN-seperated list of line segments,
+% rather than individual segments.
+
+X = [ xy(i,1) xy(j,1) NaN*ones(size(i))]';
+Y = [ xy(i,2) xy(j,2) NaN*ones(size(i))]';
+X = X(:);
+Y = Y(:);
+
+xymin = min(xy);
+xymax = max(xy);
+plot (X, Y, meshcolor,'erasemode','none');
+axis([xymin(1) xymax(1) xymin(2) xymax(2)]);
+axis('equal');
+axis('off');
+hold on;
+
+% Highlight sep set.
+
+B = A(sep,sep);
+xB = xy(sep,1);
+yB = xy(sep,2);
+[i,j] = find(B);
+X = [xB(i) xB(j)]';
+Y = [yB(i) yB(j)]';
+plot (X, Y, highcolor,'erasemode','none');
+if n < 1200
+   handle = plot(xB,yB,[highcolor(1) '.'],'erasemode','none');
+   set(handle,'markersize',dotsize)
+end
+
+hold off;

mp3/Project_3_Maggioni_Claudio/src/Visualization/spypart.m

@@ -0,0 +1,39 @@
+function spypart(S, rp, cp);
+%SPYPART Spy plot with partitioning.
+%   SPYPART(S,rp,cp) plots the sparsity pattern of a matrix S,
+%   with lines marking a block partition described by 
+%   rp (rows) and cp (columns).
+%   If S is square, cp may be omitted and defaults to rp.
+%
+%   Partitions are specified as in the output of DMPERM:
+%   There are length(rp)-1 row blocks, of sizes diff(rp), with rp(1)=1.
+
+
+if (nargin < 3), cp = rp; end
+
+clf
+colordef(gcf,'black')
+[m,n] = size(S);
+if max(m,n) > 100 
+    spy(S,3)
+else
+    spy(S,5)
+end
+hold on
+
+if length(rp) > 2
+   k = length(rp)-2;
+   X = [zeros(1,k); n+ones(1,k)];
+   Y = rp(2:k+1) - 0.5;
+   Y = [Y; Y];
+   plot(X,Y,'w-')
+end
+if length(cp) > 2
+   k = length(cp)-2;
+   X = cp(2:k+1) - .5;
+   X = [X; X];
+   Y = [zeros(1,k); m+ones(1,k)];
+   plot(X,Y,'w-') 
+end
+axis('ij')
+hold off

mp3/Project_3_Maggioni_Claudio/src/addpaths_GP.m

@@ -0,0 +1,8 @@
+% add necessary paths for the Partitioning Toolbox
+
+addpath ../datasets/
+addpath ../datasets/Mesh_generation/
+addpath ../datasets/2d_meshes
+addpath ../datasets/Roads/
+addpath ../external/
+addpath Visualization/

mp3/Project_3_Maggioni_Claudio/src/adjacency_to_incidence.m

@@ -0,0 +1,54 @@
+function [Inc] = adjacency_to_incidence(Adj)
+
+A = Adj;
+for i = 1:size(A,1)    
+    A(i,i) = 0;    
+end
+
+Inc = adj_2_inc(A);
+% Inc = adj2inc(A);
+Inc = Inc';
+end
+
+function Ic = adj_2_inc(A)
+
+% adjacency2incidence - convert an adjacency matrix to an incidence matrix
+%
+%   Ic = adjacency2incidence(A);
+%
+%   A(i,j) = 1 iff (i,j) is an edge of the graph.
+%   For each edge number k of the graph linking (i,j)
+%       Ic(i,k)=1 and Ic(j,k)=-1 
+%
+%   Ic is a sparse matrix.
+%   Ic is also known as the graph gradient.
+%
+%   Copyright (c) 2006 Gabriel Peyre
+
+%% compute list of edges
+[i,j,s] = find(sparse(A));
+I = find(i<=j);
+i = i(I);
+j = j(I);
+% number of edges
+n = length(i);
+% number of vertices
+nverts = size(A,1);
+
+%% build sparse matrix
+s = [ones(n,1); -ones(n,1)];
+is = [(1:n)'; (1:n)'];
+js = [i(:); j(:)];
+Ic = sparse(is,js,s,n,nverts);
+Ic = Ic';
+
+% fix self-linking problem (0)
+a = find(i==j);
+if not(isempty(a))
+    for t=a'
+        Ic(i(t),t) = 1;
+    end
+end
+
+end
+

mp3/Project_3_Maggioni_Claudio/src/bisection_coordinate.m

@@ -0,0 +1,33 @@
+function [part1,part2] = bisection_coordinate(A,xy,picture)
+% COORDPART : Coordinate bisection partition of a mesh.
+%
+% [part1,part2] = bisection_coordinate(A,xy,picture) returns a list of the vertices on one side of a partition
+%     obtained by bisection perpendicular to a coordinate axis.  We try every
+%     coordinate axis and return the best cut.
+%     Input A is the adjacency matrix of the mesh; 
+%     each row of xy is the coordinates of a point in d-space.
+%
+% coordpart(A,xy,1) also draws a picture.
+
+d = size(xy,2);
+best_cut = inf;
+for dim = 1:d
+    v = zeros(d,1);
+    v(dim) = 1;
+    [p1,p2] = partition(xy,v);
+    this_cut   = cutsize(A,p1);
+    if this_cut < best_cut
+        best_cut = this_cut;
+        part1 = p1;
+        part2 = p2;
+    end
+end
+
+if picture
+    clf reset
+    gplotpart(A,xy,part1);
+    title('Coordinate bisection')
+end
+
+
+end

mp3/Project_3_Maggioni_Claudio/src/bisection_inertial.m

@@ -0,0 +1,55 @@
+function [part1,part2] = bisection_inertial(A,xy,picture)
+% INERTPART : Inertial partition of a graph.
+%
+% p = inertpart(A,xy) returns a list of the vertices on one side of a partition
+%     obtained by bisection with a line or plane normal to a moment of inertia
+%     of the vertices, considered as points in Euclidean space.
+%     Input A is the adjacency matrix of the mesh (used only for the picture!);
+%     each row of xy is the coordinates of a point in d-space.
+%
+% inertpart(A,xy,1) also draws a picture.
+%
+% See also PARTITION
+
+
+%disp(' ');
+%disp(' Numerical Computing @ USI Lugano:   ');
+%disp(' Implement inertial bisection');
+%disp(' ');
+
+
+% Steps
+% 1. Calculate the center of mass.
+m = sum(xy, 1) / size(xy, 1);
+xm = m(1);
+ym = m(2);
+
+% 2. Construct the matrix M.
+%  (Consult the pdf of the assignment for the creation of M) 
+h = [(xy(:, 1) - xm) (xy(:, 2) - ym)];
+X1 = sum(h(:, 1) .^ 2);
+X2 = sum(h(:, 1) .* h(:, 2));
+X3 = sum(h(:, 2) .^ 2);
+
+M = [X1 X2; X2 X3];
+
+% 3. Calculate the smallest eigenvector of M.
+[u, ~] = eigs(M,1,'smallestabs');
+a = u(1);
+b = u(2);
+
+% 4. Find the line L on which the center of mass lies.
+S = -b * (xy(:, 1) - xm) + a * (xy(:, 2) - ym);
+Sm = median(S);
+
+% 5. Partition the points around the line L.
+%   (you may use the function partition.m)
+part1 = find(S < Sm);
+part2 = find(S >= Sm);
+
+if picture == 1
+    gplotpart(A,xy,part1);
+    title('Inertial bisection using the Fiedler Eigenvector');
+end
+
+end

mp3/Project_3_Maggioni_Claudio/src/bisection_metis.m

@@ -0,0 +1,25 @@
+function [part1,part2] = bisection_metis(A,xy,picture)
+% METISPART : Partition a graph using Metis default method.
+%
+% p = metispart(A) returns a list of the vertices on one side of
+%     a partition obtained by Metis 4.0 applied to the graph of A.
+%     
+% Optional arguments:
+%   metispart(A,xy) draws a picture of the partitioned graph,
+%                   using the rows of xy as vertex coordinates.
+%   [p1,p2] = metispart(...)   also returns the list of vertices
+%                              on the other side of the partition.
+%
+% See also METISMEX (which accepts all the Metis options), 
+
+
+map = metismex('PartGraphRecursive',A,2);
+[part1,part2] = other(map);
+
+if picture
+    gplotpart(A,xy,part1);
+    title('Metis bisection')
+end
+
+
+end

mp3/Project_3_Maggioni_Claudio/src/bisection_spectral.m

@@ -0,0 +1,64 @@
+function [part1,part2] = bisection_spectral(A,xy,picture)
+% bisection_spectral : Spectral partition of a graph.
+%
+% [part1,part2] = bisection_spectral(A) returns a partition of the n vertices
+%                 of A into two lists part1 and part2 according to the
+%                 spectral bisection algorithm of Simon et al:  
+%                 Label the vertices with the components of the Fiedler vector
+%                 (the second eigenvector of the Laplacian matrix) and partition
+%                 them around the median value or 0.
+
+
+
+%disp(' ');
+%disp(' Numerical Computing @ USI Lugano:   ');
+%disp(' Implement inertial bisection');
+%disp(' ');
+
+% Steps
+
+% 1. Construct the Laplacian.
+[i, j, ~] = find(triu(A));
+G = graph(i, j);
+L = laplacian(G);
+
+% 2. Calculate its eigensdecomposition.
+[V, ~] = eigs(L,2,'smallestabs');
+W = V(:,2);
+
+% 3. Label the vertices with the components of the Fiedler vector.
+[values, indexes] = sort(W);
+
+% 4. Partition them around their median value, or 0.
+cut = size(values(values < 0), 1);
+part1 = indexes(1:cut);
+part2 = indexes(cut+1:end);
+
+if picture == 1
+    gplotpart(A,xy,part1); 
+    title('Spectral bisection using the Fiedler Eigenvector');
+elseif picture == 2
+    drawgraph(A, xy, W, part1, part2)
+    % 3d plot
+    title('Spectral bisection using the Fiedler Eigenvector');
+end
+end
+
+function drawgraph(A, xy, W, part1, part2)
+    scatter3(xy(:, 1), xy(:, 2), W, 30, W, 'filled')
+    colorbar
+    
+    hold on
+    
+    gplot(A(part1, part1), xy(part1, :), '-k');
+    
+    set(gca,'ColorOrder', ones(size(part2, 2), 3) .* 0.3)
+    gplot(A(part2, part2), xy(part2, :));
+    
+    P1P2 = zeros(size(A));
+    P1P2(part1, part2) = A(part1, part2);
+    P1P2(part2, part1) = A(part2, part1);
+    gplot(P1P2, xy, '-r');
+    
+    hold off
+end

mp3/Project_3_Maggioni_Claudio/src/blockdiags.m

@@ -0,0 +1,86 @@
+function A = blockdiags(B,d,m,n)
+% BLOCKDIAGS : Create sparse block diagonal matrices.
+%
+% A = blockdiags(B,d,m,n).  
+%
+%   Blockdiags, which generalizes the function "spdiags", 
+%   produces a sparse matrix with specified block diagonals.
+%
+%   A is an m*k-by-n*k matrix, or an m-by-n matrix of k-by-k blocks.  
+%       The nonzero blocks of A are located on p block diagonals.  
+%   B is a min(m,n)*k-by-p*k matrix whose k-by-k block columns
+%       are the block diagonals of A.  
+%   (Alternatively, B is k-by-p*k, and then A is block Toeplitz.)
+%   d is a vector of p integers in the range -m+1 : n-1, 
+%       specifying which block diagonals in A are to be nonzero.
+%   The values of p and k are determined from the dimensions of B and d.
+%
+%   For k=1 this is exactly the same as A = spdiags(B,d,m,n); see spdiags 
+%   for examples of use.  For k>1 this is conceptually the same as spdiags,
+%   but k-by-k blocks replace matrix elements everywhere.
+%
+%   For example, the following code sets A to the n^2-by-n^2 matrix of 
+%   the Laplacian on an n-by-n square grid; the matrix is block tridiagonal, 
+%   and the nonzero blocks themselves are tridiagonal or the identity.
+%
+%        a = blockdiags ([-1 4 -1], -1:1, n, n);
+%        I = speye (n, n);
+%        A = blockdiags ([-I a -I], -1:1, n, n);
+
+
+if nargin ~= 4
+   error ('Usage: A = blockdiags(B,d,m,n)');
+end
+
+k = length(d);
+[nrB,ncB] = size(B);
+p = ncB/k;
+
+% Check for reasonable input.
+
+if any(size(m)>1) | any(size(n)>1)
+   error ('blockdiags(B,d,m,n): m or n not scalar');
+end
+if min(size(d)~=1) 
+   error ('blockdiags(B,d,m,n): d not a vector');
+end
+if any(rem(size(B),p))
+   error ('blockdiags(B,d,m,n): block size does not divide size of B');
+end
+
+
+% Process A in compact form.
+
+[i,j,a] = find(B);
+nzB = length(a);
+
+% Duplicate rows of B if A is to be block Toeplitz.
+
+if nrB == p
+   nzA = m*nzB;
+   range = [0:nzA-1]';
+   r = 1 + rem(range,nzB);
+   a = a(r);
+   i = i(r);
+   j = j(r);
+   a = a(:);
+   i = i(:) + p * floor(range/nzB);
+   j = j(:);
+end
+
+% Rows of B are rows of A.  Shift columns appropriately.
+
+dd = d(1+floor((j-1)/p));
+dd = dd(:);
+j = 1 + rem(j-1,p) + p * (dd + floor((i-1)/p));
+
+% Clip columns that shifted off the matrix.
+
+f = find(j>=1 & j<=n*p);
+a = a(f);
+i = i(f);
+j = j(f);
+
+% Put A back into sparse form.
+
+A = sparse (i,j,a,m*p,n*p);

mp3/Project_3_Maggioni_Claudio/src/cutsize.m

@@ -0,0 +1,26 @@
+function [ne,edges] = cutsize(A,map);
+% CUTSIZE : Edges cut by a vertex partition.
+%
+% ne = cutsize(A,part)        : part is a list of vertices of A;
+%                               ne is the number of edges crossing the cut.
+% [ne,edges] = cutsize(A,part): Same, edges is the crossing edges as rows [i j].
+%
+% "part" may also be a map, with one entry per vertex.
+
+
+% Take input in either order unless A is 1 by 1.
+if size(map,1) == size(map,2) & size(map,1) > 1
+    t = A;
+    A = map;
+    map = t;
+end;
+
+% Convert input to a map in any case.
+if length(map) ~= length(A) | max(map) == length(A)
+    [ignore,map] = other(map,A);
+end;
+
+[i,j] = find(A);
+f = find(map(i) > map(j));
+ne = length(f);
+edges = [i(f) j(f)];

mp3/Project_3_Maggioni_Claudio/src/ex2_bisection_table.m

@@ -0,0 +1,32 @@
+load Toy_meshes;
+addpaths_GP;
+
+[W1,c1] = grid5rec(10,100);
+[W2,c2] = grid5rec(100,10);
+[W3,c3] = grid5recRotate(100,10,-45);
+[W4,c4] = gridt(40);
+[W5,c5] = grid9(30);
+
+W = {W1, W2, W3, W4, W5, Smallmesh, Tapir, Eppstein};
+c = {c1, c2, c3, c4, c5, Smallmesh_coords, Tapir_coords, Eppstein_coords};
+
+fprintf('implementation\tgrid5rec(10,100)\tgrid5rec(100,10)\tgrid5recRotate(100,10,-45)\tgridt(40)\tgrid9(30)\tSmallmesh\tTapir\tEppstein\n');
+fprintf('coordinate\t');
+runtest(@(a, b) bisection_coordinate(a, b, 0));
+fprintf('metis     \t');
+runtest(@(a, b) bisection_metis(a, b, 0));
+fprintf('spectral\t');
+runtest(@(a, b) bisection_spectral(a, b, 0));
+fprintf('inertial\t');
+runtest(@(a, b) bisection_inertial(a, b, 0));
+
+function runtest(implementation)
+    global W
+    global c
+    for i=1:8
+        out = implementation(W{i}, c{i});
+        ne = cutsize(W{i},out);
+        fprintf('%d\t', ne);
+    end
+    fprintf('\n');
+end

mp3/Project_3_Maggioni_Claudio/src/ex3_plots.m

@@ -0,0 +1,14 @@
+load Toy_meshes;
+addpaths_GP;
+
+[A, xy] = grid9(50);
+bisection_spectral(A, xy, 2);
+%matlab2tikz('../../ex3_grid9.tex');
+
+figure;
+bisection_spectral(Smallmesh, Smallmesh_coords, 2);
+%matlab2tikz('../../ex3_small.tex');
+
+figure;
+bisection_spectral(Eppstein, Eppstein_coords, 2);
+%matlab2tikz('../../ex3_eppstein.tex');

mp3/Project_3_Maggioni_Claudio/src/other.m

@@ -0,0 +1,26 @@
+function [out1,out2] = other(in1,in2);
+% OTHER : Find the other part of a partition, or
+%         convert a partition to the other representation.
+%
+% part2 = other(part1,n)     : part1 is a list of subscripts in 1:n;
+%                              part2 becomes the complementary list.
+% part2 = other(part1,A)     : Same, except n is taken as max(size(A)).
+% [part2,p] = other(part1,n) : Also returns 0/1 partition vector of length n.
+% [part2,p] = other(part1,A) : Same.
+% [part1,part2] = other(p)   : Converts 0/1 vector to lists of subscripts.
+
+
+if nargin == 1
+
+    out1 = find(in1 == 0);
+    out2 = find(in1 ~= 0);
+
+else
+
+    if max(size(in2)) > 1
+        in2 = max(size(in2));
+    end;
+    out2 = full(spones(sparse(1,in1,1,1,in2)));
+    out1 = find(out2 == 0);
+
+end;

mp3/Project_3_Maggioni_Claudio/src/partition.m

@@ -0,0 +1,39 @@
+function [a,b] = partition(xyz,v)
+% PARTITION : Partition points by a plane.
+%
+% [a,b] = partition(xyz,v):
+% Each row of xyz is an input point in d-space.
+% Input v is a vector, giving a direction normal to the partitioning plane.
+%
+% The output is two vectors of integers, 
+% the indices of the points on each side of the plane.  
+% Points on the plane are assigned to balance the cut.
+
+
+[n,d] = size(xyz);
+
+v = v(:); % Make v a column vector
+
+if length(v) ~= d
+    error('v must be a d-vector')
+end;
+
+dotprod = xyz * v;
+split = median(dotprod);
+a = find(dotprod <  split);
+b = find(dotprod > split);
+c = find(dotprod == split);
+nc = length(c);
+if nc
+    na = length(a);
+    nca = max([ceil(n/2)-na, 0]);
+    nca = min(nca,nc);
+    if nca > 0
+        a = [a; c(1:nca)];
+    end;
+    if nca < nc
+        b = [b; c(nca+1:nc)];
+    end;
+end;
+a = a';
+b = b';

mp3/Project_3_Maggioni_Claudio/src/rec_bisection.m

@@ -0,0 +1,61 @@
+function [map,sepij,sepA] = rec_bisection(method,levels,A,varargin)
+% DICE Separate a graph recursively.
+%
+%   [map,sepij,sepA] = DICE(method,levels,A,arg...) partitions the
+%   mesh or graph A recursively by a specified method. 'method' is the name
+%   of the 2-way edge separator function to call. levels is the number of
+%   levels of partitioning. A is the adjacency matrix of the graph. arg2,
+%   arg3, arg4 are optional additional arguments to the 2-way function. arg2
+%   is special: If it has the same number of rows as A, then the recursive
+%   calls use the appropriate subset of arg2 as well as of A. This is useful
+%   if the partitioner uses xy coords. map is a vector of integers from 0 to
+%   2^levels-1, indexed by vertex, giving the partition number for each
+%   vertex. sepij is a list of separating edges; each row is an edge [i j].
+%   sepA is the adjacency matrix of the graph less the separating edges.
+%
+%       For example,  dice('geopart',7,A,xy,0,ntries) splits A into 2^7=128
+%       parts, using a call at each stage that looks like
+%       geopart(A,xy,0,ntries).
+
+minpoints = 8;   % Don't separate pieces smaller than this.
+
+nargcall = nargin - 2;
+
+n = size(A,1);
+
+if n < minpoints || levels < 1
+    map = zeros(1,n);
+else
+    % Call the partitioner to split the graph, giving one part.
+    [p1,p2] = feval(method, A, varargin{:});
+    
+    % Call recursively.
+    vararginp1 = varargin;
+    vararginp2 = varargin;
+    if nargcall >= 2
+        if size(varargin{1},1) == n
+            vararginp1{1} = varargin{1}(p1,:);
+            vararginp2{1} = varargin{1}(p2,:);
+        end
+    end
+    
+    mapa = rec_bisection(method,levels-1,A(p1,p1),vararginp1{:});
+    mapb = rec_bisection(method,levels-1,A(p2,p2),vararginp2{:});
+    
+    % Set up the whole map.
+    map = zeros(1,n);
+    mapb = mapb + max(mapa) + 1;
+    map(p1) = mapa;
+    map(p2) = mapb;
+end
+
+% Set up the separating edge list and separated graph.
+if nargout >= 2
+    [i,j] = find(A);
+    f = find(map(i) > map(j));
+    sepij = [i(f) j(f)];
+    f = find(map(i) == map(j));
+    sepA = sparse(i(f), j(f), 1, n, n);
+end
+
+end

mp3/Project_3_Maggioni_Claudio/src/setfilter.m

@@ -0,0 +1,15 @@
+function B = setfilter(A)
+% SETFILTER : Sort and remove duplicates.
+%
+% B = setfilter(A) returns a row vector with one occurrence
+%     of each different element of A, in sorted order.
+
+
+if length(A) == 0
+    B = []; 
+    return;
+end;
+
+B = sort(A(:));
+B(find(diff(B)==0)) = [];
+B = B.';

mp3/assignment.sty

@@ -0,0 +1,95 @@
+\usepackage{ifthen}
+\usepackage[utf8]{inputenc}
+\usepackage{graphics}
+\usepackage{graphicx}
+\usepackage{hyperref}
+
+\pagestyle{plain}
+\voffset -5mm
+\oddsidemargin  0mm
+\evensidemargin -11mm
+\marginparwidth 2cm
+\marginparsep 0pt
+\topmargin 0mm
+\headheight 0pt
+\headsep 0pt
+\topskip 0pt        
+\textheight 255mm
+\textwidth 165mm
+
+\newcommand{\duedate} {}
+\newcommand{\setduedate}[1]{%
+\renewcommand\duedate {Due date:~ #1}}
+\newcommand\isassignment {false}
+\newcommand{\setassignment}{\renewcommand\isassignment {true}}
+\newcommand{\ifassignment}[1]{\ifthenelse{\boolean{\isassignment}}{#1}{}}
+\newcommand{\ifnotassignment}[1]{\ifthenelse{\boolean{\isassignment}}{}{#1}}
+
+\newcommand{\assignmentpolicy}{
+\begin{table}[h]
+\begin{center}
+\scalebox{0.8} {%
+\begin{tabular}{|p{0.02cm}p{16cm}|}
+\hline
+&\\
+\multicolumn{2}{|c|}{\Large\textbf{Numerical Computing 2020 ---  Submission Instructions}}\\
+\multicolumn{2}{|c|}{\large\textbf{(Please, notice that following instructions are mandatory: }}\\
+\multicolumn{2}{|c|}{\large\textbf{submissions that don't comply with, won't be considered)}}\\
+&\\
+\textbullet & Assignments must be submitted to \href{https://www.icorsi.ch/course/view.php?id=10018}{iCorsi} (i.e. in electronic format).\\
+\textbullet & Provide both executable package and sources (e.g. C/C++ files, Matlab). 
+If you are using libraries, please add them in the file. Sources must be organized in directories called:\\
+\multicolumn{2}{|c|}{\textit{Project\_number\_lastname\_firstname}}\\
+& and  the  file must be called:\\
+\multicolumn{2}{|c|}{\textit{project\_number\_lastname\_firstname.zip}}\\
+\multicolumn{2}{|c|}{\textit{project\_number\_lastname\_firstname.pdf}}\\
+\textbullet &  The TAs will grade your project by reviewing your project write-up, and looking at the implementation 
+                 you attempted, and benchmarking your code's performance.\\
+
+\textbullet & You are allowed to discuss all questions with anyone you like; however: (i) your submission must list anyone you discussed problems with and (ii) you must write up your submission independently.\\
+\hline
+\end{tabular}
+}
+\end{center}
+\end{table}
+}
+\newcommand{\punkte}[1]{\hspace{1ex}\emph{\mdseries\hfill(#1~\ifcase#1{Points}\or{Points}\else{Points}\fi)}}
+
+
+\newcommand\serieheader[6]{
+\thispagestyle{empty}%
+\begin{flushleft}
+\includegraphics[width=0.4\textwidth]{usi_inf}
+\end{flushleft}
+  \noindent%
+  {\large\ignorespaces{\textbf{#1}}\hspace{\fill}\ignorespaces{ \textbf{#2}}}\\ \\%
+  {\large\ignorespaces #3 \hspace{\fill}\ignorespaces #4}\\
+  \noindent%
+  \bigskip
+  \hrule\par\bigskip\noindent%
+  \bigskip {\ignorespaces {\Large{\textbf{#5}}}
+  \hspace{\fill}\ignorespaces \large \ifthenelse{\boolean{\isassignment}}{\duedate}{#6}}
+  \hrule\par\bigskip\noindent%  \linebreak
+ }
+
+\makeatletter
+\def\enumerateMod{\ifnum \@enumdepth >3 \@toodeep\else
+      \advance\@enumdepth \@ne
+      \edef\@enumctr{enum\romannumeral\the\@enumdepth}\list
+      {\csname label\@enumctr\endcsname}{\usecounter
+        {\@enumctr}%%%? the following differs from "enumerate"
+	\topsep0pt%
+	\partopsep0pt%
+	\itemsep0pt%
+	\def\makelabel##1{\hss\llap{##1}}}\fi}
+\let\endenumerateMod =\endlist
+\makeatother
+
+
+
+
+\usepackage{textcomp}
+
+
+
+

mp3/project_3_Maggioni_Claudio.tex

@@ -0,0 +1,114 @@
+\documentclass[unicode,11pt,a4paper,oneside,numbers=endperiod,openany]{scrartcl}
+\usepackage{graphicx}
+\usepackage{subcaption}
+\usepackage{amsmath}
+\input{assignment.sty}
+\usepackage{pgfplots}
+\pgfplotsset{compat=newest}
+\usetikzlibrary{plotmarks}
+\usetikzlibrary{arrows.meta}
+\usepgfplotslibrary{patchplots}
+\usepackage{grffile}
+\usepackage{amsmath}
+\usepackage{subcaption}
+
+\hyphenation{PageRank}
+\hyphenation{PageRanks}
+
+\begin{document}
+
+
+\setassignment
+\setduedate{Wednesday, 4 November 2020, 11:55 PM}
+
+\serieheader{Numerical Computing}{2020}{Student: Claudio Maggioni}{Discussed with: --}{Solution for Project 3}{}
+\newline
+
+\assignmentpolicy
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{Install METIS 5.0.2, and the corresponding Matlab mex interface}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\section{Implement various graph partitioning algorithms \punkte{60}}
+I summarize the various benchmark results in table \ref{table:bisection}. Please note that this table
+can be generated at will with the script \texttt{ex2\_bisection\_table.m}.
+
+\section{Visualize the Fiedler eigenvector\punkte{10}}
+
+In figure \ref{fig:run1} there are graph outputs respectively from \textit{Grid9}, \textit{Small}, and \textit{Eppstein}.
+
+\begin{figure}[h]
+\begin{subfigure}{0.5\textwidth}
+\centering
+\includegraphics[trim=50 200 50 200,clip,width=\textwidth]{ex3_grid9.pdf}
+	\caption{Plot for \textit{Grid9}}
+\end{subfigure}
+\begin{subfigure}{0.5\textwidth}
+\centering
+\includegraphics[trim=50 200 50 200,clip,width=\textwidth]{ex3_small.pdf}
+	\caption{Plot for \textit{Small}}
+\end{subfigure}
+\begin{subfigure}{0.5\textwidth}
+\includegraphics[trim=50 200 50 200,clip,width=\textwidth]{ex3_eppstein.pdf}
+	\caption{Plot for \textit{Eppstein}}
+\end{subfigure}
+\caption{Graph outputs for the 3 adjacency matrices.}
+\label{fig:run1}
+\end{figure}
+
+\section{Recursively bisecting meshes \punkte{20}}
+
+Summarize your results in table \ref{table:Rec_bisection}.
+
+\section{Compare recursive bisection to direct $k$-way partitioning\punkte{10}}
+
+Summarize your results in table \ref{table:Compare_Metis}.
+
+\begin{table}[h]
+	\caption{Bisection results}
+	\centering
+	\begin{tabular}{|l|r|r|r|r|} \hline\hline
+		Mesh             &  Coordinate & Metis 5.0.2  & Spectral & Inertial  \\ \hline
+		grid5rect(10,100)&   10 &        10 &           10  &      10 \\
+		grid5rect(100,10)&   10 &        10 &           10  &      10 \\
+       grid5recRotate(100,10,-45)&   18 &        10 &           10  &      10 \\
+		gridt(40)        &   58 &        58 &           58  &      58 \\
+		grid9(30)        &   88 &        92 &           104 &      88 \\
+		Smallmesh        &   25 &        13 &           12  &      30 \\
+		Tapir            &   55 &        34 &           18  &      49 \\
+		Eppstein         &   42 &        48 &           45  &      45 \\  \hline \hline
+	\end{tabular}
+	\label{table:bisection}
+\end{table}
+
+
+\begin{table}[h]
+	\caption{Edge-cut results for recursive bi-partitioning.}
+	\centering
+	\begin{tabular}{|l|r|r|r|r|r|} \hline\hline
+		Case            &  Spectral             &  Metis 5.0.2    & Coordinate & Inertial  \\ \hline
+		mesh3e1         &                       &                 &            &           \\
+		airfoil1        &                       &                 &            &           \\
+		3elt            &                       &                 &            &           \\
+		barth4          &                       &                 &            &           \\
+		crack           &                       &                 &            &           \\ \hline \hline
+	\end{tabular}
+	\label{table:Rec_bisection}
+\end{table}
+
+
+\begin{table}[h]
+\caption{Comparing the number of cut edges for recursive bisection and direct multiway partitioning in Metis 5.0.2.}
+\centering
+\begin{tabular}{|l|r|r|r|r|} \hline\hline
+Partitions       &   crack          & airfoil1  \\ \hline
+ 16              &                       &             \\
+ 32              &                       &             \\ \hline \hline
+\end{tabular}
+\label{table:Compare_Metis}
+\end{table}
+
+
+
+\end{document}