function ks=spfmex_utils(varargin);

% Gateway function for inclusion of the SPFMEX solver into the 
% OpenFEM ofact object
%
% Standard Factor building method
%
%   ofact('spfmex','FactBuild', ...
%    'ks=spfmex_utils(''fact'',k,ks,[1 0 size(k,1)/32  0.01*size(k,1)]);');
%
% Call to the external SDT_Server. 
%  ofact('spfmex','FactBuild', ...
%   'ks=spfmex_utils(''servfact'',k,ks,[1 0 size(k,1)/32  0.01*size(k,1)]);');
%
%  To close the server down : pnet('closeall')
%
% You can set host and port using
%  setpref('SDT','ServerHost','localhost');
%  setpref('SDT','ServerPort',888);
%
%       $Revision: 1.14 $  $Date: 2006/05/18 15:39:45 $

Cam=varargin{1};carg=2;
if comstr(Cam,'method')

   ks.name='spfmex';
   ks.header='SDT sparse LDLt solver';
   ks.SymRenumber='';
   ks.FactBuild='ks=spfmex_utils(''fact'',k,ks,[1 0 size(k,1)/32  0.01*size(k,1)]);';
%   ks.Solve='q=spfmex(''solve'',k.ty(2),full(b));';
   ks.Solve='q=spfmex_utils(''solve'',k,full(b));';
   ks.Clear='spfmex_utils(''clear'',ks);';
   ks.Available=exist('spfmex')==3;
   ks.HandlesComplex=1;

% ---------------------------------------
elseif comstr(Cam,'serv');[CAM,Cam]=comstr(Cam,5);

% 'ServFact server factorization' ---------------------------
if comstr(Cam,'fact')

 k=varargin{carg};carg=carg+1;
 ki=varargin{carg};carg=carg+1;

 method=ki.method;
 method.host=getpref('SDT','ServerHost','localhost');
 method.port=getpref('SDT','ServerPort',8888);
 
 gf=findall(0,'type','figure','tag','rwarn');
 if isempty(gf)
 gf=figure('Numbertitle','off', ...
     'IntegerHandle','off','units','pixels','tag','rwarn',...
     'CloseRequestFcn','cingui(''rwarn_HIDE'');','visible','off');
 end
 go=findall(gf,'tag','sdt_serv');
 if isempty(go)
   go=uicontrol('Style','list','String','','tag','sdt_serv', ...
        'value',1,'parent',gf,  ...
        'callback','','visible','off','HandleVisibility','off');
 end
 uo=get(go,'userdata');
 if ~isfield(uo,'con')|eval('pnet(uo.con,''status'')<=0','1')
  eval('con=pnet(''tcpconnect'',method.host,method.port);');
 else con=uo.con;
 end

 if con==-1, % If server not running try starting it 
  if isunix 
   st=which(sprintf('sdt_serv_%s.exe',strrep(mexext,'mex','')));
   st=sprintf('unix(''xterm -T SDT_Server -e %s &'')',st);
  else   st=which('sdt_serv.exe');
         st=sprintf('dos(''%s &'',''-echo'');',st);
  end
  st=getpref('SDT','ServerStart',st);
  eval(st);
  pause(1);  eval('con=pnet(''tcpconnect'',method.host,method.port);');
 end
 if con==-1 error 'Bad url or server down.....'; end

 uo.con=con; set(go,'userdata',uo);
 eval('pnet(con,''setwritetimeout'',10);');
 eval('pnet(con,''setreadtimeout'',10);');

 eval('pnet(con,''printf'',''SPOOLESFACT\n'');')

 k=k+spalloc(size(k,1),size(k,1),0);
 eval('pnet(con,''write'',int32([size(k) nnz(k) 0]),''intel'');')
 eval('pnet(con,''sparsewrite'',k,''intel'');')
 fprintf('Matrix transmitted');

 try;
  st='';
  while isempty(st)  
   pause(1); % Try to do some SMART GUESS of how long to wait
   eval('st=pnet(con,''readline'',''noblock'');')
   stat=pnet(con,'status'); if stat<=0 error('The server has diconnected');end
  end
 end
  if ~strcmp(st,'DONE') error('Problem during factorization');end
  i1=double(pnet(con,'Read',1,'uint32')); % factor number

 pnet(con,'printf','EMPTY\n');
 method.con=con;
 method.Solve='q=spfmex_utils(''servsolve'',k,full(b));';
 method.Clear='spfmex_utils(''servclear'',ks);';
 ki.method=method;
 ks=ki; ks.ty(2:3)=[i1 size(k,1)]; ks.ind='sdt_serv';

% 'ServSolve server solve' ---------------------------
elseif comstr(Cam,'solve') % -------------------------------------

 ki=varargin{carg};carg=carg+1;
 b=varargin{carg};carg=carg+1;

 con=ki.method.con;

  pnet(con,'printf','SPOOLESSOLVE\n'); 
 
  cF=0;
  pnet(con,'write',int32([size(b) ki.ty(2) 0]),'intel');
  eval('pnet(con,''write'',b,''intel'');')

  dataclass='';
  while isempty(dataclass)  
   pause(1); % Try to do some SMART GUESS of how long to wait
   dataclass=pnet(con,'readline',1024,'noblock');
   stat=pnet(con,'status'); if stat<=0 error('The server has diconnected');end
  end

  datadims=double(pnet(con,'Read',1,'uint32','intel'));
  datasize=double(pnet(con,'Read',datadims,'uint32','intel'));
  fprintf('Receiving %i %i %s\n',datasize,dataclass);
  ks=pnet(con,'Read',datasize,'double','intel');

  if size(ks,1)~=size(b,1) error('Improper return');end
  pnet(con,'printf','EMPTY\n');


% 'ServSolve server Clear' ---------------------------
elseif comstr(Cam,'clear') % -------------------------------------

  ki=varargin{carg};carg=carg+1;
  try;
   con=ki.method.con;
   pnet(con,'printf','SPOOLESCLEAR\n'); 
   pnet(con,'write',int32([ki.ty(2)]),'intel');
   pnet(con,'printf','EMPTY\n');
  catch
   pnet('closeall');
  end
% 'ServOn server On' ---------------------------
elseif comstr(Cam,'on')

ks=ofact('spfmex','FactBuild', ...
 'ks=spfmex_utils(''servfact'',k,ks,[1 1 size(k,1)/32  0.01*size(k,1)]);',...
 'Solve','q=spfmex_utils(''servsolve'',k,full(b));', ...
 'Clear','spfmex_utils(''servclear'',ks)');

% 'ServOn server Off' ---------------------------
elseif comstr(Cam,'off')

ks=ofact('spfmex','FactBuild', ...
 'ks=spfmex_utils(''fact'',k,ks,[1 0 size(k,1)/32  0.01*size(k,1)]);', ...
 'Solve','q=spfmex_utils(''solve'',k,full(b));', ...
 'Clear','spfmex_utils(''clear'',ks);');

else sdtw('''Serv%s'' unknown',Cam);
end

elseif comstr(Cam,'fact') % -----------------------------------------

 k=varargin{carg};carg=carg+1;
 ki=varargin{carg};carg=carg+1;
 if carg>nargin opt=[1 0 size(k,1)/32  0.01*size(k,1)];
 else opt=varargin{carg};carg=carg+1;
 end

 if ~issparse(k) k = sparse(k); end
 t=cputime;
 i1=find(diag(k)==0);
 if ~isempty(i1) 
     if length(i1)<10
      error(sprintf('Zero terms on diagonal :%i\n',i1)); 
     else error(sprintf('%i zero terms on diagonal \n',length(i1))); end   
     return
 end
 % opt=[matlabmsglvl msglvl maxdomainsize maxzeros]

 if spfmex==5.0001 i1=spfmex('fact',k); 
 elseif spfmex==5.0003 [i1,ind]=spfmex('fact',k,opt); 
 else i1=spfmex('fact',k,opt); end
 if opt(1) fprintf(' %.2f s',cputime-t); end
 ks=ki; ks.ty(2:3)=[i1 size(k,1)]; ks.ind='spfmex';
 if spfmex==5.0003 ks.dinv=ind; end

elseif comstr(Cam,'solve')

 k=varargin{carg};carg=carg+1;
 b=varargin{carg};carg=carg+1;
 
 if size(b,1)~=size(k,2) ; error('RHS has a bad size'); end

 if ~isreal(b) | ~spfmex('isreal',k.ty(2))
  if isempty(k.dinv) % complex numbering outside
    ks=spfmex('solve',k.ty(2),complex(full(b)));
  else % complex numbering here
    if ~isreal(b)
      ks(k.dinv,:)=spfmex('solve',k.ty(2),full(b(k.dinv,:)));
    else
      ks(k.dinv,:)=spfmex('solve',k.ty(2),complex(full(b(k.dinv,:))));
    end
  end
 elseif isempty(k.dinv)  % real numbered outside
  ks=spfmex('solve',k.ty(2),full(b));
 else % real numbering here
  if size(k.dinv,1)~=k.ty(3) error('Bad size of RHS'); end
  ks=zeros(size(b,1),size(b,2));
  %ks(k.dinv,:)=spfmex('solve',k.ty(2),full(b(k.dinv,:)));
  b=full(b(k.dinv,:)); spfmex('solve',k.ty(2),b,ks,int32(k.dinv-1));
 end
 

elseif comstr(Cam,'clear')

 ks=varargin{carg};carg=carg+1;
 if isempty(ks)  spfmex('clear',-1);
 elseif isa(ks,'double') 
    try; spfmex('clear',ks);catch; disp(lasterr);end
 else
    try;  spfmex('clear',ks.ty(2)); 
    catch;disp(lasterr);
    end
 end

elseif comstr(Cam,'silent') % - - - - - - - - - - - - - - - - - -

 k=varargin{carg};carg=carg+1;
 ki=varargin{carg};carg=carg+1;
 ki.FactBuild='ks=spfmex_utils(''fact'',k,ks,[0 0 size(k,1)/32  0.01*size(k,1)]);';
 ks=ki;

end
