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

%SKYLINE Skyline matrix class constructor
%
%	ks = skyline(k)  builds a skyline object from sparse matrix k
%           ks\b can then be used to solve problem k q = b in an optimized
%           fashion, including repeated solves.  
%       skyline                      lists availlable methods for this system
%       skyline('method')            returns the current method
%       skyline('method MethodName') sets current method 
%
%       OBSOLETE methods supported still supported by skyline are 
%
%       ks = skyline(k,ind) builds a skyline matrix from a block of k. The
%	   indices are supposed to be ordered and are a subset of all indices
%	ks = skyline(k.data,k.ind) builds a skyline matrix from its components
%       ks = skyline(k,'lu') stores sparse 'lu' factors for static computations
%       ks = skyline(k,'ch') stores sparse Cholesky factors for static
%	   computations
%	ks = skyline(k,'de') uses the 'lu' factorization for small and complex
%	   matrices and the skyline form otherwise
%	
%	See also help fe_reduc, fe_eig
%		 doc  skyline


%       Etienne Balmes  07/28/92, 16/03/02
%       Copyright (c) 2001-2002 by INRIA and SDTools

%       Use under OpenFEM trademark.html license and LGPL.txt library license
%       All Rights Reserved.

persistent methods

if isempty(methods) % initialize method selection

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

end

if nargin==2&isstr(ind) a=skyline(['method' ind]); end

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

elseif isstr(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 ~isempty(met.Clear)
    if strcmp(class(ind),'skyline') 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(lower(k),'symrenumber') % - - - - - - - - - - - - - - - -

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

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

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


   ks=[];
   ks=spfmex_utils('method',k0);
  

   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;

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

   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;


   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 SKYLINE 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 skyline(''method MethodName'') to select.\n');
    
    clear ks; return;
   % - - - - - - - - - - - - - - - - - - - - - - - -
   elseif length(i1)==1&ks(i1).Available 

    ks=ks(i1(1));
    if nargin>2 & isstr(ind)
     args={ind varargin{:}};carg=1;
     while carg<length(args)&isfield(ks,args{carg}) % set a field
       ks=setfield(ks,args{carg},args{carg+1});carg=carg+2;
      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);clear ks; end

 end

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

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


if isa(k,'skyline') ks = k;
  if nargin==2 error('k is already a skyline 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 skyline 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=skyline('method sp_util');
 ks=class(ks,'skyline');
 
elseif size(k,1)~=size(k,2) & nargin==1

  error('Skyline matrices are always square');

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

  error('Skyline 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=skyline(['method de'],k);
  end
  ks = class(ks,'skyline');
  eval(methods.FactBuild);

elseif nargin==2 % standard call with method selection

  if isstr(ind)

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

  elseif isnumeric(ind) & min(size(k))==1  % OBSOLETE skyline call
    ks=struct('ty',[0 0 0], ... 
      'ind',ind,'data',k,'dinv',[],'l',[],'u',[]);
    ks.method=skyline('method sp_util');
    if ~isfield(ks,'dinv'); ks.dinv=[]; end 
    ks = class(ks,'skyline'); 
  elseif isnumeric(ind) % OBSOLETE skyline 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,'skyline'); 
  end


end 
  

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

  gf=findall(0,'tag','skyline');
  if isempty(gf)
   gf=figure('name','skyline','visible','off','tag','skyline', ...
     'integerhandle','off','HandleVisibility','off', ...
     'CloseRequestFcn','skyline(''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,'skyline'); 

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
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,'skyline'); 

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


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

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

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

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
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,'skyline'); 
