Prev: find roots of a spline
Next: Reading of Simulink block output value
From: Daniele on 22 Apr 2010 08:51
Hi,
I would like to use the wavelet toolbox to speed-up the solution of sparse linear systems.
The idea is the following: given a linear system Ax=b, decompose A and b in a Wavelet space, solve the linear system in this space and recompose the solution of the linear system.
My code works for a 1 level decomposition, but not for a two level.

%XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
clear all
NN=4;
matrix=rand(2^NN,2^NN)+3*eye(2^NN);
rhs=rand(2^NN,1);
x=matrix\rhs

%XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

for level=1:2
fprintf('level=%d\n',level);

disp('Version 1')
[wavelet_rhs,lung] = wavedec(rhs,level,'db1');
wavelet_matrix=matrix;
for j=1:level
[c,s] = wavedec2(wavelet_matrix,1,'db1');
cad = appcoef2(c,s,'db1',1);
[chd,cvd,cdd] = detcoef2('all',c,s,1);
wavelet_matrix=[cad cvd; chd cdd];
end
wavelet_x=wavelet_matrix\wavelet_rhs;
x_reconstruct=waverec(wavelet_x,lung,'db1')

disp('Version 2')
[wavelet_rhs,lung] = wavedec(rhs,level,'db1');
[c,s] = wavedec2(matrix,level,'db1');
cad = appcoef2(c,s,'db1',1);
[chd,cvd,cdd] = detcoef2('all',c,s,1);
wavelet_matrix=[cad cvd; chd cdd];
wavelet_x=wavelet_matrix\wavelet_rhs;
x_reconstruct=waverec(wavelet_x,lung,'db1')

disp('Version 3')
lung=zeros(3,2);
for i=1:level
[rhs,lung(:,i)] = wavedec(rhs,1,'db1');
end
[c,s] = wavedec2(matrix,level,'db1');
cad = appcoef2(c,s,'db1',1);
[chd,cvd,cdd] = detcoef2('all',c,s,1);
wavelet_matrix=[cad cvd; chd cdd];
xr=wavelet_matrix\rhs;
for i=level:-1:1
xr=waverec(xr,lung(:,i),'db1');
end
x_reconstruct=xr
end
%XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

In the code there are three versions: they perform the same operations and work for the level 1 case (x_reconstruct = x). They do not work for the two level case: what is the correct way to procede in the 2 level case?

Thank you very much.

Daniele
From: Daniele on 26 Apr 2010 07:24
None can help me?

I post a new code
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1- create rhs and matrix
rhs=rand(8,1);
matrix=rand(8,8)+5*eye(8); matrix(3,2)=7;

% 2- solve the linear system
x=matrix\rhs;

% transform into the wavelet space
ris=2;
dwtmode('per','nodisp');

%--- rhs
% [c,lung] = wavedec(rhs,ris,'db5')
% rhs= appcoef(c,lung,'db5',ris);
% for j=ris:-1:1
% d = detcoef(c,lung,j);
% rhs=[rhs; d];
% end
[rhs,lung] = wavedec(rhs,ris,'db5');

%--- matrix
[c,s] = wavedec2(matrix,ris,'db5');
matrix = appcoef2(c,s,'db5',ris);
for j=ris:-1:1
[chd,cvd,cdd] = detcoef2('all',c,s,j);
matrix=[matrix cvd; chd cdd];
end
x_w=matrix\rhs;
x_r=waverec(x_w,lung,'db5');

%--- compare result
abs(x-x_r)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
If the transform into the wavelet space is exact, I'll expect an error of 10e-16, i.e. 0.
However the code works for ris=1 (level 1 transform), not for level=2. Why?

Thanks in advance.

Daniele
From: Daniele on 29 Apr 2010 04:07
Now it works, but only using the 1D routine.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
rhs=rand(2^4,1);
matrix=rand(2^4,2^4)+5*eye(2^4);
x=matrix\rhs;

dwtmode('per');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Ax=b --> PA(P^T) Px = Pb, P=wavelet trasform (P^T P=I)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ris=2
[rhsW,lung] = wavedec(rhs,ris,'db2');
% M=PA
for i=1:2^4
matrixW(:,i)=wavedec(matrix(:,i),ris,'db2');
end
% MP^T=(PM^T)^T
for i=1:2^4
matrixW(i,:)=wavedec(matrixW(i,:),ris,'db2');
end

xW=matrixW\rhsW;
xW = waverec(xW,lung,'db2');

abs(x-xW)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Is it necessary to use the for loops and the 1D routine?
 | 
Pages: 1
Prev: find roots of a spline
Next: Reading of Simulink block output value