pmcont

 parallel multi continuation routine 
 for -c lap(u)+au-f-b grad u=0 
 Do p.nstep multisteps of (initial) size i*p.ds, i=1:p.mst,

0001 function p=pmcont(p)
0002 % parallel multi continuation routine
0003 % for -c lap(u)+au-f-b grad u=0
0004 % Do p.nstep multisteps of (initial) size i*p.ds, i=1:p.mst,
0005 cstop=0; is=0; p0=p; m=p.mst; scount=p.count; % save counter for timing
0006 if(p.pnamesw==1) % check if dir exists
0007     dum1=p.count; dum2=p.branch; dum3=p.bcount; 
0008     [p,ok]=setfn(p,inputname(1)); 
0009     p.count=dum1; p.branch=dum2; p.bcount=dum3; 
0010     if ok~=1; return; end
0011 end
0012 p.ntime=0; p.tstime=0; % timers for pmnewtonloop and multisteps!
0013 totime=tic; % total time timer
0014 if (matlabpool('size') == 0) matlabpool open; end %open cores
0015 p.headfu(p); 
0016 for k=1:p.nsteps   
0017   stime=tic; 
0018   if(p.err>p.errbound)                            
0019     if(p.errchecksw==2); p=meshadac(p,'eb',p.errbound); end 
0020     if(p.errchecksw>2); p=meshref(p,'eb',p.errbound); end          
0021   end 
0022   if ((length(p.tau)~=p.neq*p.np+1)||(p.restart>0)) % initial step
0023       [p,iok]=inistep(p); if(iok==0); return; end 
0024       plotsol(p,1,1,p.pstyle); is=is+2;
0025   end
0026   ineg0=p.ineg; % use last ineg for bifdetec consistency check
0027   disp(['multi continuation step ' num2str(k)] );
0028   p.pmimax=p0.pmimax;   % set p.pmimax to orignal one
0029   foso=0; %foso=found solution
0030   while foso==0     
0031      ntime=tic;  % timer for mst newton-loops!
0032      [uold,lamold,tauold,inegp,resp,dsp,f]=pmnewtonloop(p,p0); % main call!
0033      p.ntime=p.ntime+toc(ntime); % time for p.mst newtonloops
0034      for(i=1:p.mst) if (resp(i)<p.tol) foso=foso+1; end; end;
0035      if (foso==0 && abs(p.ds)<p.dsmin*(p.mst+1) && p.imax>p.pmimax) % sscontrol
0036          p.pmimax=p.pmimax+1; disp(['increasing pmimax to ' num2str(p.pmimax)]); end
0037      if (foso==0 && abs(p.ds)>=p.dsmin*(p.mst+1)) p.ds=p.ds/(p.mst+1); 
0038          disp(['decreasing ds to ' num2str(p.ds)]); end
0039      if (foso==0 && abs(p.ds)<p.dsmin*(p.mst+1) && p.imax==p.pmimax) 
0040          disp('no solution found'); foso=-1; end
0041    end %while foso
0042    p.imax=p0.imax; % "postprocessing", check bif, plot, calc. tangent
0043    msucc=0; 
0044   for i=1:m % check bif, plot, calc. tangent
0045     if resp(i)<p.tol        
0046        msucc=msucc+1; p.restart=0; % step accepted!
0047        Gu=f(i).Gu; Glam=f(i).Glam; p.iter=f(i).iter; ds=dsp(i); p.res=resp(i);
0048        lam1=f(i).lam; u11=f(i).u; p.meth=f(i).meth; p.ineg=inegp(i); 
0049        p.err=f(i).err; % HU
0050        amat=[[Gu Glam]; [p.xi*p.tau(1:p.neq*p.np)' (1-p.xi)*p.tau(p.neq*p.np+1)]];
0051        tau=p.blss(amat,[zeros(p.neq*p.np,1);1],p,lam1); % tau at new point
0052        tau=tau/xinorm(tau,p.xi);
0053        if(p.bifchecksw>0 ) % check for bif via |det(A)|
0054          p.u=uold; p.lam=lamold; p.tau=tauold; 
0055          if(i>1); ineg0=inegp(i-1); end % HU
0056          p=bifdetec(p,u11,lam1,tau,Gu,Glam,ineg0); % HU
0057        end
0058        p.u=u11; p.lam=lam1; p.tau=tau; p.lamd=p.tau(p.neq*p.np+1);      
0059        if(mod(is,p.brmod)==0) % put step on branches
0060          brout=p.outfu(p,p.u,p.lam); % userfu to append to bif-branches
0061          p.branch=[p.branch [p.count; p.ineg; p.lam; p.err; brout]]; % put on branch, HU
0062          figure(p.brfig); hold on; % plot point
0063          if(p.ineg<=0) plot(lam1,brout(p.bpcmp),'*'); 
0064          else plot(lam1,brout(p.bpcmp),'+'); end 
0065        end 
0066        if(p.count>0 && mod(p.count,p.smod)==0) % save to file
0067          fname=[p.pname,sprintf('%i',p.count),'.mat']; save(fname,'p'); 
0068        end
0069        if(mod(p.count,p.pmod)==0); plotsol(p,1,1,p.pstyle); axis tight;  end % plot sol
0070        uold=u11; lamold=lam1; tauold=tau; %for bifdetec
0071        cstop=p.ufu(p,brout,ds);
0072        p.count=p.count+1; is=is+1; 
0073        if cstop==1; break; end
0074     end %res<p.tol
0075     p.imax=p0.imax; %because is used in Newton loop and bifdetect
0076   end % i=1:m in postprocessing
0077   if(p.smod~=0); fname=[p.pname,sprintf('%i',p.count-1),'.mat']; save(fname,'p'); end % save last point
0078   if cstop==1 break; end
0079   if (foso==p.mst && abs(p.ds)<=p.dsmax/p.dsincfac) p.ds=p.dsincfac*p.ds; 
0080      disp(['increasing ds to ' num2str(p.ds)]); end
0081   if foso==-1 break; end  
0082   if(mod(k,p.amod)==0); fprintf('adapting mesh\n'); p=meshadac(p); end 
0083   cstime=toc(stime); p.tstime=p.tstime+cstime;
0084   if p.timesw>0; fprintf('time for %2i successful steps: %g\n',msucc,cstime);end
0085 end %nstep loop
0086 if (matlabpool('size') ~= 0); matlabpool close force; end  %close cores
0087 p.totime=toc(totime); 
0088 if p.timesw>0 
0089     fprintf('Total time=%g, av.time per succ.step=%g, av.newton-time=%g\n', ... 
0090         p.totime, p.totime/(p.count-scount),p.ntime/(p.count-scount)); 
0091 end