nag_sparse_complex_herm_precon_ilu_solve (f11jp) solves a system of complex linear equations involving the incomplete Cholesky preconditioning matrix generated by
nag_sparse_complex_herm_precon_ilu (f11jn).
nag_sparse_complex_herm_precon_ilu_solve (f11jp) solves a system of linear equations
involving the preconditioning matrix
, corresponding to an incomplete Cholesky decomposition of a complex sparse Hermitian matrix stored in symmetric coordinate storage (SCS) format (see
Symmetric coordinate storage (SCS) format in the F11 Chapter Introduction), as generated by
nag_sparse_complex_herm_precon_ilu (f11jn).
In the above decomposition
is a complex lower triangular sparse matrix with unit diagonal,
is a real diagonal matrix and
is a permutation matrix.
and
are supplied to
nag_sparse_complex_herm_precon_ilu_solve (f11jp) through the matrix
which is a lower triangular
by
complex sparse matrix, stored in SCS format, as returned by
nag_sparse_complex_herm_precon_ilu (f11jn). The permutation matrix
is returned from
nag_sparse_complex_herm_precon_ilu (f11jn) via the array
ipiv.
nag_sparse_complex_herm_precon_ilu_solve (f11jp) may also be used in combination with
nag_sparse_complex_herm_precon_ilu (f11jn) to solve a sparse complex Hermitian positive definite system of linear equations directly (see
nag_sparse_complex_herm_precon_ilu (f11jn)). This is illustrated in
Example.
None.
The computed solution
is the exact solution of a perturbed system of equations
, where
is a modest linear function of
, and
is the
machine precision.
The time taken for a call to
nag_sparse_complex_herm_precon_ilu_solve (f11jp) is proportional to the value of
nnzc returned from
nag_sparse_complex_herm_precon_ilu (f11jn).
This example reads in a complex sparse Hermitian positive definite matrix
and a vector
. It then calls
nag_sparse_complex_herm_precon_ilu (f11jn), with
and
, to compute the
complete Cholesky decomposition of
:
Finally it calls
nag_sparse_complex_herm_precon_ilu_solve (f11jp) to solve the system
function f11jp_example
fprintf('f11jp example results\n\n');
n = int64(9);
nz = int64(23);
a = zeros(3*nz,1);
irow = zeros(3*nz, 1, 'int64');
icol = zeros(3*nz, 1, 'int64');
a(1:nz) = [ 6 + 0.i; -1 + 1.i; 6 + 0.i; 0 + 1.i;
5 + 0.i; 5 + 0.i; 2 - 2.i; 4 + 0.i;
1 + 1.i; 2 + 0.i; 6 + 0.i; -4 + 3.i;
0 + 1.i; -1 + 0.i; 6 + 0.i; -1 - 1.i;
0 - 1.i; 9 + 0.i; 1 + 3.i; 1 + 2.i;
-1 + 0.i; 1 + 4.i; 9 + 0.i];
b = [ 8 + 54i;-10 - 92i; 25 + 27i; 26 - 28i;
54 + 12i; 26 - 22i; 47 + 65i; 71 - 57i;
60 + 70i];
irow(1:nz) = int64([1;2;2;3;3;4;5;5;6;6;6;7;7;7;7;8;8;8;9;9;9;9;9]);
icol(1:nz) = int64([1;1;2;2;3;4;1;5;3;4;6;2;5;6;7;4;6;8;1;5;6;8;9]);
lfill = int64(0);
dtol = 0;
mic = 'N';
dscale = 0;
ipiv = zeros(n, 1, 'int64');
[a, irow, icol, ipiv, istr, nnzc, npivm, ifail] = ...
f11jn( ...
nz, a, irow, icol, lfill, dtol, mic, dscale, ipiv);
method = 'CG ';
precon = 'Preconditioned';
tol = (x02aj)^(3/8);
maxitn = int64(20);
anorm = 0;
sigmax = 0;
maxits = int64(9);
monit = int64(2);
[lwreq, work, ifail] = ...
f11gr( ...
method, precon, int64(n), tol, maxitn, anorm, sigmax, ...
maxits, monit, 'sigcmp', 's', 'norm_p', '1');
irevcm = int64(0);
u = complex(zeros(n,1));
v = b;
wgt = zeros(n,1);
while (irevcm ~= 4)
[irevcm, u, v, work, ifail] = ...
f11gs( ...
irevcm, u, v, wgt, work);
if (irevcm == 1)
[v, ifail] = f11xs( ...
a(1:nz), irow(1:nz), icol(1:nz), 'N', u);
elseif (irevcm == 2)
[v, ifail] = f11jp( ...
a, irow, icol, ipiv, istr, 'N', u);
elseif (irevcm == 3)
[itn, stplhs, stprhs, anorm, sigmax, its, sigerr, ifail] = ...
f11gt(work);
fprintf('\nMonitoring at iteration number %2d\n',itn);
fprintf('residual norm: %14.4e\n', stplhs);
fprintf('\n Solution Vector\n');
disp(u);
fprintf('\n Residual Vector\n');
disp(v);
end
end
[itn, stplhs, stprhs, anorm, sigmax, its, sigerr, ifail] = ...
f11gt(work);
fprintf('\nNumber of iterations for convergence: %4d\n', itn);
fprintf('Residual norm: %14.4e\n', stplhs);
fprintf('Right-hand side of termination criteria: %14.4e\n', stprhs);
fprintf('i-norm of matrix a: %14.4e\n', anorm);
fprintf('\n Solution Vector\n');
disp(u);
fprintf('\n Residual Vector\n');
disp(v);