Home > K25104 > PDEs > hockney.m

hockney

PURPOSE ^

HOCKNEY implements the hockney algorithm without FFT to solve the block

SYNOPSIS ^

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

DESCRIPTION ^

 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.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

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