function F = factorize (A,strategy,burble)
%FACTORIZE an object-oriented method for solving linear systems
% and least-squares problems, and for representing operations with the
% inverse of a square matrix or the pseudo-inverse of a rectangular matrix.
%
% F = factorize(A) returns an object F that holds the factorization of A.
% x=F\b then solves a linear system or a least-squares problem.  S=inverse(F)
% or S=inverse(A) returns a factorized representation of the inverse of A so
% that inverse(A)*b is mathematically equivalent to pinv(A)*b, but the former
% does not actually compute the inverse or pseudo-inverse of A.
%
% Example
%
%   F = factorize(A) ;      % LU, QR, or Cholesky factorization of A
%   x = F\b ;               % solve A*x=b; same as x=A\b
%   S = inverse (F) ;       % S represents the factorization of inv(A)
%   x = S*b ;               % same as x = A\b.
%   E = A-B*inverse(D)*C    % efficiently computes the Schur complement
%   E = A-B*inv(D)*C        % bad method for computing the Schur complement
%   S = inverse(A) ; S(:,1) % compute just the first column of inv(A),
%                           % without computing inv(A)
%
%   F = factorize (A, strategy, burble) ;   % optional 2nd and 3rd inputs
%
% A string can be specified as a second input parameter to select the strategy
% used to factorize the matrix.  The first two are meta-strategies:
%
%   'default'   if rectangular
%                   use QR for sparse A or A' (whichever is tall and thin);
%                   use COD for full A
%               else
%                   if symmetric
%                       if positive real diagonal: try CHOL
%                       else (or if CHOL fails): try LDL
%                   end
%                   if not yet factorized: try LU (fails if rank-deficient)
%               end
%               if all else fails, or if QR or LU report that the matrix
%               is singular (or nearly so): use COD
%               This strategy mimics backslash, except that backslash never
%               uses COD.  Backslash also exploits other solvers, such as
%               specialized tridiagonal and banded solvers.
%
%   'symmetric' as 'default' above, but assumes A is symmetric without
%               checking, which is faster if you already know A is symmetric.
%               Uses tril(A) and assumes triu(A) is its transpose.  Results
%               will be incorrect if A is not symmetric.  If A is rectangular,
%               the 'default' strategy is used instead.
%
%   'unsymmetric'  as 'default', but assumes A is unsymmetric.
%
% The next "strategies" just select a single method, listed in decreasing order
% of generality and increasing order of speed and memory efficiency.  All of
% them except the SVD can exploit sparsity.
%
%   'svd'       use SVD.  Never fails ... unless it runs out of time or memory.
%                   Coerces a sparse matrix A to full.
%
%   'cod'       use COD.  Almost as accurate as SVD, and much faster.  Based
%                   on dense or sparse QR with rank estimation.  Handles
%                   rank-deficient matrices, as long as it correctly estimates
%                   the rank.  If the rank is ill-defined, use the SVD instead.
%                   Sparse COD requires the SPQR package to be installed
%                   (see http://www.suitesparse.com).
%
%   'qr'        use QR.  Reports a warning if A is singular.
%
%   'lu'        use LU.  Fails if A is rectangular; warning if A singular.
%
%   'ldl'       use LDL.  Fails if A is rank-deficient or not symmetric, or if
%                   A is sparse and complex.  Uses tril(A) and assumes triu(A)
%                   is the transpose of tril(A).
%
%   'chol'      use CHOL.  Fails if A is rank-deficient or not symmetric
%                   positive definite.  If A is sparse, it uses tril(A) and
%                   assumes triu(A) is the transpose of tril(A).  If A is dense,
%                   triu(A) is used instead.
%
% A third input, burble, can be provided to tell this function to print what
% methods it tries (if burble is nonzero).
%
% For a demo type "fdemo" in the Demo directory or see the Doc/ directory.
%
% See also inverse, slash, linsolve, spqr.

% Copyright 2011-2012, Timothy A. Davis, http://www.suitesparse.com

assert (ndims (A) == 2, 'Matrix must be 2D.') ;

if (nargin < 2 || isempty (strategy))
    strategy = 'default' ;
end
if (nargin < 3)
    burble = 0 ;
end

if (burble)
    fprintf ('\nfactorize: strategy %s, A has size %d-by-%d, ', ...
        strategy, size (A)) ;
    if (issparse (A))
        fprintf ('sparse with %d nonzeros.\n', nnz (A)) ;
    else
        fprintf ('full.\n') ;
    end
end

switch strategy

    case 'default'
        [F, me] = backslash_mimic (A, burble, 0) ;

    case 'symmetric'
        [F, me] = backslash_mimic (A, burble, 1) ;

    case 'unsymmetric'
        [F, me] = backslash_mimic (A, burble, 2) ;

    case 'svd'
        [F, me] = factorize_svd (A, burble) ;

    case 'cod'
        [F, me] = factorize_cod (A, burble) ;

    case 'qr'
        [F, me] = factorize_qr (A, burble, 0) ;

    case 'lu'
        % do not report a failure if the matrix is singular
        [F, me] = factorize_lu (A, burble, 0) ;

    case 'ldl'
        [F, me] = factorize_ldl (A, burble) ;

    case 'chol'
        [F, me] = factorize_chol (A, burble) ;

    otherwise
        error ('FACTORIZE:invalidStrategy', 'unrecognized strategy.') ;

end

if (~isobject (F))
    throw (me) ;
end

%-------------------------------------------------------------------------------

function [F, me] = backslash_mimic (A, burble, strategy)
%BACKSLASH_MIMIC automatically select a method to factorize A.
F = [ ] ;
me = [ ] ;
[m, n] = size (A) ;

% If the following condition is true, then the QR, QRT, or LU factorizations
% will report a failure if A is singular (or nearly so).  This allows COD
% or COD_SPARSE to then be used instead.  COD_SPARSE relies on the SPQR
% mexFunction in SuiteSparse, which might not be installed.  In this case,
% QR, QRT, and LU do not report failures for sparse matrices that are singular
% (or nearly so), since there is no COD_SPARSE to fall back on.
fail_if_singular = ~issparse (A) || (exist ('spqr') == 3) ;                 %#ok
try_cod = true ;

if (m ~= n)
    if (issparse (A))
        % Use QR for the sparse rectangular case (ignore 'strategy' argument).
        [F, me] = factorize_qr (A, burble, fail_if_singular) ;
    else
        % Use COD for the full rectangular case (ignore 'strategy' argument).
        % If this fails, there's no reason to retry the COD below.  If A has
        % full rank, then COD is the same as QR with column pivoting (with the
        % same cost in terms of run time and memory).  Backslash in MATLAB uses
        % QR with column pivoting alone, so this is just as fast as x=A\b in
        % the full-rank case, but gives a more reliable result in the rank-
        % deficient case.
        try_cod = false ;
        [F, me] = factorize_cod (A, burble) ;
    end
else
    % square case: Cholesky, LDL, or LU factorization of A
    switch strategy
        case 0
            is_symmetric = (nnz (A-A') == 0) ;
        case 1
            is_symmetric = true ;
        case 2
            is_symmetric = false ;
    end
    if (is_symmetric)
        % A is symmetric (or assumed to be so)
        d = diag (A) ;
        if (all (d > 0) && nnz (imag (d)) == 0)
            % try a Cholesky factorization
            [F, me] = factorize_chol (A, burble) ;
        end
        if (~isobject (F) && (~issparse (A) || isreal (A)))
            % try an LDL factorization.
            % complex sparse LDL does not yet exist in MATLAB
            [F, me] = factorize_ldl (A, burble) ;
        end
    end
    if (~isobject (F))
        % use LU if Cholesky and/or LDL failed, or were skipped.
        [F, me] = factorize_lu (A, burble, fail_if_singular) ;
    end
end
if (~isobject (F) && try_cod)
    % everything else failed, matrix is rank-deficient.  Use COD
    [F, me] = factorize_cod (A, burble) ;
end


%-------------------------------------------------------------------------------

function [F, me] = factorize_qr (A, burble, fail_if_singular)
% QR fails if the matrix is rank-deficient.
F = [ ] ;
me = [ ] ;
try
    [m, n] = size (A) ;
    if (m >= n)
        if (burble)
            fprintf ('factorize: try QR of A ... ') ;
        end
        if (issparse (A))
            F = factorization_qr_sparse (A, fail_if_singular) ;
        else
            F = factorization_qr_dense (A, fail_if_singular) ;
        end
    else
        if (burble)
            fprintf ('factorize: try QR of A'' ... ') ;
        end
        if (issparse (A))
            F = factorization_qrt_sparse (A, fail_if_singular) ;
        else
            F = factorization_qrt_dense (A, fail_if_singular) ;
        end
    end
    if (burble)
        fprintf ('OK.\n') ;
    end
catch me
    if (burble)
        fprintf ('failed.\nfactorize: %s\n', me.message) ;
    end
end


%-------------------------------------------------------------------------------

function [F, me] = factorize_chol (A, burble)
% LDL fails if the matrix is rectangular, rank-deficient, or not positive
% definite.  Only the lower triangular part of A is used.
F = [ ] ;
me = [ ] ;
try
    if (burble)
        fprintf ('factorize: try CHOL ... ') ;
    end
    if (issparse (A))
        F = factorization_chol_sparse (A) ;
    else
        F = factorization_chol_dense (A) ;
    end
    if (burble)
        fprintf ('OK.\n') ;
    end
catch me
    if (burble)
        fprintf ('failed.\nfactorize: %s\n', me.message) ;
    end
end


%-------------------------------------------------------------------------------

function [F, me] = factorize_ldl (A, burble)
% LDL fails if the matrix is rectangular or rank-deficient.
% As of MATLAB R2012a, ldl does not work for complex sparse matrices.
% Only the lower triangular part of A is used.
F = [ ] ;
me = [ ] ;
try
    if (burble)
        fprintf ('factorize: try LDL ... ') ;
    end
    if (issparse (A))
        F = factorization_ldl_sparse (A) ;
    else
        F = factorization_ldl_dense (A) ;
    end
    if (burble)
        fprintf ('OK.\n') ;
    end
catch me
    if (burble)
        fprintf ('failed.\nfactorize: %s\n', me.message) ;
    end
end


%-------------------------------------------------------------------------------

function [F, me] = factorize_lu (A, burble, fail_if_singular)
% LU fails if the matrix is rectangular or rank-deficient.
F = [ ] ;
me = [ ] ;
try
    if (burble)
        fprintf ('factorize: try LU ... ') ;
    end
    if (issparse (A))
        F = factorization_lu_sparse (A, fail_if_singular) ;
    else
        F = factorization_lu_dense (A, fail_if_singular) ;
    end
    if (burble)
        fprintf ('OK.\n') ;
    end
catch me
    if (burble)
        fprintf ('failed.\nfactorize: %s\n', me.message) ;
    end
end


%-------------------------------------------------------------------------------

function [F, me] = factorize_cod (A, burble)
% COD only fails when it runs out of memory.
F = [ ] ;
me = [ ] ;
try
    if (burble)
        fprintf ('factorize: try COD ... ') ;
    end
    if (issparse (A))
        F = factorization_cod_sparse (A) ;
    else
        F = factorization_cod_dense (A) ;
    end
    if (burble)
        fprintf ('OK.\n') ;
    end
catch me
    if (burble)
        fprintf ('failed.\nfactorize: %s\n', me.message) ;
    end
end

%-------------------------------------------------------------------------------

function [F, me] = factorize_svd (A, burble)
% SVD only fails when it runs out of memory.
F = [ ] ;
me = [ ] ;
try
    if (burble)
        fprintf ('factorize: try SVD ... ') ;
    end
    F = factorization_svd (A) ;
    if (burble)
        fprintf ('OK.\n') ;
    end
catch me
    if (burble)
        fprintf ('failed.\nfactorize: %s\n', me.message) ;
    end
end