HOCKNEY implements the hockney algorithm without 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 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( a,b,c,d,f ) 0002 % HOCKNEY implements the hockney algorithm without 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 arguments: 0009 % a,b,c,d, Specifying diagonal entries 0010 % f, column vector holding the data on the right hand side. 0011 % Output arguments: 0012 % u, vector holding the solution. 0013 0014 [n,M]=size(f); 0015 if M~=1; 0016 error('f must be a column vector'); 0017 end 0018 N=sqrt(n); 0019 0020 % now form the NxN orthogonal matrix Q 0021 Q=zeros(N); 0022 for i=1:N; 0023 for j=1:N; 0024 Q(i,j)=sin(i*j*pi/(N+1)); 0025 end 0026 end 0027 Q=Q*sqrt(2/(N+1)); 0028 0029 % form the N matrix vector products c=Qb. we will hold these in an NxN 0030 % array rather than one long vector as this will make column/row reordering 0031 % more intuitive 0032 C=[]; % initialize array 0033 for k=1:N; 0034 % here we naievely form the products Qf to demonstrate the superiority 0035 % of the FFT 0036 for i=1:N; 0037 temp(i,1)=0; 0038 for j=1:N; 0039 temp(i,1)=temp(i,1)+Q(i,j)*f((k-1)*N+j); 0040 end 0041 end 0042 C=horzcat(C,temp); 0043 end 0044 0045 % rearrange columns to rows, this is as simple as applying a transpose 0046 C=C'; 0047 0048 % we now need to solve the N uncoupled systems 0049 y=[]; % initialise y to an empty array 0050 for k=1:N; 0051 % first form the required TST matrix 0052 T=(a+2*b*cos(k*pi/(N+1)))*diag(ones(N,1))+... 0053 (c+2*d*cos(k*pi/(N+1)))*diag(ones(N-1,1),1)+... 0054 (c+2*d*cos(k*pi/(N+1)))*diag(ones(N-1,1),-1); 0055 0056 L=eye(N);U=T; % initialise L,U 0057 % find LU factorization of T by banded Gaussian elimination 0058 for j=1:N-1; 0059 temp=(U(j+1,j)/U(j,j)); % eliminate elements below the diagonal 0060 U(j+1,:)=U(j+1,:)-U(j,:)*temp; % by subtracting a linear combination 0061 % of rows 0062 0063 L(j+1,j)=temp; % construct L 0064 end 0065 % solve the linear system by back substitution utilising the LU 0066 % factorization 0067 U = U - tril(U,-1); % remove rounding errors 0068 L = L - triu(L,1); 0069 v=Forward(L,C(:,k)); 0070 temp=Backward(U,v); 0071 y=horzcat(y,temp); 0072 end 0073 0074 % rearrange rows to columns 0075 y=y'; 0076 0077 u=[]; % initialise variable to hold the solution 0078 for k=1:N; 0079 temp=Q*y(:,k); 0080 u=vertcat(u,temp); 0081 end 0082 % note that u is a vector, in the context of fast poisson solvers it can be 0083 % interpreted as a natural ordering by columns of the solution on the grid. 0084 end 0085 0086 0087 0088 0089 0090 0091 0092 0093 0094 0095