File: diagainv_example.m

package info (click to toggle)
mumps 5.1.2-5
  • links: PTS, VCS
  • area: main
  • in suites: sid
  • size: 15,704 kB
  • sloc: fortran: 310,672; ansic: 12,364; xml: 521; makefile: 469
file content (58 lines) | stat: -rw-r--r-- 1,272 bytes parent folder | download | duplicates (3)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
%Example of using MUMPS in matlab to compute diagonal of inverse of A

% Change to true to test example complex arithmetic
complex_arithmetic = false;

% initialization of a matlab MUMPS structure
id = initmumps;
if (complex_arithmetic)
  id = zmumps(id);
else
  id = dmumps(id);
end
load lhr01;
mat = Problem.A;

if (complex_arithmetic)
  % To test complex version
  mat = mat + i * speye(size(mat,1),size(mat,1)); 
end

% JOB = 4 means analysis+factorization
id.JOB = 4;
if (complex_arithmetic)
  id = zmumps(id,mat);
else
  id = dmumps(id,mat);
end

% Set the right hand side structure to requested entries of A-1
id.RHS = speye(size(mat,1),size(mat,1)); % Sparse format required

%call MUMPS solution phase to compute diagonal entries of A-1
id.ICNTL(30)=1; % Ask for A-1 entries
id.JOB=3;


if (complex_arithmetic)
  id = zmumps(id,mat);
else
  id = dmumps(id,mat);
end

% diagonal values have been computed in
% the (sparse) matrix id.SOL, which has
% the same structure as id.RHS

% Compare diagonal of inverse computed by Mumps and by matlab
disp(' ');
disp('Computing 2-norm of error on diagonal of inverse:');
norm(diag( diag(diag(inv(mat)))-id.SOL ),2)

% destroy mumps instance
id.JOB = -2;
if (complex_arithmetic)
  id = zmumps(id)
else
  id = dmumps(id)
end