## File: schur_example.m

package info (click to toggle)
mumps 5.1.2-5
 `1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192` ``````%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') ``````