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')
|
