#include <stdlib.h> /* for getenv and atoi */
#include "pdsp_defs.h"

void
pdgstrf(pdgstrf_options_t *pdgstrf_options, SuperMatrix *A, int *perm_r,
	SuperMatrix *L, SuperMatrix *U,	Gstat_t *Gstat, int *info)
{
/*
 * -- SuperLU MT routine (version 1.0) --
 * Univ. of California Berkeley, Xerox Palo Alto Research Center,
 * and Lawrence Berkeley National Lab.
 * August 15, 1997
 *
 * Purpose
 * =======
 *
 * pdgstrf() computes an LU factorization of a general sparse nrow-by-ncol
 * matrix A using partial pivoting with row interchanges. The factorization
 * has the form
 *     Pr * A = L * U
 * where Pr is a row permutation matrix, L is lower triangular with unit
 * diagonal elements (lower trapezoidal if A->nrow > A->ncol), and U is
 * upper triangular (upper trapezoidal if A->nrow < A->ncol).
 *
 * Arguments
 * =========
 * 
 * pdgstrf_options (input) pdgstrf_options_t*
 *        The structure defines the parameters to control how the sparse
 *        LU factorization is performed. The following fields must be set 
 *        by the user:
 *
 *        o nprocs (int)
 *          Number of processes to be spawned and used for factorization.
 *
 *        o refact (yes_no_t)
 *          Specifies whether this is first time or subsequent factorization.
 *          = NO:  this factorization is treated as the first one;
 *          = YES: it means that a factorization was performed prior to this
 *                 one. Therefore, this factorization will re-use some
 *                 existing data structures, such as L and U storage, column
 *                 elimination tree, and the symbolic information of the
 *                 Householder matrix.
 *
 *        o panel_size (int)
 *          A panel consists of at most panel_size consecutive columns.
 *
 *        o relax (int)
 *          Degree of relaxing supernodes. If the number of nodes (columns)
 *          in a subtree of the elimination tree is less than relax, this 
 *          subtree is considered as one supernode, regardless of the row
 *          structures of those columns.
 *
 *        o diag_pivot_thresh (double)
 *	    Diagonal pivoting threshold. At step j of Gaussian elimination,
 *          if abs(A_jj) >= diag_pivot_thresh * (max_(i>=j) abs(A_ij)),
 *          use A_jj as pivot. 0 <= diag_pivot_thresh <= 1. The default
 *          value is 1.0, corresponding to partial pivoting.
 *
 *        o usepr (yes_no_t)
 *          Whether the pivoting will use perm_r specified by the user.
 *          = YES: use perm_r; perm_r is input, unchanged on exit.
 *          = NO:  perm_r is determined by partial pivoting, and is output.
 *
 *        o drop_tol (double) (NOT IMPLEMENTED)
 *	    Drop tolerance parameter. At step j of the Gaussian elimination,
 *          if abs(A_ij)/(max_i abs(A_ij)) < drop_tol, drop entry A_ij.
 *          0 <= drop_tol <= 1. The default value of drop_tol is 0,
 *          corresponding to not dropping any entry.
 *
 *        o perm_c (int*)
 *	    Column permutation vector of size A->ncol, which defines the 
 *          permutation matrix Pc; perm_c[i] = j means column i of A is 
 *          in position j in A*Pc.
 *
 *        o perm_r (int*)
 *	    Column permutation vector of size A->nrow.
 *          If pdgstrf_options->usepr = NO, this is an output argument.
 *
 *        o work (void*) of size lwork
 *          User-supplied work space and space for the output data structures.
 *          Not referenced if lwork = 0;
 *
 *        o lwork (int)
 *          Specifies the length of work array.
 *            = 0:  allocate space internally by system malloc;
 *            > 0:  use user-supplied work array of length lwork in bytes,
 *                  returns error if space runs out.
 *            = -1: the routine guesses the amount of space needed without
 *                  performing the factorization, and returns it in
 *                  superlu_memusage->total_needed; no other side effects.
 *
 * A      (input) SuperMatrix*
 *	  Original matrix A, permuted by columns, of dimension
 *        (A->nrow, A->ncol). The type of A can be:
 *        Stype = NCP; Dtype = _D; Mtype = GE.
 *
 * perm_r (input/output) int*, dimension A->nrow
 *        Row permutation vector which defines the permutation matrix Pr,
 *        perm_r[i] = j means row i of A is in position j in Pr*A.
 *        If pdgstrf_options->usepr = NO, perm_r is output argument;
 *        If pdgstrf_options->usepr = YES, the pivoting routine will try 
 *           to use the input perm_r, unless a certain threshold criterion
 *           is violated. In that case, perm_r is overwritten by a new
 *           permutation determined by partial pivoting or diagonal 
 *           threshold pivoting.
 *
 * L      (output) SuperMatrix*
 *        The factor L from the factorization Pr*A=L*U; use compressed row 
 *        subscripts storage for supernodes, i.e., L has type: 
 *        Stype = SCP, Dtype = _D, Mtype = TRLU.
 *
 * U      (output) SuperMatrix*
 *	  The factor U from the factorization Pr*A*Pc=L*U. Use column-wise
 *        storage scheme, i.e., U has types: Stype = NCP, Dtype = _D,
 *        Mtype = TRU.
 *
 * Gstat  (output) Gstat_t*
 *        Record all the statistics about the factorization; 
 *        See Gstat_t structure defined in util.h.
 *
 * info   (output) int*
 *        = 0: successful exit
 *        < 0: if info = -i, the i-th argument had an illegal value
 *        > 0: if info = i, and i is
 *             <= A->ncol: U(i,i) is exactly zero. The factorization has
 *                been completed, but the factor U is exactly singular,
 *                and division by zero will occur if it is used to solve a
 *                system of equations.
 *             > A->ncol: number of bytes allocated when memory allocation
 *                failure occurred, plus A->ncol.
 *
 */
    pdgstrf_threadarg_t *pdgstrf_threadarg;
    pxgstrf_shared_t pxgstrf_shared;
    register int nprocs = pdgstrf_options->nprocs;
    register int i, iinfo;
    double    *utime = Gstat->utime;
    double    usrtime, wtime;
    double    usertimer_();
#if ( MACH==SUN )
    thread_t  *thread_id;
#elif ( MACH==DEC || MACH==PTHREAD )
    pthread_t *thread_id;
    void      *status;
#endif
    void      *pdgstrf_thread(void *);


    /* --------------------------------------------------------------
       Initializes the parallel data structures for pdgstrf_thread().
       --------------------------------------------------------------*/
    pdgstrf_threadarg = pdgstrf_thread_init(A, L, U, pdgstrf_options,
					    &pxgstrf_shared, Gstat, info);
    if ( *info ) return;

    /* Start timing factorization. */
    usrtime = usertimer_();
    wtime = SuperLU_timer_(); 

    /* ------------------------------------------------------------
       On a SUN multiprocessor system, use Solaris thread.
       ------------------------------------------------------------*/
#if ( MACH==SUN )
    
    /* Create nproc threads for concurrent factorization. */
    thread_id = (thread_t *) SUPERLU_MALLOC(nprocs * sizeof(thread_t));
    
    for (i = 1; i < nprocs; ++i) {
#if ( PRNTlevel==1 )
	printf(".. Create unbound threads: i %d, nprocs %d\n", i, nprocs);
#endif
	if ( (iinfo = thr_create(NULL, 0,
			       pdgstrf_thread, &(pdgstrf_threadarg[i]),
			       0, &thread_id[i])) )
	{
	    fprintf(stderr, "thr_create: %d\n", iinfo);
	    ABORT("thr_creat()");
	}
    }
	 
    pdgstrf_thread( &(pdgstrf_threadarg[0]) );
    
    /* Wait for all threads to terminate. */
    for (i = 1; i < nprocs; i++) thr_join(thread_id[i], 0, 0);
    SUPERLU_FREE (thread_id);
/* _SOLARIS_2 */


    /* ------------------------------------------------------------
       On a DEC multiprocessor system, use pthread.
       ------------------------------------------------------------*/
#elif ( MACH==DEC )	/* Use DECthreads ... */
    
    /* Create nproc threads for concurrent factorization. */
    thread_id = (pthread_t *) SUPERLU_MALLOC(nprocs * sizeof(pthread_t));
    
    for (i = 0; i < nprocs; ++i) {
	if ( iinfo = pthread_create(&thread_id[i],
				    NULL,
				    pdgstrf_thread, 
				    &(pdgstrf_threadarg[i])) ) {
	    fprintf(stderr, "pthread_create: %d\n", iinfo);
	    ABORT("pthread_create()");
	}
	/*	pthread_bind_to_cpu_np(thread_id[i], i);*/
    }
	 
    /* Wait for all threads to terminate. */
    for (i = 0; i < nprocs; i++)
	pthread_join(thread_id[i], &status);
    SUPERLU_FREE (thread_id);
/* _DEC */

    
    /* ------------------------------------------------------------
       On a SGI Power Challenge or Origin multiprocessor system,
       use parallel C.
       ------------------------------------------------------------*/
#elif ( MACH==SGI || MACH==ORIGIN )	/* Use parallel C ... */

    if ( getenv("MP_SET_NUMTHREADS") ) {
        i = atoi(getenv("MP_SET_NUMTHREADS"));
	if ( nprocs > i ) {
	    printf("nprocs=%d > environment allowed: MP_SET_NUMTHREADS=%d\n",
		   nprocs, i);
	    exit(-1);
	}
    }
#pragma parallel
#pragma shared  (pdgstrf_threadarg)
/*#pragma numthreads (max = nprocs)*/
#pragma numthreads (nprocs)
    {
      pdgstrf_thread( pdgstrf_threadarg );
    }
/* _SGI or _ORIGIN */

    /* ------------------------------------------------------------
       On a Cray PVP multiprocessor system, use microtasking.
       ------------------------------------------------------------*/
#elif ( MACH==CRAY_PVP )       /* Use C microtasking. */

    if ( getenv("NCPUS") ) {
        i = atoi(getenv("NCPUS"));
	if ( nprocs > i ) {
	    printf("nprocs=%d > environment allowed: NCPUS=%d\n",
		   nprocs, i);
	    exit(-1);
	}
    }
#pragma _CRI taskloop private (i,nprocs)  shared (pdgstrf_threadarg)
    /* Stand-alone task loop */
    for (i = 0; i < nprocs; ++i) {
	pdgstrf_thread( &(pdgstrf_threadarg[i]) );
    }
/* _CRAY_PVP */

    /* ------------------------------------------------------------
       Use POSIX threads.
       ------------------------------------------------------------*/
#elif ( MACH==PTHREAD )	/* Use pthread ... */
    
    /* Create nproc threads for concurrent factorization. */
    thread_id = (pthread_t *) SUPERLU_MALLOC(nprocs * sizeof(pthread_t));
    
    for (i = 0; i < nprocs; ++i) {
	if ( iinfo = pthread_create(&thread_id[i],
				    NULL,
				    pdgstrf_thread, 
				    &(pdgstrf_threadarg[i])) ) {
	    fprintf(stderr, "pthread_create: %d\n", iinfo);
	    ABORT("pthread_create()");
	}
    }
	 
    /* Wait for all threads to terminate. */
    for (i = 0; i < nprocs; i++)
	pthread_join(thread_id[i], &status);
    SUPERLU_FREE (thread_id);
/* _PTHREAD */

    /* ------------------------------------------------------------
       On all other systems, use single processor.
       ------------------------------------------------------------*/
#else

    printf("pdgstrf() is not parallelized on this machine.\n");
    printf("pdgstrf() will be run on single processor.\n");
    pdgstrf_thread( &(pdgstrf_threadarg[0]) );
    
#endif    

    wtime = SuperLU_timer_() - wtime;
    usrtime = usertimer_() - usrtime;
    utime[FACT] = wtime;

#if ( PRNTlevel==1 )
    printf(".. pdgstrf_thread() returns info %d, usrtime %.2f, wtime %.2f\n",
           *info, usrtime, wtime);
#endif

    /* check_mem_leak("after pdgstrf_thread()"); */
    
    /* ------------------------------------------------------------
       Clean up and free storage after multithreaded factorization.
       ------------------------------------------------------------*/
    pdgstrf_thread_finalize(pdgstrf_threadarg, &pxgstrf_shared, 
			    A, perm_r, L, U);

}