utils/matlab/distri.m

% build and distribute matrix of interior,interior
%-------------------------------------------------
% input  : nx,ny,npx,npy,A,as2s,aw2d,cent
% output : sumij

snx = nx/npx;
sny = ny/npy;

ninterf = snx+sny-1;
n1      = (snx-1)*(sny-1);

ka = 1:snx;
kb = snx+1:snx+sny-1;

%------ Cholesky factorization of matrix A
% compute sum( B' inv(A) B )

clear ACH;
clear ACHS;
clear sumij;

sumij = spalloc( ninterf*npx*npy, ninterf*npx*npy, 20*ninterf*npx*npy );
for i = 1:npx,
  ii=(i-1)*snx;
  for j = 1:npy,
%    [i,j]
    jj=(j-1)*sny;

    ad=zeros(snx-1,sny-1,snx-1,sny-1);
    for m2 = 2:sny,
     for l2 = 2:snx,
      for m1 = max(2,m2-1):min(sny,m2+1),
       for l1 = max(2,l2-1):min(snx,l2+1),
        ad(l1-1,m1-1,l2-1,m2-1) = -A(ii+l1+(jj+m1-1)*nx,ii+l2+(jj+m2-1)*nx);
       end
      end
     end
    end

    bi = zeros(snx-1,sny-1,ninterf,npx,npy);
%y=1
    for ix = 2:snx,
      bi(ix-1,   1,      ix,i,j) = -as2d(ii+ix,jj+ 2);
    end
%x=1
    for jy = 2:sny,
      bi(   1,jy-1,jy+snx-1,i,j) = -aw2d(ii+ 2,jj+jy);
    end

%y=1
    jp1 = j+1;
    if jp1 <= npy
     for ix = 2:snx,
       bi( ix-1,sny-1,     ix,i,jp1) = -as2d(ii+ix, (jp1-1)*sny+ 1);
     end
    end
%x=1
    ip1 = 1+mod(i  ,npx);
    for jy = 2:sny,
      bi(snx-1, jy-1,jy+snx-1,ip1,j) = -aw2d((ip1-1)*snx+ 1,jj+jy);
    end
%
    bi = reshape(bi,(snx-1)*(sny-1),ninterf*npx*npy);
%    'bi is there'

    ACHS = chol(sparse( reshape(ad,(snx-1)*(sny-1),(snx-1)*(sny-1)) ));
%    'ACHS is there'

    bi   = sparse(bi' / ACHS);
%    'bi solved'
    sumij = sumij - sparse(bi * bi');
  end
end
clear bi
clear ACHS
clear ad
clear A

%------ distributed matrix of (interface,interface)
qq = zeros(ninterf,npx,npy,ninterf,npx,npy);
for i = 1:npx,
  ii=(i-1)*snx;
  xx=ii+1:ii+snx;
  for j = 1:npy,
    jj=(j-1)*sny;
    yy=jj+2:jj+sny;

    qq(    1,i,j,snx+1,i,j) = -as2d(1+ii,jj+2);
    qq(snx+1,i,j,    1,i,j) = -as2d(1+ii,jj+2);

    for k = 1:snx
      qq(k,i,j,k,i,j) = -cent(ii+k,1+jj);
    end
    for k = 1:snx-1
      qq(k+1,i,j,k,i,j) = -aw2d(ii+k+1,1+jj);
    end
    for k = 2:snx
      qq(k-1,i,j,k,i,j) = -aw2d(ii+k,1+jj);
    end

    for k = 2:sny
      qq(snx-1+k,i,j,snx-1+k,i,j) = -cent(ii+1,jj+k);
    end
    for k = 2:sny-1
      qq(snx-1+k+1,i,j,snx-1+k,i,j) = -as2d(ii+1,jj+k+1);
    end
    for k = 3:sny
      qq(snx-1+k-1,i,j,snx-1+k,i,j) = -as2d(ii+1,jj+k);
    end

    im1 = 1+mod(i-2,npx);
    jm1 = j-1;

    qq(1,i,j,snx,im1,j) = -aw2d( ii+1,1+jj);
    if jm1 > 0
      qq(1,i,j,snx+sny-1,i,jm1) = -as2d( ii+1,1+jj);
    end

    ip1 = 1+mod(i  ,npx);
    jp1 = j+1;
    iif=(ip1-1)*snx;
    jjf=(jp1-1)*sny;

    qq(snx,i,j,1,ip1,j) = -aw2d(iif+1,1+jj);
    if jp1 <= npy
      qq(snx+sny-1,i,j,1,i,jp1) = -as2d(ii+1,1+jjf);
    end
  end
end

%------ compute matrix S = Q - sum(B' inv(A) B)
sumij = sumij + sparse(reshape(qq,ninterf*npx*npy,ninterf*npx*npy));

clear qq
1	% build and distribute matrix of interior,interior
2	%-------------------------------------------------
3	% input : nx,ny,npx,npy,A,as2s,aw2d,cent
4	% output : sumij
5
6	snx = nx/npx;
7	sny = ny/npy;
8
9	ninterf = snx+sny-1;
10	n1 = (snx-1)*(sny-1);
11
12	ka = 1:snx;
13	kb = snx+1:snx+sny-1;
14
15	%------ Cholesky factorization of matrix A
16	% compute sum( B' inv(A) B )
17
18	clear ACH;
19	clear ACHS;
20	clear sumij;
21
22	sumij = spalloc( ninterfnpxnpy, ninterfnpxnpy, 20ninterfnpx*npy );
23	for i = 1:npx,
24	ii=(i-1)*snx;
25	for j = 1:npy,
26	% [i,j]
27	jj=(j-1)*sny;
28
29	ad=zeros(snx-1,sny-1,snx-1,sny-1);
30	for m2 = 2:sny,
31	for l2 = 2:snx,
32	for m1 = max(2,m2-1):min(sny,m2+1),
33	for l1 = max(2,l2-1):min(snx,l2+1),
34	ad(l1-1,m1-1,l2-1,m2-1) = -A(ii+l1+(jj+m1-1)nx,ii+l2+(jj+m2-1)nx);
35	end
36	end
37	end
38	end
39
40	bi = zeros(snx-1,sny-1,ninterf,npx,npy);
41	%y=1
42	for ix = 2:snx,
43	bi(ix-1, 1, ix,i,j) = -as2d(ii+ix,jj+ 2);
44	end
45	%x=1
46	for jy = 2:sny,
47	bi( 1,jy-1,jy+snx-1,i,j) = -aw2d(ii+ 2,jj+jy);
48	end
49
50	%y=1
51	jp1 = j+1;
52	if jp1 <= npy
53	for ix = 2:snx,
54	bi( ix-1,sny-1, ix,i,jp1) = -as2d(ii+ix, (jp1-1)*sny+ 1);
55	end
56	end
57	%x=1
58	ip1 = 1+mod(i ,npx);
59	for jy = 2:sny,
60	bi(snx-1, jy-1,jy+snx-1,ip1,j) = -aw2d((ip1-1)*snx+ 1,jj+jy);
61	end
62	%
63	bi = reshape(bi,(snx-1)(sny-1),ninterfnpx*npy);
64	% 'bi is there'
65
66	ACHS = chol(sparse( reshape(ad,(snx-1)(sny-1),(snx-1)(sny-1)) ));
67	% 'ACHS is there'
68
69	bi = sparse(bi' / ACHS);
70	% 'bi solved'
71	sumij = sumij - sparse(bi * bi');
72	end
73	end
74	clear bi
75	clear ACHS
76	clear ad
77	clear A
78
79	%------ distributed matrix of (interface,interface)
80	qq = zeros(ninterf,npx,npy,ninterf,npx,npy);
81	for i = 1:npx,
82	ii=(i-1)*snx;
83	xx=ii+1:ii+snx;
84	for j = 1:npy,
85	jj=(j-1)*sny;
86	yy=jj+2:jj+sny;
87
88	qq( 1,i,j,snx+1,i,j) = -as2d(1+ii,jj+2);
89	qq(snx+1,i,j, 1,i,j) = -as2d(1+ii,jj+2);
90
91	for k = 1:snx
92	qq(k,i,j,k,i,j) = -cent(ii+k,1+jj);
93	end
94	for k = 1:snx-1
95	qq(k+1,i,j,k,i,j) = -aw2d(ii+k+1,1+jj);
96	end
97	for k = 2:snx
98	qq(k-1,i,j,k,i,j) = -aw2d(ii+k,1+jj);
99	end
100
101	for k = 2:sny
102	qq(snx-1+k,i,j,snx-1+k,i,j) = -cent(ii+1,jj+k);
103	end
104	for k = 2:sny-1
105	qq(snx-1+k+1,i,j,snx-1+k,i,j) = -as2d(ii+1,jj+k+1);
106	end
107	for k = 3:sny
108	qq(snx-1+k-1,i,j,snx-1+k,i,j) = -as2d(ii+1,jj+k);
109	end
110
111	im1 = 1+mod(i-2,npx);
112	jm1 = j-1;
113
114	qq(1,i,j,snx,im1,j) = -aw2d( ii+1,1+jj);
115	if jm1 > 0
116	qq(1,i,j,snx+sny-1,i,jm1) = -as2d( ii+1,1+jj);
117	end
118
119	ip1 = 1+mod(i ,npx);
120	jp1 = j+1;
121	iif=(ip1-1)*snx;
122	jjf=(jp1-1)*sny;
123
124	qq(snx,i,j,1,ip1,j) = -aw2d(iif+1,1+jj);
125	if jp1 <= npy
126	qq(snx+sny-1,i,j,1,i,jp1) = -as2d(ii+1,1+jjf);
127	end
128	end
129	end
130
131	%------ compute matrix S = Q - sum(B' inv(A) B)
132	sumij = sumij + sparse(reshape(qq,ninterfnpxnpy,ninterfnpxnpy));
133
134	clear qq