femwind/prol_rest_err.m

   1 function err=prol_rest_err(hzc,icl3,X,params)
   2 % err=prol_rest_err(hzc,icl3,X,params)
   3 % test prolongation and restriction
   4 % err=0,return
   5 n = size(X{1});
   6 [icl1,icl2]=hzc2icl(hzc,n);
   7 nc = [length(icl1),length(icl2),length(icl3)];
   8 P =coarsening_P(hzc,icl3,X,params);
   9 uc = rand(nc);
  10 u0 = reshape(P*uc(:),n);
  11 u1 = prolongation(uc,hzc,icl3,X,params);
  12 err1 = big(u0-u1);
  13 tol = eps(single(1))*15*big(u0);
  14 if err1 > tol
  15     err1,tol
  16     warning('prol_restr_err: prolongation error vs. matrix P too large')
  17 end
  18 if exist('fortran/prolongation_test.exe') && isfield(params,'test_fortran') && params.test_fortran
  19     disp('testing if prolongation same result in fortran')
  20     u2 = prolongation_fortran(uc,hzc,icl3,X,params);
  21     err2 = big(u0-u2);
  22     if err2 > tol
  23     err2,tol
  24         warning('prol_restr_err: prolongation error fortran too large')
  25     end
  26 else
  27     err2=0;
  28 end
  29
  30 u = rand(n);
  31 uc0 = reshape(P'*u(:),nc);
  32 uc1 = restriction(u,hzc,icl3,X,params);
  33 errc = big(uc0-uc1);
  34 if errc > tol
  35     errc,tol
  36     warning('prol_restr_err: restriction error vs. matrix P too large')
  37 end
  38 if exist('fortran/restriction_test.exe')  && isfield(params,'test_fortran') && params.test_fortran
  39     disp('testing if restriction same result in fortran')
  40     uc2 = restriction_fortran(u,hzc,icl3,X,params);
  41     errc2 = big(uc0-uc2);
  42     if errc2 > tol
  43     errc2,tol
  44         warning('prol_restr_err: prolongation error fortran too large')
  45     end
  46 else
  47     errc2=0;
  48 end
  49 err=max([err1,err2,errc,errc2]);
  50 if err>tol,
  51     error('prol_rest_err: error too large')
  52 else
  53     fprintf('prol_rest_err %g OK\n',err)
  54 end
  55
  56 end