File: schur_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 (92 lines) | stat: -rw-r--r-- 2,188 bytes parent folder | download | duplicates (6)
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
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
%Example of using MUMPS in matlab with schur option

% initialization of a matlab MUMPS structure
id = initmumps;
id = dmumps(id);
load lhr01;
mat = Problem.A;
themax = max(max(abs(mat)));
n = size(mat,1);
mat = mat+sparse(1:n,1:n,3*themax*ones(n,1));

% initialization of Schur option
id.VAR_SCHUR = [n-9:n];

% JOB = 6 means analysis+facto+solve
id.JOB = 6;
id.RHS = ones(size(mat,1),1);
%call to mumps
id = dmumps(id,mat);
disp('*** check solution restricted to mat(1:n-10,1:n-10)');
if(norm(mat(1:n-10,1:n-10)*id.SOL(1:n-10) - ones(n-10,1),'inf') > sqrt(eps))
	disp('WARNING : precision may not be OK');
else
	disp('SCHUR SOLUTION CHECK1 OK');
end
norm(mat(1:n-10,1:n-10)*id.SOL(1:n-10) - ones(n-10,1),'inf')


% we want to use Schur complement to solve 
% A * sol = rhs
% with sol = x   and rhs = rhs1
%            y             rhs2
%
% check that the complete solution verify
% y = S^(-1) * (rhs2 - A_{2,1} * A_{1,1}^(-1) * rhs1)
% and 
% x = A_{1,1}^(-1) * rhs1) - A_{1,2} * y
%
sol1 = id.SOL(1:n-10);
rhsy = ones(10,1)-mat(n-9:n,1:n-10)*sol1;

%%%%%%%%%%%%%%%%%%%
% TO CHANGE :
% usually the resolution below is replaced by an iterative scheme
y = id.SCHUR \ rhsy;
%%%%%%%%%%%%%%%%%%%%

rhsx = mat(1:n-10,n-9:n)*y;
id.JOB = 3;
id.RHS(1:n-10) = rhsx;
id = dmumps(id,mat);
rhsx = id.SOL(1:n-10);
x = sol1-rhsx;
sol = [x;y];
r = mat*sol - ones(n,1);
disp('*** check complete solution');
if( norm(r,'inf') > sqrt(eps))
	disp('WARNING : precision may not be OK');
else
	disp('SCHUR SOLUTION CHECK2 OK');
end
norm(r,'inf')



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  NOW TRY REDUCED RHS FUNCTIONALITY 
%  (easier to use than previous
%   computations)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

id.JOB=3;
% Do forward solution step to obtain a reduced RHS
id.ICNTL(26)=1;
RHS=mat*ones(n,1);
id.RHS=RHS;
id = dmumps(id,mat);
% Solve the problem on the interface
id.REDRHS = id.SCHUR \ id.REDRHS;

% Do backward solution stage to expand the solution
id.ICNTL(26)=2;
id = dmumps(id,mat);
r = mat*id.SOL-RHS;
disp('*** check solution when REDRHS is used');
if( norm(r,'inf') > sqrt(eps))
	disp('WARNING : precision may not be OK');
else
	disp('SCHUR SOLUTION CHECK3 OK');
end
norm(r,'inf')