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.
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