examplecomplex

 This is an example script demonstrating how PARDISO works on a fairly
 large sparse, symmetric and complex matrix.

0001 % This is an example script demonstrating how PARDISO works on a fairly
0002 % large sparse, symmetric and complex matrix.
0003 clear
0004 
0005 dkepath = load_structures_yp('dkepath','','');
0006 
0007 % Script parameters.
0008 % -----------------
0009 verbose = false;
0010 l       = 2;
0011 k       = l * pi;
0012 r0      = 0.5;
0013 R       = 2;
0014 M       = 15 * l;
0015 N       = 120 * l; 
0016 
0017 fprintf('Grid size: %d x %d\n',M,N);
0018 
0019 % Construct the complex matrix C.
0020 % ------------------------------
0021 n  = M*N; 
0022 r  = linspace(r0,R,M);
0023 dr = (R - r0)/(M - 1);
0024 
0025 dtheta = 2 * pi/N;
0026 theta  = 0 : dtheta : (2*pi - dtheta);
0027 A      = sparse([],[],[],n,n,5*n);
0028 b      = zeros(n,1);
0029 
0030 ratio = ones(M-1,1);
0031 for j = 2:M-1
0032   alpha0 = 1/(dr^2) - 1/(2*dr*r(j));
0033   alpha1 = 1/(r(j)*dtheta)^2;
0034   alpha2 = k^2 - 2/(dr^2) - 2/(r(j)^2 * dtheta^2);
0035 
0036   if j > 2
0037     ratio(j-1) = ratio(j-2) * alpha0 / alpha3;
0038   end
0039 
0040   alpha3 = 1/(dr^2) + 1/(2*dr*r(j));
0041 
0042   for j2 = (j-2)*N+1:(j-1)*N
0043      A(j2+N,j2) = alpha0;
0044   end
0045   
0046   j2 = (j-1)*N + 1;
0047   A(j2:j2+1,j2) = [ alpha2
0048                     alpha1 ];
0049 
0050   for j2 = (j-1)*N+2:j*N-1        
0051     A(j2-1:j2+1,j2) = [ alpha1
0052                         alpha2
0053                         alpha1 ];
0054   end                    
0055   
0056   j2=j*N;                
0057   A(j2-1:j2,j2) = [ alpha1
0058                     alpha2 ];    
0059   
0060   for j2 = j*N+1:(j+1)*N
0061     A(j2-N,j2) = alpha3;
0062   end
0063   
0064   A((j-1)*N+1,j*N) = alpha1; 
0065   A(j*N,(j-1)*N+1) = alpha1; 
0066 end
0067 
0068 for j = 1:N
0069     A(j,j)=1;
0070 end
0071 
0072 for j = 1:N
0073   b(j) = -exp(i*k*r(1)*cos(theta(j)));
0074 end
0075 
0076 beta0      = 2/(dr^2);
0077 beta1      = 1/(r(end)*dtheta)^2;
0078 beta2      = k^2 - 2/(dr^2) - 2/(r(end)*dtheta)^2 + ...
0079              2*i*k*dr*(1/dr^2 + 1/(2*dr*r(end)));
0080 ratio(M-1) = ratio(M-2) * beta0 / alpha3;
0081     
0082 for j = (M-2)*N+1:(M-1)*N
0083   A(j+N,j) = beta0;
0084 end
0085 
0086 j2            = (M-1)*N+1;      
0087 A(j2:j2+1,j2) = [ beta2
0088                   beta1 ];    
0089 
0090 for j2 = (M-1)*N+2:M*N-1        
0091   A(j2-1:j2+1,j2) = [ beta1
0092                       beta2
0093                       beta1 ];
0094 end        
0095 
0096 j2                = M*N;
0097 A(j2-1:j2,j2)     = [ beta1
0098                       beta2 ];
0099 A((M-1)*N+1, M*N) = beta1;
0100 A(M*N,(M-1)*N+1)  = beta1;
0101 
0102 B        = A;
0103 B(1:N,:) = [];
0104 B(:,1:N) = [];
0105 c        = b;
0106 c(1:N)   = [];
0107 
0108 for j = 1:N
0109   c(1:N) = c(1:N) - A(N+1:2*N,j) * b(j);
0110 end
0111 
0112 dL = sqrt(abs(1./ratio));
0113 dR = sqrt(abs(ratio)) .* sign(ratio);
0114 DL = kron(spdiags(dL,0,M-1,M-1),speye(N));
0115 DR = kron(spdiags(dR,0,M-1,M-1),speye(N));
0116 C  = DL * B * DR;
0117 C  = (C + C.')/2;
0118 c  = DL * c;
0119 
0120 % Solve for v in linear system C * v = c.
0121 % --------------------------------------
0122 % Initialize the PARDISO internal data structures. We've told PARDISO to
0123 % handle complex matrices using the sparse direct solver.
0124 info = pardisoinit(6,0);
0125 
0126 % Analyze the matrix and compute a symbolic factorization.
0127 info = pardisoreorder(tril(C),info,verbose);
0128 fprintf('The factors have %d nonzero entries.\n',info.iparm(18));
0129 
0130 % Compute the numeric factorization.
0131 info = pardisofactor(tril(C),info,verbose);
0132 
0133 % Compute the solution v using the symbolic factorization.
0134 [v info] = pardisosolve(tril(C),c,info,verbose);
0135 fprintf('PARDISO performed %d iterative refinement steps.\n',info.iparm(7));
0136 
0137 % Compute the residuals.
0138 r = abs(C*v - c);
0139 fprintf('The maximum residual for the solution is %0.3g.\n',max(r));
0140 
0141 % Free the PARDISO data structures.
0142 pardisofree(info);
0143 clear info