femwind/rb_schwarz_solve.m

   1 function [x,it]=rb_schwarz_solve(K,F,X)
   2 % x=rb_line_gs_solve(K,F,X)
   3 disp('solver red-black horizonal schwarz')
   4 n = size(X{1})
   5 nn = size(F,1)
   6 if nn~=prod(n)
   7     error('rb_line_gs_solve: inconsistent sizes')
   8 end
   9
  10 % parameters
  11 maxit=1000
  12 tol = 1e-5*big(F)/big(K)
  13 sub_size=[6,6]
  14 sub_step=[4,4]
  15 overlap=sub_size-sub_step
  16 % index sets
  17
  18 % divide into intervals in dimensions 1 and 2
  19 sub_nodes=cell(1,2);
  20 for i=1:2
  21     msub=ceil(n(i)/sub_step(i));
  22     % sub_nodes{i}=cell(msub,1);
  23     for j=1:msub
  24         from = (j-1)*sub_step(i)+1;
  25         to = min(from+sub_size(i)-1,n(i));
  26         sub_nodes{i}{j}=[from:to];
  27         disp(['direction ',num2str(i),' group ',num2str(j),' indices ',num2str(sub_nodes{i}{j})])
  28         if to >= n(i)
  29             nsub(i)=j;
  30             break
  31         end
  32      end
  33 end
  34
  35 check_cover(n(1),sub_nodes{1}(:))
  36 check_cover(n(2),sub_nodes{2}(:))
  37
  38
  39 % global node indices for the subdomains
  40 IX = cell(nsub);
  41 for i1=1:nsub(1)
  42     for i2=1:nsub(2)
  43         idx1 = sub_nodes{1}{i1}; % horizontal rectangle
  44         idx2 = sub_nodes{2}{i2};
  45         idx3 = 1:n(3); % vertical all
  46         disp(['subdomain ',num2str([i1,i2]),' nodes ',num2str(idx1),' by ',num2str(idx2)])
  47         [ix1,ix2,ix3]=ndgrid(idx1,idx2,idx3);
  48         ix=sub2ind(n,ix1(:),ix2(:),ix3(:));
  49         IX{i1,i2}=ix;
  50     end
  51 end
  52
  53 % check indices
  54 check_cover(nn,IX(:))
  55
  56 disp('decomposing subdomain matrices')
  57 R=cell(nsub);  % cholesky factor
  58 P=cell(nsub);  % permutation
  59 for i1=1:nsub(1)
  60     for i2=1:nsub(2)
  61         idx1 = sub_nodes{1}{i1}; % horizontal rectangle
  62         idx2 = sub_nodes{2}{i2};
  63         idx3 = 1:n(3); % vertical all
  64         ix = IX{i1,i2};
  65         [R{i1,i2},ierr,P{i1,i2}]=chol(K(ix,ix),'vector');
  66         if ierr
  67             error(['Subdomain matrix is not positive definite, error ',num2str(ierr)])
  68         end
  69         p = P{i1,i2}; r = R{i1,i2}; k = K(ix,ix); err(1) = big(k(p,p)-r'*r);
  70         b = rand(length(p),1); x = r\(r'\b(p)); x(p)=x; err(2) = big(k*x-b);
  71         fprintf('subdomain %g %g size %g %g %g nodes %g nonzeros %g Cholesky %g error %g %g\n',...
  72             i1,i2,length(idx1),length(idx2),length(idx3),length(ix),...
  73             nnz(K(ix,ix)),nnz(R{i1,i2}),err)
  74     end
  75 end
  76
  77 x = zeros(nn,1);
  78 % xx = zeros(nn,1);
  79 for it=1:maxit
  80     for rb1=1:2
  81         for rb2=1:2
  82             for i1=rb1:2:nsub(1)
  83                 for i2=rb2:2:nsub(2)
  84                     % solving horizontal location i1 i2 and vertical line
  85                     ix = IX{i1,i2};
  86                     % x(ix) = x(ix) + K(ix,ix)\(F(ix) - K(:,ix)'*x);
  87                     p = P{i1,i2}; % Ksub(p,p)= R'*R
  88                     res = F(ix) - K(:,ix)'*x;
  89                     sol = R{i1,i2}\(R{i1,i2}'\res(p));
  90                     sol(p) = sol;
  91                     % err = big(res - K(ix,ix)*sol)
  92                     x(ix) = x(ix) + sol;
  93                     % err = big(xx(ix)-x(ix))
  94                 end
  95             end
  96         end
  97     end
  98     res= norm(K*x-F);
  99     fprintf('iteration %g residual %g tolerance %g\n',it,res,tol)
 100     if res<tol,
 101         break
 102     end
 103 end
 104
 105 end
 106
 107 function check_cover(m,x)
 108 % check if cell array x contains indices that cover 1 to m
 109 c=zeros(m,1);
 110 for i=1:length(x)
 111    c(x{i})=c(x{i})+1;
 112 end
 113 if any(c==0),
 114     error('subdomains do not cover all nodes')
 115 end
 116 end
 117