function [W,L,U,pcol,P] = fzmult(A,V,doTranspose,L,U,pcol,P)
%

%FZMULT Multiplication with fundamental nullspace basis.
%   W = FZMULT(A,V) computes the product W of matrix Z with matrix V, that
%   is, W = Z*V, where Z is a fundamental basis for the nullspace of matrix
%   A. A must be a sparse m-by-n matrix where m < n, rank(A) = m, and
%   rank(A(1:m,1:m)) = m.  V must be p-by-q, where p = n-m. If V is sparse
%   W is sparse, else W is full.
%
%   W = FZMULT(A,V,true) computes the product of the transpose of
%   the fundamental basis times V, that is, W = Z'*V. V must be p-by-q,
%   where q = n-m. 
%
%   [W,L,U,pcol,P]  = FZMULT(A,V) returns the sparse LU-factorization of
%   matrix A(1:m,1:m), that is, A1 = A(1:m,1:m) and P*A1(:,pcol) = L*U.
%
%   W = FZMULT(A,V,doTranspose,L,U,pcol,P) uses the precomputed sparse LU
%   factorization of matrix A(1:m,1:m), that is, A1 = A(1:m,1:m) and
%   P*A1(:,pcol) = L*U.
%
%   The nullspace basis matrix Z is not formed explicitly. An implicit
%   representation is used based on the sparse LU factorization of
%   A(1:m,1:m).

%   Copyright 1990-2020 The MathWorks, Inc.

% NOTE: Since these are now private functions, only used by trust-region
% and projected CG sub-problems, we have removed all error checking for
% compatible size, and number of inputs.
%
% The expectation is that the caller will ensure that sizes agree, and the
% right number (and type) of inputs are provided.
%
% Specifically, these are the requirements of the caller
% - A is sparse, m x n
% - n > m
% - V is mv x k
% - mv == (n - m) if doTranspose is false
% - mv == n if doTranspose is true

% Import the warning free mldivide to avoid disabling / enabling the nearly
% singular warnings
import matlab.internal.math.nowarn.mldivide

% Initialization
[m,n] = size(A); 

% Check to see if the factors, reordering, etc. are given. If not,
% compute them.
if nargin < 4
    % A1 is square
    A1 = A(:,1:m);      
    pcol = colamd(A1);
    [L,U,P] = lu(A1(:,pcol));
end

A2 = A(:,m+1:n); 

if ~doTranspose
    W = -A2*V;
    WW = L\P*W;
    W = U\WW;
    WW(pcol,:) = W;
    W = [WW;V];
else
    V1 = V(1:m,:);
    V2 = V(m+1:n,:);
    WW = -V1(pcol,:);
    W = U'\WW;
    WW = L'\W;
    W = P'*WW;
    WW = A2'*W;
    W = WW + V2;
end

% Follow regular rules of sparse/full multiplication.
% A is always sparse. Thus V determines if W is sparse or full.
if ~issparse(V)
    W = full(W);
else
    W = sparse(W);
end
