Home > K25104 > PDEs > hockney_FFT.m

hockney_FFT

PURPOSE ^

Implements the Hockney algorithm with FFT to solve the block

SYNOPSIS ^

function [ u ] = hockney_FFT( a,b,c,d,f )

DESCRIPTION ^

 Implements the Hockney algorithm with FFT to solve the block 
 TST system Au=f.  
 Input arguments a and b prescribe the diagonal and off diagonal 
 terms respectively of NXN TST matrix B, c and d prescribe the diagonal  
 and off diagonal terms respectively of NXN TST matrix C. A is then the 
 (N^2)X(N^2) block TST matrix with B on the diagonal and C off diagonal.
 Input argument f is a column vector holding the data on the right hand
 side.
 Input arguments:
   a,b,c,d, Specifying diagonal entries
   f, column vector holding the data on the right hand side. 
 Output arguments:
   u, vector holding the solution.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [ u ] = hockney_FFT( a,b,c,d,f )
0002 % Implements the Hockney algorithm with FFT to solve the block
0003 % TST system Au=f.
0004 % Input arguments a and b prescribe the diagonal and off diagonal
0005 % terms respectively of NXN TST matrix B, c and d prescribe the diagonal
0006 % and off diagonal terms respectively of NXN TST matrix C. A is then the
0007 % (N^2)X(N^2) block TST matrix with B on the diagonal and C off diagonal.
0008 % Input argument f is a column vector holding the data on the right hand
0009 % side.
0010 % Input arguments:
0011 %   a,b,c,d, Specifying diagonal entries
0012 %   f, column vector holding the data on the right hand side.
0013 % Output arguments:
0014 %   u, vector holding the solution.
0015 
0016 [n,M]=size(f);
0017 if M~=1;
0018     error('f must be a column vector');
0019 end
0020 N=sqrt(n);
0021 
0022 % form the N matrix vector products c=Qb using FFT. we will hold these in an NxN
0023 % array rather than one long vector as this will make column/row reordering
0024 % more intuitive
0025 C=[];    % initialize array
0026 for k=1:N;
0027     u=f((k-1)*N+1:k*N)';   
0028     u=horzcat(0,u,zeros(1,N+1));  % extend the vector to make it consistent
0029                                   % with the Fast Fourier Transform
0030     temp=inv_FFT(u);              % apply FFT
0031     temp=temp(2:N+1);             % cut off unwanted parts
0032     temp=sqrt(2/(N+1))*imag(temp);  
0033     temp=temp';
0034     C=horzcat(C,temp);
0035     clear temp;
0036 end
0037 
0038 % rearrange columns to rows, this is as simple as applying a traspose
0039 C=C';
0040 
0041 % we now need to solve the N uncoupled systems
0042 y=[];   % initialize y to an empty array
0043 for k=1:N;
0044     % first form the required TST matrix
0045     T=(a+2*b*cos(k*pi/(N+1)))*diag(ones(N,1))+...
0046         (c+2*d*cos(k*pi/(N+1)))*diag(ones(N-1,1),1)+...
0047         (c+2*d*cos(k*pi/(N+1)))*diag(ones(N-1,1),-1);
0048     
0049     L=eye(N);U=T;   % initialize L,U
0050     % find LU factorization of T by banded Gaussian elimination
0051     for j=1:N-1;
0052         temp=(U(j+1,j)/U(j,j));       % eliminate elements below the diagonal
0053         U(j+1,:)=U(j+1,:)-U(j,:)*temp;  % by subtracting a linear combination
0054                                     % of rows
0055                                      
0056         L(j+1,j)=temp;        % construct L
0057     end 
0058     % solve the linear system by back substitution utilizing the LU
0059     % factorization
0060     U = U - tril(U,-1); % remove rounding errors
0061     L = L - triu(L,1);
0062     v=Forward(L,C(:,k));   
0063     temp=Backward(U,v);    
0064     y=horzcat(y,temp);
0065 end
0066 
0067 
0068 u=[];   % initialize array to hold solution
0069 for k=1:N;
0070     in=y(k,:);   
0071     in=horzcat(0,in,zeros(1,N+1));  % extend the vector to make it consistent
0072                                   % with the Fast Fourier Transform
0073     temp=inv_FFT(in);             % apply FFT
0074     temp=temp(2:N+1);             % cut off unwanted parts
0075     temp=sqrt(2/(N+1))*imag(temp);  
0076     temp=temp';
0077     u=vertcat(u,temp);
0078     clear temp;
0079 end
0080 % note that u is a vector, in the context of fast Poisson solvers it can be
0081 % interpreted as a natural ordering by columns of the solution on the grid.
0082 end

Generated on Mon 18-Jan-2016 10:25:49 by m2html © 2005