%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% STAR-P Demo for 18.337 Lecture 2, Spring 2004
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

np % Number of processes

% Big matrix multiplications
n=500*p
a=randn(n)
b=randn(n)
tic, c=a*b; toc

n=n*2, a=randn(n); b=randn(n);
tic, c=a*b; toc

n=n*2, a=randn(n); b=randn(n);
tic, c=a*b; toc

n=n*2, a=randn(n); b=randn(n);
tic, c=a*b; toc

% Big matrix inverse
n=500*p
a=randn(n); tic, b=inv(a); toc
n=n*2, a=randn(n); tic, b=inv(a); toc
n=n*2, a=randn(n); tic, b=inv(a); toc

e=eye(n)
c=a*b-e
norm(c,inf)

% Hilbert matrix
type hilb
hilb(3)
hilb(7)
c=hilb(6*p)
c(:,:)
c(:,:)-hilb(6)
c3 = c*c*c;
c3(:,:)-hilb(6)^3

% Big dense linear system
n=2000*p;
b=randn(n,1)
a=randn(n)
c=a\b;
residual = a*c-b
norm(residual)

%%%%%%%%%%%%%%%%%%%%%%%
% MM-Mode - Compute pi
%%%%%%%%%%%%%%%%%%%%%%%

% Method: atan(1)*4
format long
atan(1)*4

% Use MATLAB quadrature
quadl('4./(1+x.^2)',0,1,1e-14)

a = (0:(np-1)*p)/np
a(:)
b = a + (1/np)
b(:)
mypi = mm('quadl','4./(1+x.^2)',a,b,1e-14);
mypi(:)

mypi=sum(mypi)