%---------- HB matrix example for 1D Laplacian --- Hierachical Basis (HB)
% The Christian Wagner example (p.25, 1999) fully illustrated:
%        www.iwr.uni-heidelberg.de/iwr/techsim/chris
help HB1
I=ones(4); B=band(I,1,0); B=-B(:,1:3);
Agiven=[2*eye(4) B; B' 2*eye(3)];  % Wagner matrix p.25
 HB200 = [1 5 2 6 3 7 4];
 A_Lexico  = Agiven(HB200,HB200)  

 HB2LX = [1 5 2 7 3 6 4];
for j=1:6, fprintf('[%d]---',LX2HB(j)); end,fprintf('[%d] HB2LEXICO\n',LX2HB(7));
 LX2HB = [1 3 5 7 2 6 4];
 A = A_Lexico(LX2HB,LX2HB)  % Permute to the 2-level HB ordering (C -> F)
disp('Step 1: finest std FEM (1 2 3 4 5 6 7) --> HB level 1 (HB 1 2 3 4 FEM 5 6 7)')

for j=1:6,fprintf('[%d]---',j); end,  fprintf('[%d] Lexicographically\n',7);
for j=1:6,fprintf('[%d]---',HB2LX(j)); end,fprintf('[%d] HB ordering\n',HB2LX(7));
   disp('     Recall:  \phi_C = \phi_f + 0.5*\phi_neighbours (finest level)')
J1=[eye(4) -B/2; 0*B' eye(3)]; %  % First HB transform

 o=[1:5 7 6]; J1=J1(:,o); J1(6:7,6:7)=eye(2);
A1=J1'*A*J1;
mat_prt(J1,[],' ','J1 to work out coarse \phi at 1 2 3 4')
disp('Pause ...');pause
mat_prt(A1,[],' ','Result A1=J1'' A J1 (Next coarse 5 6)')
disp('Step 2: HB level 1 (HB 1 2 3 4 FEM 5 6 7) --> HB basis (HB 1 2 3 4 5 6 FEM 7)')
disp('Pause ...');pause

J2=eye(7); J2(5,7)=1/2; J2(6,7)=1/2; % Second HB transform
for j=1:6
if HB2LX(j)>4, fprintf('[%d]---',HB2LX(j)); else,
               fprintf('[%s]---',' '); end, end, fprintf('[ ] HB ordering\n')
   disp('     Recall:  \phi_C = \phi_f + 0.5*\phi_neighbour (fine level)')
A2=J2'*A1*J2;  %------------- Final diagonal matrix in HB basis
mat_prt(J2,[],' ','J2 to work out coarse \phi at 5 6')
mat_prt(A2,[],' ','Result A2=J2'' A1 J2 (Final HB matrix)')
disp('Pause ...');pause

J = J1*J2;

disp('1D Laplician (std FEM to HB FEM)========================= Display ======')
mat_prt(A_Lexico,A,'|','Original Matrix A',...
                        'Original in HB ordering (Coarse to Fine)')
mat_prt(J1,A1,':',      'The First HB transform J1', ...
                        'The First HB transformed A1=J1''*A*J1')
mat_prt(J2,J,'|',      'The Second HB transform J2', ...
                         'The HB transform J=J1*J2 for A2=J''*A*J')
mat_prt(A2,r,'+', 'The Second HB transformed A2=J2''*A1*J2', ...
                        'The HB nodal ordering (integers of course!)')
disp('HB1.m==================================================== Display ======')
