/*! @file util.c

 * \brief Utility functions

 * 

 * 

 * -- SuperLU routine (version 4.1) --

 * Univ. of California Berkeley, Xerox Palo Alto Research Center,

 * and Lawrence Berkeley National Lab.

 * November, 2010

 *

 * Copyright (c) 1994 by Xerox Corporation.  All rights reserved.

 *

 * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY

 * EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.

 *

 * Permission is hereby granted to use or copy this program for any

 * purpose, provided the above notices are retained on all copies.

 * Permission to modify the code and to distribute modified code is

 * granted, provided the above notices are retained, and a notice that

 * the code was modified is included with the above copyright notice.

 * 

 */

#include <math.h></math.h>

#include "slu_ddefs.h"

/*! \brief Global statistics variale

 */

void superlu_abort_and_exit(char* msg)

{

    fprintf(stderr, msg);

    exit (-1);

}

/*! \brief Set the default values for the options argument.

 */

void set_default_options(superlu_options_t *options)

{

    options->Fact = DOFACT;

    options->Equil = YES;

    options->ColPerm = COLAMD;

    options->Trans = NOTRANS;

    options->IterRefine = NOREFINE;

    options->DiagPivotThresh = 1.0;

    options->SymmetricMode = NO;

    options->PivotGrowth = NO;

    options->ConditionNumber = NO;

    options->PrintStat = YES;

}

/*! \brief Set the default values for the options argument for ILU.

 */

void ilu_set_default_options(superlu_options_t *options)

{

    set_default_options(options);

    /* further options for incomplete factorization */

    options->DiagPivotThresh = 0.1;

    options->RowPerm = LargeDiag;

    options->ILU_DropRule = DROP_BASIC | DROP_AREA;

    options->ILU_DropTol = 1e-4;

    options->ILU_FillFactor = 10.0;

    options->ILU_Norm = INF_NORM;

    options->ILU_MILU = SILU;

    options->ILU_MILU_Dim = 3.0; /* -log(n)/log(h) is perfect */

    options->ILU_FillTol = 1e-2;

}

/*! \brief Print the options setting.

 */

void print_options(superlu_options_t *options)

{

    printf(".. options:\n");

    printf("\tFact\t %8d\n", options->Fact);

    printf("\tEquil\t %8d\n", options->Equil);

    printf("\tColPerm\t %8d\n", options->ColPerm);

    printf("\tDiagPivotThresh %8.4f\n", options->DiagPivotThresh);

    printf("\tTrans\t %8d\n", options->Trans);

    printf("\tIterRefine\t%4d\n", options->IterRefine);

    printf("\tSymmetricMode\t%4d\n", options->SymmetricMode);

    printf("\tPivotGrowth\t%4d\n", options->PivotGrowth);

    printf("\tConditionNumber\t%4d\n", options->ConditionNumber);

    printf("..\n");

}

/*! \brief Print the options setting.

 */

void print_ilu_options(superlu_options_t *options)

{

    printf(".. ILU options:\n");

    printf("\tDiagPivotThresh\t%6.2e\n", options->DiagPivotThresh);

    printf("\ttau\t%6.2e\n", options->ILU_DropTol);

    printf("\tgamma\t%6.2f\n", options->ILU_FillFactor);

    printf("\tDropRule\t%0x\n", options->ILU_DropRule);

    printf("\tMILU\t%d\n", options->ILU_MILU);

    printf("\tMILU_ALPHA\t%6.2e\n", MILU_ALPHA);

    printf("\tDiagFillTol\t%6.2e\n", options->ILU_FillTol);

    printf("..\n");

}

/*! \brief Deallocate the structure pointing to the actual storage of the matrix. */

void

Destroy_SuperMatrix_Store(SuperMatrix *A)

{

    SUPERLU_FREE ( A->Store );

}

void

Destroy_CompCol_Matrix(SuperMatrix *A)

{

    SUPERLU_FREE( ((NCformat *)A->Store)->rowind );

    SUPERLU_FREE( ((NCformat *)A->Store)->colptr );

    SUPERLU_FREE( ((NCformat *)A->Store)->nzval );

    SUPERLU_FREE( A->Store );

}

void

Destroy_CompRow_Matrix(SuperMatrix *A)

{

    SUPERLU_FREE( ((NRformat *)A->Store)->colind );

    SUPERLU_FREE( ((NRformat *)A->Store)->rowptr );

    SUPERLU_FREE( ((NRformat *)A->Store)->nzval );

    SUPERLU_FREE( A->Store );

}

void

Destroy_SuperNode_Matrix(SuperMatrix *A)

{

    SUPERLU_FREE ( ((SCformat *)A->Store)->rowind );

    SUPERLU_FREE ( ((SCformat *)A->Store)->rowind_colptr );

    SUPERLU_FREE ( ((SCformat *)A->Store)->nzval );

    SUPERLU_FREE ( ((SCformat *)A->Store)->nzval_colptr );

    SUPERLU_FREE ( ((SCformat *)A->Store)->col_to_sup );

    SUPERLU_FREE ( ((SCformat *)A->Store)->sup_to_col );

    SUPERLU_FREE ( A->Store );

}

/*! \brief A is of type Stype==NCP */

void

Destroy_CompCol_Permuted(SuperMatrix *A)

{

    SUPERLU_FREE ( ((NCPformat *)A->Store)->colbeg );

    SUPERLU_FREE ( ((NCPformat *)A->Store)->colend );

    SUPERLU_FREE ( A->Store );

}

/*! \brief A is of type Stype==DN */

void

Destroy_Dense_Matrix(SuperMatrix *A)

{

    DNformat* Astore = A->Store;

    SUPERLU_FREE (Astore->nzval);

    SUPERLU_FREE ( A->Store );

}

/*! \brief Reset repfnz[] for the current column 

 */

void

resetrep_col (const int nseg, const int *segrep, int *repfnz)

{

    int i, irep;

    for (i = 0; i < nseg; i++) {

	irep = segrep[i];

	repfnz[irep] = EMPTY;

    }

}

/*! \brief Count the total number of nonzeros in factors L and U,  and in the symmetrically reduced L. 

 */

void

countnz(const int n, int *xprune, int *nnzL, int *nnzU, GlobalLU_t *Glu)

{

    int          nsuper, fsupc, i, j;

    int          nnzL0, jlen, irep;

    int          *xsup, *xlsub;

    xsup   = Glu->xsup;

    xlsub  = Glu->xlsub;

    *nnzL  = 0;

    *nnzU  = (Glu->xusub)[n];

    nnzL0  = 0;

    nsuper = (Glu->supno)[n];

    if ( n <= 0 ) return;

    /* 

     * For each supernode

     */

    for (i = 0; i <= nsuper; i++) {

	fsupc = xsup[i];

	jlen = xlsub[fsupc+1] - xlsub[fsupc];

	for (j = fsupc; j < xsup[i+1]; j++) {

	    *nnzL += jlen;

	    *nnzU += j - fsupc + 1;

	    jlen--;

	}

	irep = xsup[i+1] - 1;

	nnzL0 += xprune[irep] - xlsub[irep];

    }

    /* printf("\tNo of nonzeros in symm-reduced L = %d\n", nnzL0);*/

}

/*! \brief Count the total number of nonzeros in factors L and U.

 */

void

ilu_countnz(const int n, int *nnzL, int *nnzU, GlobalLU_t *Glu)

{

    int          nsuper, fsupc, i, j;

    int          jlen, irep;

    int          *xsup, *xlsub;

    xsup   = Glu->xsup;

    xlsub  = Glu->xlsub;

    *nnzL  = 0;

    *nnzU  = (Glu->xusub)[n];

    nsuper = (Glu->supno)[n];

    if ( n <= 0 ) return;

    /*

     * For each supernode

     */

    for (i = 0; i <= nsuper; i++) {

	fsupc = xsup[i];

	jlen = xlsub[fsupc+1] - xlsub[fsupc];

	for (j = fsupc; j < xsup[i+1]; j++) {

	    *nnzL += jlen;

	    *nnzU += j - fsupc + 1;

	    jlen--;

	}

	irep = xsup[i+1] - 1;

    }

}

/*! \brief Fix up the data storage lsub for L-subscripts. It removes the subscript sets for structural pruning,	and applies permuation to the remaining subscripts.

 */

void

fixupL(const int n, const int *perm_r, GlobalLU_t *Glu)

{

    register int nsuper, fsupc, nextl, i, j, k, jstrt;

    int          *xsup, *lsub, *xlsub;

    if ( n <= 1 ) return;

    xsup   = Glu->xsup;

    lsub   = Glu->lsub;

    xlsub  = Glu->xlsub;

    nextl  = 0;

    nsuper = (Glu->supno)[n];

    /* 

     * For each supernode ...

     */

    for (i = 0; i <= nsuper; i++) {

	fsupc = xsup[i];

	jstrt = xlsub[fsupc];

	xlsub[fsupc] = nextl;

	for (j = jstrt; j < xlsub[fsupc+1]; j++) {

	    lsub[nextl] = perm_r[lsub[j]]; /* Now indexed into P*A */

	    nextl++;

  	}

	for (k = fsupc+1; k < xsup[i+1]; k++) 

	    	xlsub[k] = nextl;	/* Other columns in supernode i */

    }

    xlsub[n] = nextl;

}

/*! \brief Diagnostic print of segment info after panel_dfs().

 */

void print_panel_seg(int n, int w, int jcol, int nseg, 

		     int *segrep, int *repfnz)

{

    int j, k;

    for (j = jcol; j < jcol+w; j++) {

	printf("\tcol %d:\n", j);

	for (k = 0; k < nseg; k++)

	    printf("\t\tseg %d, segrep %d, repfnz %d\n", k, 

			segrep[k], repfnz[(j-jcol)*n + segrep[k]]);

    }

}

void

StatInit(SuperLUStat_t *stat)

{

    register int i, w, panel_size, relax;

    panel_size = sp_ienv(1);

    relax = sp_ienv(2);

    w = SUPERLU_MAX(panel_size, relax);

    stat->panel_histo = intCalloc(w+1);

    stat->utime = (double *) SUPERLU_MALLOC(NPHASES * sizeof(double));

    if (!stat->utime) ABORT("SUPERLU_MALLOC fails for stat->utime");

    stat->ops = (flops_t *) SUPERLU_MALLOC(NPHASES * sizeof(flops_t));

    if (!stat->ops) ABORT("SUPERLU_MALLOC fails for stat->ops");

    for (i = 0; i < NPHASES; ++i) {

        stat->utime[i] = 0.;

        stat->ops[i] = 0.;

    }

    stat->TinyPivots = 0;

    stat->RefineSteps = 0;

    stat->expansions = 0;

#if ( PRNTlevel >= 1 )

    printf(".. parameters in sp_ienv():\n");

    printf("\t 1: panel size \t %4d \n"

           "\t 2: relax      \t %4d \n"

           "\t 3: max. super \t %4d \n"

           "\t 4: row-dim 2D \t %4d \n"

           "\t 5: col-dim 2D \t %4d \n"

           "\t 6: fill ratio \t %4d \n",

	   sp_ienv(1), sp_ienv(2), sp_ienv(3), 

	   sp_ienv(4), sp_ienv(5), sp_ienv(6));

#endif

}

void

StatPrint(SuperLUStat_t *stat)

{

    double         *utime;

    flops_t        *ops;

    utime = stat->utime;

    ops   = stat->ops;

    printf("Factor time  = %8.2f\n", utime[FACT]);

    if ( utime[FACT] != 0.0 )

      printf("Factor flops = %e\tMflops = %8.2f\n", ops[FACT],

	     ops[FACT]*1e-6/utime[FACT]);

    printf("Solve time   = %8.2f\n", utime[SOLVE]);

    if ( utime[SOLVE] != 0.0 )

      printf("Solve flops = %e\tMflops = %8.2f\n", ops[SOLVE],

	     ops[SOLVE]*1e-6/utime[SOLVE]);

    printf("Number of memory expansions: %d\n", stat->expansions);

}

void

StatFree(SuperLUStat_t *stat)

{

    SUPERLU_FREE(stat->panel_histo);

    SUPERLU_FREE(stat->utime);

    SUPERLU_FREE(stat->ops);

}

flops_t

LUFactFlops(SuperLUStat_t *stat)

{

    return (stat->ops[FACT]);

}

flops_t

LUSolveFlops(SuperLUStat_t *stat)

{

    return (stat->ops[SOLVE]);

}

/*! \brief Fills an integer array with a given value.

 */

void ifill(int *a, int alen, int ival)

{

    register int i;

    for (i = 0; i < alen; i++) a[i] = ival;

}

/*! \brief Get the statistics of the supernodes 

 */

#define NBUCKS 10

static 	int	max_sup_size;

void super_stats(int nsuper, int *xsup)

{

    register int nsup1 = 0;

    int          i, isize, whichb, bl, bh;

    int          bucket[NBUCKS];

    max_sup_size = 0;

    for (i = 0; i <= nsuper; i++) {

	isize = xsup[i+1] - xsup[i];

	if ( isize == 1 ) nsup1++;

	if ( max_sup_size < isize ) max_sup_size = isize;	

    }

    printf("    Supernode statistics:\n\tno of super = %d\n", nsuper+1);

    printf("\tmax supernode size = %d\n", max_sup_size);

    printf("\tno of size 1 supernodes = %d\n", nsup1);

    /* Histogram of the supernode sizes */

    ifill (bucket, NBUCKS, 0);

    for (i = 0; i <= nsuper; i++) {

        isize = xsup[i+1] - xsup[i];

        whichb = (float) isize / max_sup_size * NBUCKS;

        if (whichb >= NBUCKS) whichb = NBUCKS - 1;

        bucket[whichb]++;

    }

    printf("\tHistogram of supernode sizes:\n");

    for (i = 0; i < NBUCKS; i++) {

        bl = (float) i * max_sup_size / NBUCKS;

        bh = (float) (i+1) * max_sup_size / NBUCKS;

        printf("\tsnode: %d-%d\t\t%d\n", bl+1, bh, bucket[i]);

    }

}

float SpaSize(int n, int np, float sum_npw)

{

    return (sum_npw*8 + np*8 + n*4)/1024.;

}

float DenseSize(int n, float sum_nw)

{

    return (sum_nw*8 + n*8)/1024.;;

}

/*! \brief Check whether repfnz[] == EMPTY after reset.

 */

void check_repfnz(int n, int w, int jcol, int *repfnz)

{

    int jj, k;

    for (jj = jcol; jj < jcol+w; jj++) 

	for (k = 0; k < n; k++)

	    if ( repfnz[(jj-jcol)*n + k] != EMPTY ) {

		fprintf(stderr, "col %d, repfnz_col[%d] = %d\n", jj,

			k, repfnz[(jj-jcol)*n + k]);

		ABORT("check_repfnz");

	    }

}

/*! \brief Print a summary of the testing results. */

void

PrintSumm(char *type, int nfail, int nrun, int nerrs)

{

    if ( nfail > 0 )

	printf("%3s driver: %d out of %d tests failed to pass the threshold\n",

	       type, nfail, nrun);

    else

	printf("All tests for %3s driver passed the threshold (%6d tests run)\n", type, nrun);

    if ( nerrs > 0 )

	printf("%6d error messages recorded\n", nerrs);

}

int print_int_vec(char *what, int n, int *vec)

{

    int i;

    printf("%s\n", what);

    for (i = 0; i < n; ++i) printf("%d\t%d\n", i, vec[i]);

    return 0;

}

int slu_PrintInt10(char *name, int len, int *x)

{

    register int i;

    printf("%10s:", name);

    for (i = 0; i < len; ++i)

    {

	if ( i % 10 == 0 ) printf("\n\t[%2d-%2d]", i, i + 9);

	printf("%6d", x[i]);

    }

    printf("\n");

    return 0;

}