--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/UMFPACK/UMFPACK/Source/umfpack_solve.c	Fri Feb 25 19:55:28 2005 +0000
@@ -0,0 +1,245 @@
+/* ========================================================================== */
+/* === UMFPACK_solve ======================================================== */
+/* ========================================================================== */
+
+/* -------------------------------------------------------------------------- */
+/* UMFPACK Version 4.4, Copyright (c) 2005 by Timothy A. Davis.  CISE Dept,   */
+/* Univ. of Florida.  All Rights Reserved.  See ../Doc/License for License.   */
+/* web: http://www.cise.ufl.edu/research/sparse/umfpack                       */
+/* -------------------------------------------------------------------------- */
+
+/*
+    User-callable.  Solves a linear system using the numerical factorization
+    computed by UMFPACK_numeric.  See umfpack_solve.h for more details.
+
+    For umfpack_*_solve:
+	Dynamic memory usage:  UMFPACK_solve calls UMF_malloc twice, for
+	workspace of size c*n*sizeof(double) + n*sizeof(Int), where c is
+	defined below.  On return, all of this workspace is free'd via UMF_free.
+
+    For umfpack_*_wsolve:
+	No dynamic memory usage.  Input arrays are used for workspace instead.
+	Pattern is a workspace of size n Integers.  The double array W must be
+	at least of size c*n, where c is defined below.
+
+    If iterative refinement is requested, and Ax=b, A'x=b or A.'x=b is being
+    solved, and the matrix A is not singular, then c is 5 for the real version
+    and 10 for the complex version.  Otherwise, c is 1 for the real version and
+    4 for the complex version.
+*/
+
+#include "umf_internal.h"
+#include "umf_valid_numeric.h"
+#include "umf_solve.h"
+
+#ifndef WSOLVE
+#include "umf_malloc.h"
+#include "umf_free.h"
+#ifndef NDEBUG
+PRIVATE Int init_count ;
+#endif
+#endif
+
+GLOBAL Int
+#ifdef WSOLVE
+UMFPACK_wsolve
+#else
+UMFPACK_solve
+#endif
+(
+    Int sys,
+    const Int Ap [ ],
+    const Int Ai [ ],
+    const double Ax [ ],
+#ifdef COMPLEX
+    const double Az [ ],
+#endif
+    double Xx [ ],
+#ifdef COMPLEX
+    double Xz [ ],
+#endif
+    const double Bx [ ],
+#ifdef COMPLEX
+    const double Bz [ ],
+#endif
+    void *NumericHandle,
+    const double Control [UMFPACK_CONTROL],
+    double User_Info [UMFPACK_INFO]
+#ifdef WSOLVE
+    , Int Pattern [ ],
+    double W [ ]
+#endif
+)
+{
+    /* ---------------------------------------------------------------------- */
+    /* local variables */
+    /* ---------------------------------------------------------------------- */
+
+    double Info2 [UMFPACK_INFO], stats [2] ;
+    double *Info ;
+    NumericType *Numeric ;
+    Int n, i, irstep, status ;
+#ifndef WSOLVE
+    Int *Pattern, wsize ;
+    double *W ;
+#endif
+
+    /* ---------------------------------------------------------------------- */
+    /* get the amount of time used by the process so far */
+    /* ---------------------------------------------------------------------- */
+
+    umfpack_tic (stats) ;
+
+#ifndef WSOLVE
+#ifndef NDEBUG
+    init_count = UMF_malloc_count ;
+#endif
+#endif
+
+    /* ---------------------------------------------------------------------- */
+    /* get parameters */
+    /* ---------------------------------------------------------------------- */
+
+    irstep = GET_CONTROL (UMFPACK_IRSTEP, UMFPACK_DEFAULT_IRSTEP) ;
+
+    if (User_Info != (double *) NULL)
+    {
+	/* return Info in user's array */
+	Info = User_Info ;
+	/* clear the parts of Info that are set by UMFPACK_solve */
+	for (i = UMFPACK_IR_TAKEN ; i <= UMFPACK_SOLVE_TIME ; i++)
+	{
+	    Info [i] = EMPTY ;
+	}
+    }
+    else
+    {
+	/* no Info array passed - use local one instead */
+	Info = Info2 ;
+	for (i = 0 ; i < UMFPACK_INFO ; i++)
+	{
+	    Info [i] = EMPTY ;
+	}
+    }
+
+    Info [UMFPACK_STATUS] = UMFPACK_OK ;
+    Info [UMFPACK_SOLVE_FLOPS] = 0 ;
+
+    Numeric = (NumericType *) NumericHandle ;
+    if (!UMF_valid_numeric (Numeric))
+    {
+	Info [UMFPACK_STATUS] = UMFPACK_ERROR_invalid_Numeric_object ;
+	return (UMFPACK_ERROR_invalid_Numeric_object) ;
+    }
+
+    Info [UMFPACK_NROW] = Numeric->n_row ;
+    Info [UMFPACK_NCOL] = Numeric->n_col ;
+
+    if (Numeric->n_row != Numeric->n_col)
+    {
+	/* only square systems can be handled */
+	Info [UMFPACK_STATUS] = UMFPACK_ERROR_invalid_system ;
+	return (UMFPACK_ERROR_invalid_system) ;
+    }
+    n = Numeric->n_row ;
+    if (Numeric->nnzpiv < n
+	|| SCALAR_IS_ZERO (Numeric->rcond) || SCALAR_IS_NAN (Numeric->rcond))
+    {
+	/* turn off iterative refinement if A is singular */
+	/* or if U has NaN's on the diagonal. */
+	irstep = 0 ;
+    }
+
+    if (!Xx || !Bx)
+    {
+	Info [UMFPACK_STATUS] = UMFPACK_ERROR_argument_missing ;
+	return (UMFPACK_ERROR_argument_missing) ;
+    }
+
+    if (sys >= UMFPACK_Pt_L)
+    {
+	/* no iterative refinement except for nonsingular Ax=b, A'x=b, A.'x=b */
+	irstep = 0 ;
+    }
+
+    /* ---------------------------------------------------------------------- */
+    /* allocate or check the workspace */
+    /* ---------------------------------------------------------------------- */
+
+#ifdef WSOLVE
+
+    if (!W || !Pattern)
+    {
+	Info [UMFPACK_STATUS] = UMFPACK_ERROR_argument_missing ;
+	return (UMFPACK_ERROR_argument_missing) ;
+    }
+
+#else
+
+#ifdef COMPLEX
+    if (irstep > 0)
+    {
+	wsize = 10*n ;		/* W, X, Z, S, Y, B2 */
+    }
+    else
+    {
+	wsize = 4*n ;		/* W, X */
+    }
+#else
+    if (irstep > 0)
+    {
+	wsize = 5*n ;		/* W, Z, S, Y, B2 */
+    }
+    else
+    {
+	wsize = n ;		/* W */
+    }
+#endif
+
+    Pattern = (Int *) UMF_malloc (n, sizeof (Int)) ;
+    W = (double *) UMF_malloc (wsize, sizeof (double)) ;
+    if (!W || !Pattern)
+    {
+	DEBUGm4 (("out of memory: solve work\n")) ;
+	Info [UMFPACK_STATUS] = UMFPACK_ERROR_out_of_memory ;
+	(void) UMF_free ((void *) W) ;
+	(void) UMF_free ((void *) Pattern) ;
+	return (UMFPACK_ERROR_out_of_memory) ;
+    }
+
+#endif	/* WSOLVE */
+
+    /* ---------------------------------------------------------------------- */
+    /* solve the system */
+    /* ---------------------------------------------------------------------- */
+
+    status = UMF_solve (sys, Ap, Ai, Ax, Xx, Bx,
+#ifdef COMPLEX
+	Az, Xz, Bz,
+#endif
+	Numeric, irstep, Info, Pattern, W) ;
+
+    /* ---------------------------------------------------------------------- */
+    /* free the workspace (if allocated) */
+    /* ---------------------------------------------------------------------- */
+
+#ifndef WSOLVE
+    (void) UMF_free ((void *) W) ;
+    (void) UMF_free ((void *) Pattern) ;
+    ASSERT (UMF_malloc_count == init_count) ;
+#endif
+
+    /* ---------------------------------------------------------------------- */
+    /* get the time used by UMFPACK_*solve */
+    /* ---------------------------------------------------------------------- */
+
+    Info [UMFPACK_STATUS] = status ;
+    if (status >= 0)
+    {
+	umfpack_toc (stats) ;
+	Info [UMFPACK_SOLVE_WALLTIME] = stats [0] ;
+	Info [UMFPACK_SOLVE_TIME] = stats [1] ;
+    }
+
+    return (status) ;
+}