1//A simple demo for the MUMPS interface, with the return of the schur complement
2//to run it, You just have to execute  the instruction within Scilab
3//		exec sparse_example.sce;
4
5
6//*********************** MATRIX INITIALISATION ***********************//
7	n=10;
8	mat=sprand(n,n,.5)+speye(n,n);
9        size_schur=3;
10
11// Right Hand side setting
12	RHS = ones(n,1);
13
14
15//****************** Initialisation of the Scilab MUMPS structure ******************//
16timer();
17[id]=initmumps();
18
19//Here Job=-1, the next call will only initialise the C and Fortran structure
20[id]=dmumps(id);
21
22id.RHS=RHS;
23id.VAR_SCHUR = [n-size_schur+1:n];
24
25//******************** CALL TO MUMPS FOR RESOLUTION ON INTERNAL PROBLEM ************//
26job=6;
27id.JOB=job;
28
29[id]=dmumps(id,mat);
30
31// verification of the solution
32solution=id.SOL;
33norm1=norm(mat(1:n-size_schur,1:n-size_schur)*solution(1:n-size_schur) - ones(n-size_schur,1),'inf');
34if norm1> 10^(-9) then
35	write(%io(2),'WARNING: solution on internal problem may not be OK');
36else
37	write(%io(2),'SOLUTION on internal problem ok');
38end
39
40
41//******************* TRY REDUCED RHS FUNCTIONALITY **************//
42id.JOB=3;
43id.ICNTL(26)=1;
44
45// Forward
46[id]=dmumps(id,mat);
47
48// Solve the problem on the Schur complement
49id.REDRHS=id.SCHUR \ id.REDRHS;
50
51// and reinject it to MUMPS
52id.ICNTL(26)=2;
53[id]=dmumps(id,mat);
54solution=id.SOL;
55norm1=norm(mat*solution-RHS,'inf')
56if norm1> 10^(-9) then
57	write(%io(2),'WARNING: solution on complete problem may not be OK');
58else
59	write(%io(2),'SOLUTION on complete problem ok');
60end
61
62
63
64//****************** DESTRUCTION OF THE MUMPS INSTANCE ******************//
65job=-2;
66id.JOB=job;
67[id]=dmumps(id);
68t=timer()
69