Mercurial > octave

view liboctave/numeric/sparse-lu.cc @ 21691:263d18409fdf

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Eliminate unused variable warnings for conditionally compiled code. We had more or less decided not to bother trying to eliminate all these warnings for cases in which external dependencies are missing. But then we get people trying to fix these in various ways, so we might as well do it for all cases and use a consistent method. * oct-conf-post.in.h (octave_unused_parameter): New function for C++ code and new macro for C code. * mk-octave-config-h.sh: Emit octave_unused_parameter function and macro for octave-config.h. * CSparse.cc, __delaunayn__.cc, __eigs__.cc, __fltk_uigetfile__.cc, __glpk__.cc, __magick_read__.cc, __osmesa_print__.cc, __voronoi__.cc, amd.cc, audiodevinfo.cc, audioread.cc, ccolamd.cc, cdisplay.c, colamd.cc, convhulln.cc, dSparse.cc, dmperm.cc, fftw.cc, gl-render.cc, lo-error.c, load-save.cc, ls-hdf5.cc, ls-mat5.cc, oct-hdf5-types.cc, ov-base-int.cc, ov-bool-mat.cc, ov-bool-sparse.cc, ov-bool.cc, ov-cell.cc, ov-class.cc, ov-complex.cc, ov-cx-mat.cc, ov-cx-sparse.cc, ov-fcn-handle.cc, ov-fcn-inline.cc, ov-float.cc, ov-flt-complex.cc, ov-flt-cx-mat.cc, ov-flt-re-mat.cc, ov-java.cc, ov-range.cc, ov-re-mat.cc, ov-re-sparse.cc, ov-scalar.cc, ov-str-mat.cc, ov-struct.cc, sparse-chol.cc, sparse-dmsolve.cc, sparse-lu.cc, sparse-qr.cc, sparse-util.cc, symbfact.cc: Use octave_unused_parameter to eliminate warnings for conditionally compiled code.
author John W. Eaton <jwe@octave.org>
date Fri, 13 May 2016 09:36:14 -0400
parents 40de9f8f23a6
children aba2e6293dd8
line wrap: on
line source

/*

Copyright (C) 2016 John W. Eaton
Copyright (C) 2004-2015 David Bateman
Copyright (C) 1998-2004 Andy Adler

This file is part of Octave.

Octave is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 3 of the License, or (at your
option) any later version.

Octave is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
for more details.

You should have received a copy of the GNU General Public License
along with Octave; see the file COPYING.  If not, see
<http://www.gnu.org/licenses/>.

*/

#ifdef HAVE_CONFIG_H
#  include "config.h"
#endif

#include "CSparse.h"
#include "PermMatrix.h"
#include "dSparse.h"
#include "lo-error.h"
#include "oct-locbuf.h"
#include "oct-sparse.h"
#include "oct-spparms.h"
#include "sparse-lu.h"

// Wrappers for SuiteSparse (formerly UMFPACK) functions that have
// different names depending on the sparse matrix data type.
//
// All of these functions must be specialized to forward to the correct
// SuiteSparse functions.

template <typename T>
void
umfpack_defaults (double *Control);

template <typename T>
void
umfpack_free_numeric (void **Numeric);

template <typename T>
void
umfpack_free_symbolic (void **Symbolic);

template <typename T>
octave_idx_type
umfpack_get_lunz (octave_idx_type *lnz, octave_idx_type *unz, void *Numeric);

template <typename T>
octave_idx_type
umfpack_get_numeric (octave_idx_type *Lp, octave_idx_type *Lj,
                     T *Lx, // Or Lz_packed
                     octave_idx_type *Up, octave_idx_type *Ui,
                     T *Ux, // Or Uz_packed
                     octave_idx_type *p, octave_idx_type *q,
                     double *Dz_packed, octave_idx_type *do_recip,
                     double *Rs, void *Numeric);

template <typename T>
octave_idx_type
umfpack_numeric (const octave_idx_type *Ap, const octave_idx_type *Ai,
                 const T *Ax, // Or Az_packed
                 void *Symbolic, void **Numeric,
                 const double *Control, double *Info);

template <typename T>
octave_idx_type
umfpack_qsymbolic (octave_idx_type n_row, octave_idx_type n_col,
                   const octave_idx_type *Ap, const octave_idx_type *Ai,
                   const T *Ax, // Or Az_packed
                   const octave_idx_type *Qinit, void **Symbolic,
                   const double *Control, double *Info);

template <typename T>
void
umfpack_report_control (const double *Control);

template <typename T>
void
umfpack_report_info (const double *Control, const double *Info);

template <typename T>
void
umfpack_report_matrix (octave_idx_type n_row, octave_idx_type n_col,
                       const octave_idx_type *Ap, const octave_idx_type *Ai,
                       const T *Ax, // Or Az_packed
                       octave_idx_type col_form, const double *Control);

template <typename T>
void
umfpack_report_numeric (void *Numeric, const double *Control);

template <typename T>
void
umfpack_report_perm (octave_idx_type np, const octave_idx_type *Perm,
                     const double *Control);

template <typename T>
void
umfpack_report_status (double *Control, octave_idx_type status);

template <typename T>
void
umfpack_report_symbolic (void *Symbolic, const double *Control);

#if defined (HAVE_UMFPACK)

// SparseMatrix Specialization.

template <>
inline void
umfpack_defaults<double> (double *Control)
{
  UMFPACK_DNAME (defaults) (Control);
}

template <>
inline void
umfpack_free_numeric<double> (void **Numeric)
{
  UMFPACK_DNAME (free_numeric) (Numeric);
}

template <>
inline void
umfpack_free_symbolic<double> (void **Symbolic)
{
  UMFPACK_DNAME (free_symbolic) (Symbolic);
}

template <>
inline octave_idx_type
umfpack_get_lunz<double>
  (octave_idx_type *lnz, octave_idx_type *unz, void *Numeric)
{
  octave_idx_type ignore1, ignore2, ignore3;

  return UMFPACK_DNAME (get_lunz) (lnz, unz, &ignore1, &ignore2,
                                   &ignore3, Numeric);
}

template <>
inline octave_idx_type
umfpack_get_numeric<double>
  (octave_idx_type *Lp, octave_idx_type *Lj, double *Lx,
   octave_idx_type *Up, octave_idx_type *Ui, double *Ux,
   octave_idx_type *p, octave_idx_type *q, double *Dx,
   octave_idx_type *do_recip, double *Rs, void *Numeric)
{
  return UMFPACK_DNAME (get_numeric) (Lp, Lj, Lx, Up, Ui, Ux, p, q, Dx,
                                      do_recip, Rs, Numeric);
}

template <>
inline octave_idx_type
umfpack_numeric<double>
  (const octave_idx_type *Ap, const octave_idx_type *Ai,
   const double *Ax, void *Symbolic, void **Numeric,
   const double *Control, double *Info)
{
  return UMFPACK_DNAME (numeric) (Ap, Ai, Ax, Symbolic, Numeric, Control,
                                  Info);
}

template <>
inline octave_idx_type
umfpack_qsymbolic<double>
  (octave_idx_type n_row, octave_idx_type n_col, const octave_idx_type *Ap,
   const octave_idx_type *Ai, const double *Ax,
   const octave_idx_type *Qinit, void **Symbolic,
   const double *Control, double *Info)
{
  return UMFPACK_DNAME (qsymbolic) (n_row, n_col, Ap, Ai, Ax, Qinit,
                                    Symbolic, Control, Info);
}

template <>
inline void
umfpack_report_control<double> (const double *Control)
{
  UMFPACK_DNAME (report_control) (Control);
}

template <>
inline void
umfpack_report_info<double> (const double *Control, const double *Info)
{
  UMFPACK_DNAME (report_info) (Control, Info);
}

template <>
inline void
umfpack_report_matrix<double>
  (octave_idx_type n_row, octave_idx_type n_col, const octave_idx_type *Ap,
   const octave_idx_type *Ai, const double *Ax, octave_idx_type col_form,
   const double *Control)
{
  UMFPACK_DNAME (report_matrix) (n_row, n_col, Ap, Ai, Ax, col_form, Control);
}

template <>
inline void
umfpack_report_numeric<double> (void *Numeric, const double *Control)
{
  UMFPACK_DNAME (report_numeric) (Numeric, Control);
}

template <>
inline void
umfpack_report_perm<double>
  (octave_idx_type np, const octave_idx_type *Perm, const double *Control)
{
  UMFPACK_DNAME (report_perm) (np, Perm, Control);
}

template <>
inline void
umfpack_report_status<double> (double *Control, octave_idx_type status)
{
  UMFPACK_DNAME (report_status) (Control, status);
}

template <>
inline void
umfpack_report_symbolic<double> (void *Symbolic, const double *Control)
{
  UMFPACK_DNAME (report_symbolic) (Symbolic, Control);
}

// SparseComplexMatrix specialization.

template <>
inline void
umfpack_defaults<Complex> (double *Control)
{
  UMFPACK_ZNAME (defaults) (Control);
}

template <>
inline void
umfpack_free_numeric<Complex> (void **Numeric)
{
  UMFPACK_ZNAME (free_numeric) (Numeric);
}

template <>
inline void
umfpack_free_symbolic<Complex> (void **Symbolic)
{
  UMFPACK_ZNAME (free_symbolic) (Symbolic);
}

template <>
inline octave_idx_type
umfpack_get_lunz<Complex>
  (octave_idx_type *lnz, octave_idx_type *unz, void *Numeric)
{
  octave_idx_type ignore1, ignore2, ignore3;

  return UMFPACK_ZNAME (get_lunz) (lnz, unz, &ignore1, &ignore2,
                                   &ignore3, Numeric);
}

template <>
inline octave_idx_type
umfpack_get_numeric<Complex>
  (octave_idx_type *Lp, octave_idx_type *Lj, Complex *Lz,
   octave_idx_type *Up, octave_idx_type *Ui, Complex *Uz,
   octave_idx_type *p, octave_idx_type *q, double *Dz,
   octave_idx_type *do_recip, double *Rs, void *Numeric)
{
  return UMFPACK_ZNAME (get_numeric) (Lp, Lj,
                                      reinterpret_cast<double *> (Lz),
                                      0, Up, Ui,
                                      reinterpret_cast<double *> (Uz),
                                      0, p, q,
                                      reinterpret_cast<double *> (Dz),
                                      0, do_recip, Rs, Numeric);
}

template <>
inline octave_idx_type
umfpack_numeric<Complex>
  (const octave_idx_type *Ap, const octave_idx_type *Ai,
   const Complex *Az, void *Symbolic, void **Numeric,
   const double *Control, double *Info)
{
  return UMFPACK_ZNAME (numeric) (Ap, Ai,
                                  reinterpret_cast<const double *> (Az),
                                  0, Symbolic, Numeric, Control, Info);
}


template <>
inline octave_idx_type
umfpack_qsymbolic<Complex>
  (octave_idx_type n_row, octave_idx_type n_col,
   const octave_idx_type *Ap, const octave_idx_type *Ai,
   const Complex *Az, const octave_idx_type *Qinit,
   void **Symbolic, const double *Control, double *Info)
{
  return UMFPACK_ZNAME (qsymbolic) (n_row, n_col, Ap, Ai,
                                    reinterpret_cast<const double *> (Az),
                                    0, Qinit, Symbolic, Control, Info);
}


template <>
inline void
umfpack_report_control<Complex> (const double *Control)
{
  UMFPACK_ZNAME (report_control) (Control);
}

template <>
inline void
umfpack_report_info<Complex> (const double *Control, const double *Info)
{
  UMFPACK_ZNAME (report_info) (Control, Info);
}

template <>
inline void
umfpack_report_matrix<Complex>
  (octave_idx_type n_row, octave_idx_type n_col,
   const octave_idx_type *Ap, const octave_idx_type *Ai,
   const Complex *Az, octave_idx_type col_form, const double *Control)
{
  UMFPACK_ZNAME (report_matrix) (n_row, n_col, Ap, Ai,
                                 reinterpret_cast<const double *> (Az),
                                 0, col_form, Control);
}

template <>
inline void
umfpack_report_numeric<Complex> (void *Numeric, const double *Control)
{
  UMFPACK_ZNAME (report_numeric) (Numeric, Control);
}

template <>
inline void
umfpack_report_perm<Complex>
  (octave_idx_type np, const octave_idx_type *Perm, const double *Control)
{
  UMFPACK_ZNAME (report_perm) (np, Perm, Control);
}

template <>
inline void
umfpack_report_status<Complex> (double *Control, octave_idx_type status)
{
  UMFPACK_ZNAME (report_status) (Control, status);
}

template <>
inline void
umfpack_report_symbolic <Complex> (void *Symbolic, const double *Control)
{
  UMFPACK_ZNAME (report_symbolic) (Symbolic, Control);
}

#endif

template <typename lu_type>
sparse_lu<lu_type>::sparse_lu (const lu_type& a, const Matrix& piv_thres,
                               bool scale)
  : Lfact (), Ufact (), Rfact (), cond (0), P (), Q ()
{
#ifdef HAVE_UMFPACK
  octave_idx_type nr = a.rows ();
  octave_idx_type nc = a.cols ();

  // Setup the control parameters
  Matrix Control (UMFPACK_CONTROL, 1);
  double *control = Control.fortran_vec ();
  umfpack_defaults<lu_elt_type> (control);

  double tmp = octave_sparse_params::get_key ("spumoni");
  if (! xisnan (tmp))
    Control (UMFPACK_PRL) = tmp;

  if (piv_thres.numel () == 2)
    {
      tmp = (piv_thres (0) > 1. ? 1. : piv_thres (0));
      if (! xisnan (tmp))
        Control (UMFPACK_PIVOT_TOLERANCE) = tmp;

      tmp = (piv_thres (1) > 1. ? 1. : piv_thres (1));
      if (! xisnan (tmp))
        Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
    }
  else
    {
      tmp = octave_sparse_params::get_key ("piv_tol");
      if (! xisnan (tmp))
        Control (UMFPACK_PIVOT_TOLERANCE) = tmp;

      tmp = octave_sparse_params::get_key ("sym_tol");
      if (! xisnan (tmp))
        Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
    }

  // Set whether we are allowed to modify Q or not
  tmp = octave_sparse_params::get_key ("autoamd");
  if (! xisnan (tmp))
    Control (UMFPACK_FIXQ) = tmp;

  // Turn-off UMFPACK scaling for LU
  if (scale)
    Control (UMFPACK_SCALE) = UMFPACK_SCALE_SUM;
  else
    Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;

  umfpack_report_control<lu_elt_type> (control);

  const octave_idx_type *Ap = a.cidx ();
  const octave_idx_type *Ai = a.ridx ();
  const lu_elt_type *Ax = a.data ();

  umfpack_report_matrix<lu_elt_type> (nr, nc, Ap, Ai, Ax, static_cast<octave_idx_type> (1), control);

  void *Symbolic;
  Matrix Info (1, UMFPACK_INFO);
  double *info = Info.fortran_vec ();
  int status = umfpack_qsymbolic<lu_elt_type> (nr, nc, Ap, Ai, Ax, 0, &Symbolic, control, info);

  if (status < 0)
    {
      umfpack_report_status<lu_elt_type> (control, status);
      umfpack_report_info<lu_elt_type> (control, info);

      umfpack_free_symbolic<lu_elt_type> (&Symbolic);

      (*current_liboctave_error_handler)
        ("sparse_lu: symbolic factorization failed");
    }
  else
    {
      umfpack_report_symbolic<lu_elt_type> (Symbolic, control);

      void *Numeric;
      status = umfpack_numeric<lu_elt_type> (Ap, Ai, Ax, Symbolic, &Numeric, control, info);
      umfpack_free_symbolic<lu_elt_type> (&Symbolic);

      cond = Info (UMFPACK_RCOND);

      if (status < 0)
        {
          umfpack_report_status<lu_elt_type> (control, status);
          umfpack_report_info<lu_elt_type> (control, info);

          umfpack_free_numeric<lu_elt_type> (&Numeric);

          (*current_liboctave_error_handler)
            ("sparse_lu: numeric factorization failed");
        }
      else
        {
          umfpack_report_numeric<lu_elt_type> (Numeric, control);

          octave_idx_type lnz, unz;
          status = umfpack_get_lunz<lu_elt_type> (&lnz, &unz, Numeric);

          if (status < 0)
            {
              umfpack_report_status<lu_elt_type> (control, status);
              umfpack_report_info<lu_elt_type> (control, info);

              umfpack_free_numeric<lu_elt_type> (&Numeric);

              (*current_liboctave_error_handler)
                ("sparse_lu: extracting LU factors failed");
            }
          else
            {
              octave_idx_type n_inner = (nr < nc ? nr : nc);

              if (lnz < 1)
                Lfact = lu_type (n_inner, nr, static_cast<octave_idx_type> (1));
              else
                Lfact = lu_type (n_inner, nr, lnz);

              octave_idx_type *Ltp = Lfact.cidx ();
              octave_idx_type *Ltj = Lfact.ridx ();
              lu_elt_type *Ltx = Lfact.data ();

              if (unz < 1)
                Ufact = lu_type (n_inner, nc, static_cast<octave_idx_type> (1));
              else
                Ufact = lu_type (n_inner, nc, unz);

              octave_idx_type *Up = Ufact.cidx ();
              octave_idx_type *Uj = Ufact.ridx ();
              lu_elt_type *Ux = Ufact.data ();

              Rfact = SparseMatrix (nr, nr, nr);
              for (octave_idx_type i = 0; i < nr; i++)
                {
                  Rfact.xridx (i) = i;
                  Rfact.xcidx (i) = i;
                }
              Rfact.xcidx (nr) = nr;
              double *Rx = Rfact.data ();

              P.resize (dim_vector (nr, 1));
              octave_idx_type *p = P.fortran_vec ();

              Q.resize (dim_vector (nc, 1));
              octave_idx_type *q = Q.fortran_vec ();

              octave_idx_type do_recip;
              status = umfpack_get_numeric<lu_elt_type> (Ltp, Ltj, Ltx, Up, Uj, Ux, p, q, 0, &do_recip, Rx, Numeric);

              umfpack_free_numeric<lu_elt_type> (&Numeric);

              if (status < 0)
                {
                  umfpack_report_status<lu_elt_type> (control, status);

                  (*current_liboctave_error_handler)
                    ("sparse_lu: extracting LU factors failed");
                }
              else
                {
                  Lfact = Lfact.transpose ();

                  if (do_recip)
                    for (octave_idx_type i = 0; i < nr; i++)
                      Rx[i] = 1.0 / Rx[i];

                  umfpack_report_matrix<lu_elt_type> (nr, n_inner, Lfact.cidx (), Lfact.ridx (), Lfact.data (), static_cast<octave_idx_type> (1), control);
                  umfpack_report_matrix<lu_elt_type> (n_inner, nc, Ufact.cidx (), Ufact.ridx (), Ufact.data (), static_cast<octave_idx_type> (1), control);
                  umfpack_report_perm<lu_elt_type> (nr, p, control);
                  umfpack_report_perm<lu_elt_type> (nc, q, control);
                }

              umfpack_report_info<lu_elt_type> (control, info);
            }
        }
    }

#else

  octave_unused_parameter (a);
  octave_unused_parameter (piv_thres);
  octave_unused_parameter (scale);

  (*current_liboctave_error_handler)
    ("support for UMFPACK was unavailable or disabled when liboctave was built");

#endif
}

template <typename lu_type>
sparse_lu<lu_type>::sparse_lu (const lu_type& a,
                               const ColumnVector& Qinit,
                               const Matrix& piv_thres, bool scale,
                               bool FixedQ, double droptol,
                               bool milu, bool udiag)
  : Lfact (), Ufact (), Rfact (), cond (0), P (), Q ()
{
#ifdef HAVE_UMFPACK

  if (milu)
    (*current_liboctave_error_handler)
      ("Modified incomplete LU not implemented");

  octave_idx_type nr = a.rows ();
  octave_idx_type nc = a.cols ();

  // Setup the control parameters
  Matrix Control (UMFPACK_CONTROL, 1);
  double *control = Control.fortran_vec ();
  umfpack_defaults<lu_elt_type> (control);

  double tmp = octave_sparse_params::get_key ("spumoni");
  if (! xisnan (tmp))
    Control (UMFPACK_PRL) = tmp;

  if (piv_thres.numel () == 2)
    {
      tmp = (piv_thres (0) > 1. ? 1. : piv_thres (0));
      if (! xisnan (tmp))
        Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
      tmp = (piv_thres (1) > 1. ? 1. : piv_thres (1));
      if (! xisnan (tmp))
        Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
    }
  else
    {
      tmp = octave_sparse_params::get_key ("piv_tol");
      if (! xisnan (tmp))
        Control (UMFPACK_PIVOT_TOLERANCE) = tmp;

      tmp = octave_sparse_params::get_key ("sym_tol");
      if (! xisnan (tmp))
        Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
    }

  if (droptol >= 0.)
    Control (UMFPACK_DROPTOL) = droptol;

  // Set whether we are allowed to modify Q or not
  if (FixedQ)
    Control (UMFPACK_FIXQ) = 1.0;
  else
    {
      tmp = octave_sparse_params::get_key ("autoamd");
      if (! xisnan (tmp))
        Control (UMFPACK_FIXQ) = tmp;
    }

  // Turn-off UMFPACK scaling for LU
  if (scale)
    Control (UMFPACK_SCALE) = UMFPACK_SCALE_SUM;
  else
    Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;

  umfpack_report_control<lu_elt_type> (control);

  const octave_idx_type *Ap = a.cidx ();
  const octave_idx_type *Ai = a.ridx ();
  const lu_elt_type *Ax = a.data ();

  umfpack_report_matrix<lu_elt_type> (nr, nc, Ap, Ai, Ax, static_cast<octave_idx_type> (1), control);

  void *Symbolic;
  Matrix Info (1, UMFPACK_INFO);
  double *info = Info.fortran_vec ();
  int status;

  // Null loop so that qinit is imediately deallocated when not needed
  do
    {
      OCTAVE_LOCAL_BUFFER (octave_idx_type, qinit, nc);

      for (octave_idx_type i = 0; i < nc; i++)
        qinit[i] = static_cast<octave_idx_type> (Qinit (i));

      status = umfpack_qsymbolic<lu_elt_type> (nr, nc, Ap, Ai, Ax, qinit, &Symbolic, control, info);
    }
  while (0);

  if (status < 0)
    {
      umfpack_report_status<lu_elt_type> (control, status);
      umfpack_report_info<lu_elt_type> (control, info);

      umfpack_free_symbolic<lu_elt_type> (&Symbolic);

      (*current_liboctave_error_handler)
        ("sparse_lu: symbolic factorization failed");
    }
  else
    {
      umfpack_report_symbolic<lu_elt_type> (Symbolic, control);

      void *Numeric;
      status = umfpack_numeric<lu_elt_type> (Ap, Ai, Ax, Symbolic, &Numeric, control, info);
      umfpack_free_symbolic<lu_elt_type> (&Symbolic);

      cond = Info (UMFPACK_RCOND);

      if (status < 0)
        {
          umfpack_report_status<lu_elt_type> (control, status);
          umfpack_report_info<lu_elt_type> (control, info);

          umfpack_free_numeric<lu_elt_type> (&Numeric);

          (*current_liboctave_error_handler)
            ("sparse_lu: numeric factorization failed");
        }
      else
        {
          umfpack_report_numeric<lu_elt_type> (Numeric, control);

          octave_idx_type lnz, unz;
          status = umfpack_get_lunz<lu_elt_type> (&lnz, &unz, Numeric);

          if (status < 0)
            {
              umfpack_report_status<lu_elt_type> (control, status);
              umfpack_report_info<lu_elt_type> (control, info);

              umfpack_free_numeric<lu_elt_type> (&Numeric);

              (*current_liboctave_error_handler)
                ("sparse_lu: extracting LU factors failed");
            }
          else
            {
              octave_idx_type n_inner = (nr < nc ? nr : nc);

              if (lnz < 1)
                Lfact = lu_type (n_inner, nr, static_cast<octave_idx_type> (1));
              else
                Lfact = lu_type (n_inner, nr, lnz);

              octave_idx_type *Ltp = Lfact.cidx ();
              octave_idx_type *Ltj = Lfact.ridx ();
              lu_elt_type *Ltx = Lfact.data ();

              if (unz < 1)
                Ufact = lu_type (n_inner, nc, static_cast<octave_idx_type> (1));
              else
                Ufact = lu_type (n_inner, nc, unz);

              octave_idx_type *Up = Ufact.cidx ();
              octave_idx_type *Uj = Ufact.ridx ();
              lu_elt_type *Ux = Ufact.data ();

              Rfact = SparseMatrix (nr, nr, nr);
              for (octave_idx_type i = 0; i < nr; i++)
                {
                  Rfact.xridx (i) = i;
                  Rfact.xcidx (i) = i;
                }
              Rfact.xcidx (nr) = nr;
              double *Rx = Rfact.data ();

              P.resize (dim_vector (nr, 1));
              octave_idx_type *p = P.fortran_vec ();

              Q.resize (dim_vector (nc, 1));
              octave_idx_type *q = Q.fortran_vec ();

              octave_idx_type do_recip;
              status = umfpack_get_numeric<lu_elt_type> (Ltp, Ltj, Ltx, Up, Uj, Ux, p, q, 0, &do_recip, Rx, Numeric);

              umfpack_free_numeric<lu_elt_type> (&Numeric);

              if (status < 0)
                {
                  umfpack_report_status<lu_elt_type> (control, status);

                  (*current_liboctave_error_handler)
                    ("sparse_lu: extracting LU factors failed");
                }
              else
                {
                  Lfact = Lfact.transpose ();

                  if (do_recip)
                    for (octave_idx_type i = 0; i < nr; i++)
                      Rx[i] = 1.0 / Rx[i];

                  umfpack_report_matrix<lu_elt_type> (nr, n_inner, Lfact.cidx (), Lfact.ridx (), Lfact.data (), static_cast<octave_idx_type> (1), control);
                  umfpack_report_matrix<lu_elt_type> (n_inner, nc, Ufact.cidx (), Ufact.ridx (), Ufact.data (), static_cast<octave_idx_type> (1), control);
                  umfpack_report_perm<lu_elt_type> (nr, p, control);
                  umfpack_report_perm<lu_elt_type> (nc, q, control);
                }

              umfpack_report_info<lu_elt_type> (control, info);
            }
        }
    }

  if (udiag)
    (*current_liboctave_error_handler)
      ("Option udiag of incomplete LU not implemented");

#else

  octave_unused_parameter (a);
  octave_unused_parameter (Qinit);
  octave_unused_parameter (piv_thres);
  octave_unused_parameter (scale);
  octave_unused_parameter (FixedQ);
  octave_unused_parameter (droptol);
  octave_unused_parameter (milu);
  octave_unused_parameter (udiag);

  (*current_liboctave_error_handler)
    ("support for UMFPACK was unavailable or disabled when liboctave was built");

#endif
}

template <typename lu_type>
lu_type
sparse_lu<lu_type>::Y (void) const
{
  octave_idx_type nr = Lfact.rows ();
  octave_idx_type nz = Lfact.cols ();
  octave_idx_type nc = Ufact.cols ();

  lu_type Yout (nr, nc, Lfact.nnz () + Ufact.nnz () - (nr<nz?nr:nz));
  octave_idx_type ii = 0;
  Yout.xcidx (0) = 0;

  for (octave_idx_type j = 0; j < nc; j++)
    {
      for (octave_idx_type i = Ufact.cidx (j); i < Ufact.cidx (j + 1); i++)
        {
          Yout.xridx (ii) = Ufact.ridx (i);
          Yout.xdata (ii++) = Ufact.data (i);
        }

      if (j < nz)
        {
          // Note the +1 skips the 1.0 on the diagonal
          for (octave_idx_type i = Lfact.cidx (j) + 1;
               i < Lfact.cidx (j +1); i++)
            {
              Yout.xridx (ii) = Lfact.ridx (i);
              Yout.xdata (ii++) = Lfact.data (i);
            }
        }

      Yout.xcidx (j + 1) = ii;
    }

  return Yout;
}

template <typename lu_type>
SparseMatrix
sparse_lu<lu_type>::Pr (void) const
{
  octave_idx_type nr = Lfact.rows ();

  SparseMatrix Pout (nr, nr, nr);

  for (octave_idx_type i = 0; i < nr; i++)
    {
      Pout.cidx (i) = i;
      Pout.ridx (P (i)) = i;
      Pout.data (i) = 1;
    }

  Pout.cidx (nr) = nr;

  return Pout;
}

template <typename lu_type>
ColumnVector
sparse_lu<lu_type>::Pr_vec (void) const
{
  octave_idx_type nr = Lfact.rows ();

  ColumnVector Pout (nr);

  for (octave_idx_type i = 0; i < nr; i++)
    Pout.xelem (i) = static_cast<double> (P(i) + 1);

  return Pout;
}

template <typename lu_type>
PermMatrix
sparse_lu<lu_type>::Pr_mat (void) const
{
  return PermMatrix (P, false);
}

template <typename lu_type>
SparseMatrix
sparse_lu<lu_type>::Pc (void) const
{
  octave_idx_type nc = Ufact.cols ();

  SparseMatrix Pout (nc, nc, nc);

  for (octave_idx_type i = 0; i < nc; i++)
    {
      Pout.cidx (i) = i;
      Pout.ridx (i) = Q (i);
      Pout.data (i) = 1;
    }

  Pout.cidx (nc) = nc;

  return Pout;
}

template <typename lu_type>
ColumnVector
sparse_lu<lu_type>::Pc_vec (void) const
{
  octave_idx_type nc = Ufact.cols ();

  ColumnVector Pout (nc);

  for (octave_idx_type i = 0; i < nc; i++)
    Pout.xelem (i) = static_cast<double> (Q(i) + 1);

  return Pout;
}

template <typename lu_type>
PermMatrix
sparse_lu<lu_type>::Pc_mat (void) const
{
  return PermMatrix (Q, true);
}

// Instantiations we need.

template class sparse_lu<SparseMatrix>;

template class sparse_lu<SparseComplexMatrix>;