function ks = ofact(k,ind,varargin);

%OFACT Ofact matrix class constructor
%
%	ks = ofact(k)  builds a ofact object from sparse matrix k
%           ks\b can then be used to solve problem k q = b in an optimized
%           fashion, including repeated solves.  
%           for single solves ofact(k,b) replaces k\b
%       ofact                      lists availlable methods for this system
%       ofact('method')            returns the current method
%       ofact('MethodName') sets current method 
%       ofact('silent') attempts to configure current method to be silent
%
%       OBSOLETE methods supported still supported by ofact are 
%
%       ks = ofact(k,ind) builds a ofact matrix from a block of k. The
%	   indices are supposed to be ordered and are a subset of all indices
%	ks = ofact(k.data,k.ind) builds a ofact matrix from its components
%       ks = ofact(k,'lu') stores sparse 'lu' factors for static computations
%       ks = ofact(k,'ch') stores sparse Cholesky factors for static
%	   computations
%	ks = ofact(k,'de') uses the 'lu' factorization for small and complex
%	   matrices and the ofact form otherwise
%
%       SDTools maintains links to compiled solvers for ofact at
%         http://www.sdtools.com/faq/FE_ofact.html
%	
%	See also help fe_reduc, fe_eig
%		 doc  ofact

%       Etienne Balmes
%       Copyright (c) 2001-2005 by INRIA and SDTools, All Rights Reserved.
%       $Revision: 1.17 $  $Date: 2006/05/15 14:29:21 $
%       Use under OpenFEM trademark.html license and LGPL.txt library license


persistent methods

if isempty(methods) % initialize method selection

 methods=struct('name','Undefined ofact method');
 methods=ofact('methodsde');

end

if nargin==2&ischar(ind) a=ofact(['method' ind]); end

if nargin==0  

  if nargout==0 ofact('methodlist'); 
  else ks=methods.name;
  end
  return
  %ks.ty = 0; ks.ind =[]; ks.data=[]; ks.dinv=[]; ks.l=[]; ks.u=[];

elseif ischar(k)
% ------------------------------------------------------------------------
% ------------------------------------------------------------------------
% ------------------------------------------------------------------------
% commands - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 if comstr(k,'close')|comstr(k,'clear')

  met=[];
  if nargin<2 met=methods;ind=[];
  else 
   try; met=ind.method;end
   if isempty(met) met=methods; end
  end
  
  if     comstr(methods.name,'psldlt') sky_fig('close');
  elseif ~isstruct(met)
  elseif ~isempty(met.Clear)
    if strcmp(class(ind),'ofact') ks=ind;
    else 
     if isempty(ind) ind=-1; end
     ks=empty_fact; ks.ty(2)=ind(1);
    end
    eval(met.Clear)
    clear ks
  end

 elseif comstr(comstr(k,-27),'silent') % - - - - - - - - - - - - - - - - - - -
 % Attempt to make current method silent
  try; methods=feval(sprintf('%s_utils',methods.name), ...
    'silent',empty_fact,methods);
  catch; 
   sdtw('Silent failed'); 
  end

 elseif ~isempty(strfind(comstr(k,-27),'fact')) % - - - - - - - - - - - - - - - -

   ks=empty_fact; ks.method=methods;
   ks=feval([methods.name '_utils'],k,ind,ks,varargin{:});

 elseif comstr(comstr(k,-27),'symrenumber') % - - - - - - - - - - - - - - - -

  if isfield(methods,'SymRenumber') ks=methods.SymRenumber;
  else ks='symmmd'; end

 else %if strncmpi(k,'method',6) % - - - - - - - - - - - - - - - - - - -

   k0=empty_fact; 
   i1=strfind(k,'complex'); 
   if ~isempty(i1) NeedComplex=1; k(i1+[0:6])='';else NeedComplex=0; end


   ks=[];
   ks=spfmex_utils('method',k0);ks.param=[];
 
   ks(end+1).name='lu';
   ks(end).SymRenumber='symmmd';
   ks(end).header='MATLAB sparse LU solver';
   ks(end).FactBuild='ks=lu_fact(k,ks);';
   ks(end).Available=1;
   ks(end).HandlesComplex=1;

   ks(end+1).name='chol';
   ks(end).SymRenumber='symmmd';
   ks(end).header='MATLAB sparse Cholesky solver';
   ks(end).FactBuild='ks=chol_fact(k,ks);';
   ks(end).Available=1;
   ks(end).HandlesComplex=1;

   try; eval('ks(end+1)=mtaucs_utils(''method'',k0);');end
   try; eval('ks(end+1)=pardiso_utils(''method'',k0);');end

   ks(end+1).name='sp_util';
   ks(end).SymRenumber='symrcm';
   ks(end).header='SDT skyline solver';
   ks(end).FactBuild='ks=spu_fact(k,ks);';
   ks(end).Available=(sp_util>5);
   ks(end).HandlesComplex=0;

   ks(end+1).name='psldlt';
   ks(end).header='SGI sparse solver';
   ks(end).SymRenumber='symmmd';
   ks(end).FactBuild='error(''Not availlable'')';
   ks(end).Available=exist('psldlt')==3;
   ks(end).HandlesComplex=0;

   ks(end+1).name = 'umfpack';
   ks(end).SymRenumber = '';
   ks(end).header = 'UMFPACK solver';
   ks(end).FactBuild = 'ks = umf_fact( k, ks );';
   ks(end).Available = exist('umfpack') == 3;
   ks(end).HandlesComplex = 1;
   ks(end).Solve ='q = umfpack( k.data, ''\'', full(b) );';
   ks(end).Clear ='';

   if strncmpi(k,'method',6) [k,Cam]=comstr(k,7);else [k,Cam]=comstr(k,1); end
   % Eliminate real solvers if complex matrix
   if nargin==2&isnumeric(ind) & ~isreal(ind)
     ks=ks(find([ks.HandlesComplex]));
   end

   i1=strmatch(Cam,{ks.name});

   if comstr(Cam,'list') % - - - - - - - - - - - - - - - - - - - - - - - -
    fprintf('\nAvailable factorization methods for OFACT object\n\n');
    ind=1:length(ks); for j1=1:length(ind)
      if ks(j1).Available
       if strcmp(methods.name,ks(j1).name)
        fprintf('->%8s : %s\n',ks(j1).name,ks(j1).header); 
       else
        fprintf('%10s : %s\n',ks(j1).name,ks(j1).header); 
       end
      else
        fprintf('%10s : %s (NOT AVAILABLE ON THIS MACHINE)\n', ...
         ['*' ks(j1).name],ks(j1).header); 
      end
    end

    fprintf('\nuse ofact(''method MethodName'') to select.\n');
    
    clear ks; return;
   % - - - - - - - - - - - - - - - - - - - - - - - -
   elseif length(i1)==1&ks(i1).Available 

    ks=ks(i1(1));
    if nargin>=2 & ischar(ind)
     args={ind varargin{:}};carg=1;
     while carg<=length(args)
      if isfield(ks,args{carg}) % set a field
       ks=setfield(ks,args{carg},args{carg+1});carg=carg+2;
      else
       eval(sprintf('ks=%s_utils(args{carg},k0);',ks.name));carg=carg+1;
      end
     end
    end

   % - - - - - - - - - - - - - - - - - - - - - - - -
   else
        r1=methods;
        if ~isfield(r1,'Available')|~r1.Available
          for j1=1:length(ks)
            if ks(j1).Available methods=ks(j1);methods=ks(j1);break;end
          end
        end
        ks=methods;
   end
   if nargout==0 
     methods=ks;
     disp(ks);
     if strcmp(ks.name,'pardiso') 
        warning('Pardiso is set for symmetric matrices by default')
     end
     clear ks;
   end

 end

 return;
end
% ------------------------------------------------------------------------
% ------------------------------------------------------------------------
% ------------------------------------------------------------------------
% actual object building 

if ~any(exist(methods.name)==[2 3 5 6])
  a=ofact('methodsde');
end


if isa(k,'ofact') ks = k;
  if nargin==2 error('k is already a ofact object'); end

% give the sp_util fields in the input
elseif nargin==2 & size(k,2)==1 & ind(end)==size(k,1)

 if ~isreal(k)
  error('Complex ofact matrices are not handled by sp_util');
 end
 ks.ty = 0; ks.ind=ind; ks.data=k;
 if ks.ind(end)~=length(ks.data)|ks.data(end)~=0
  warning('.ind and .data fields don''t seem to be consistent');
 end
 ks.dinv=[]; ks.l=[]; ks.u=[]; ks.method=ofact('method sp_util');
 ks=class(ks,'ofact');
 
elseif size(k,1)~=size(k,2) & nargin==1

  error('Ofact matrices are always square');

elseif ~isa(k,'double')&~isa(k,'sparse')

  error('Ofact objects are created from sparse or full matrices')

elseif nargin==1 % standard call 

  ks=struct('ty',[5 0 0], ... 
      'ind',[],'data',[],'dinv',[],'l',[],'u',[]);
  ks.method=methods;
  if ~isreal(k)&~ks.method.HandlesComplex
   ks.method=ofact(['method de'],k);
  end
  ks = class(ks,'ofact');
  eval(methods.FactBuild);

elseif nargin==2 % standard call with method selection

  if ischar(ind)& strcmp(ind,'diag')

   r1=full(diag(k)); if nnz(r1)==length(r1); ind=[];else ind=find(r1);end
   ks=struct('ty',[6 size(k)], ... 
      'ind',ind,'data',[],'dinv',1./r1(find(r1)),'l',[],'u',[], ...
      'method','diag');
   if isempty(ind) ks.ty(1)=6.1;end
   ks = class(ks,'ofact');

  elseif ischar(ind)

   ks=struct('ty',[5 0 0], ... 
      'ind',[],'data',[],'dinv',[],'l',[],'u',[]);
   ks.method=ofact(['method ' ind]);
   if ~isreal(k)&~ks.method.HandlesComplex
    ks.method=ofact(['method de'],k);
   end
   ks = class(ks,'ofact');
   eval(ks.method.FactBuild);

  elseif isnumeric(ind) & min(size(k))==1  % OBSOLETE ofact call
    ks=struct('ty',[0 0 0], ... 
      'ind',ind,'data',k,'dinv',[],'l',[],'u',[]);
    ks.method=ofact('method sp_util');
    if ~isfield(ks,'dinv'); ks.dinv=[]; end 
    ks = class(ks,'ofact'); 
  % Single solve with clear object 
  elseif isnumeric(ind)&size(ind,1)==size(k,2)

    kd=ofact(k); ks=kd\ind; ofact('clear',kd);

  elseif isnumeric(ind) % OBSOLETE ofact call
    k=k(ind,ind); sp_util('spind',k,ind); 
    ks = sp_util('sp2sky',k,rind,length(ind)); 
    if ~isfield(ks,'dinv'); ks.dinv=[]; end 
    ks = class(ks,'ofact'); 
  end


end 
  

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
function out=sky_fig(Cam)

  gf=findall(0,'tag','ofact');
  if isempty(gf)
   gf=figure('name','ofact','visible','off','tag','ofact', ...
     'integerhandle','off','HandleVisibility','off', ...
     'CloseRequestFcn','ofact(''close'')');
  end

  uf=get(gf,'userdata');
  if ~isfield(uf,'LDLT') uf.LDLT=[]; end

  if comstr(Cam,'new')
    uf.LDLT(end+1)=length(uf.LDLT)+1;
    out=uf.LDLT(end);
  elseif comstr(Cam,'close')
    for j1=1:length(uf.LDLT)
      psldlt('close',j1);
    end
  end

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
function ks=lu_fact(k,ki)

   ks.ty = 2; ks.ind=[];  ks.data=[]; ks.dinv=[];
   if issparse(k) [ks.l,ks.u] = lu(k,0); else  [ks.l,ks.u] = lu(k); end
   ks.method=ki.method;
   ks = class(ks,'ofact'); 

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
function ks=chol_fact(k,ki)

   ks.ty = 3; ks.ind=[]; ks.data=[]; ks.dinv=[]; ks.l=[];
   ks.u = chol(k); ks.l = ks.u';
   ks.method=ki.method;
   ks = class(ks,'ofact'); 

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
function ks=spu_fact(k,ki)


 ks = sp_util('sp2sky',k);
 ks.method=ofact('sp_util');
 ks = class(ks,'ofact'); 

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
function ks=empty_fact

     ks=struct('ty',[0 0 0], ... 
       'ind',0,'data',[],'dinv',[],'l',[],'u',[]);
     ks.method='';
     ks = class(ks,'ofact'); 

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
function ks=psldlt_fact(k,ki)

    i1=sky_fig('new');
    ks=struct('ty',[5 i1 size(k,1)], ...
      'ind','psldlt','data',[],'dinv',[],'l',[],'u',[]);
    psldlt('factor',k,i1);
    ks.method=ki.method;
    ks = class(ks,'ofact'); 

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
function ks=umf_fact(k,ki)

ks=struct('ty',[10, 0, size( k, 1 )], ...
	  'ind',[],'data',k,'dinv',[],'l',[],'u',[]);
ks.method=ki.method;
ks = class(ks,'ofact'); 
