/*! @file ddiagonal.c

 * \brief Auxiliary routines to work with diagonal elements

 *

 * 

 * -- SuperLU routine (version 4.0) --

 * Lawrence Berkeley National Laboratory

 * June 30, 2009

 * 

 */

#include "slu_ddefs.h"

int dfill_diag(int n, NCformat *Astore)

/* fill explicit zeros on the diagonal entries, so that the matrix is not

   structurally singular. */

{

    double *nzval = (double *)Astore->nzval;

    int *rowind = Astore->rowind;

    int *colptr = Astore->colptr;

    int nnz = colptr[n];

    int fill = 0;

    double *nzval_new;

    double zero = 0.0;

    int *rowind_new;

    int i, j, diag;

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

    {

	diag = -1;

	for (j = colptr[i]; j < colptr[i + 1]; j++)

	    if (rowind[j] == i) diag = j;

	if (diag < 0) fill++;

    }

    if (fill)

    {

	nzval_new = doubleMalloc(nnz + fill);

	rowind_new = intMalloc(nnz + fill);

	fill = 0;

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

	{

	    diag = -1;

	    for (j = colptr[i] - fill; j < colptr[i + 1]; j++)

	    {

		if ((rowind_new[j + fill] = rowind[j]) == i) diag = j;

		nzval_new[j + fill] = nzval[j];

	    }

	    if (diag < 0)

	    {

		rowind_new[colptr[i + 1] + fill] = i;

		nzval_new[colptr[i + 1] + fill] = zero;

		fill++;

	    }

	    colptr[i + 1] += fill;

	}

	Astore->nzval = nzval_new;

	Astore->rowind = rowind_new;

	SUPERLU_FREE(nzval);

	SUPERLU_FREE(rowind);

    }

    Astore->nnz += fill;

    return fill;

}

int ddominate(int n, NCformat *Astore)

/* make the matrix diagonally dominant */

{

    double *nzval = (double *)Astore->nzval;

    int *rowind = Astore->rowind;

    int *colptr = Astore->colptr;

    int nnz = colptr[n];

    int fill = 0;

    double *nzval_new;

    int *rowind_new;

    int i, j, diag;

    double s;

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

    {

	diag = -1;

	for (j = colptr[i]; j < colptr[i + 1]; j++)

	    if (rowind[j] == i) diag = j;

	if (diag < 0) fill++;

    }

    if (fill)

    {

	nzval_new = doubleMalloc(nnz + fill);

	rowind_new = intMalloc(nnz+ fill);

	fill = 0;

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

	{

	    s = 1e-6;

	    diag = -1;

	    for (j = colptr[i] - fill; j < colptr[i + 1]; j++)

	    {

		if ((rowind_new[j + fill] = rowind[j]) == i) diag = j;

		s += fabs(nzval_new[j + fill] = nzval[j]);

	    }

	    if (diag >= 0) {

		nzval_new[diag+fill] = s * 3.0;

	    } else {

		rowind_new[colptr[i + 1] + fill] = i;

		nzval_new[colptr[i + 1] + fill] = s * 3.0;

		fill++;

	    }

	    colptr[i + 1] += fill;

	}

	Astore->nzval = nzval_new;

	Astore->rowind = rowind_new;

	SUPERLU_FREE(nzval);

	SUPERLU_FREE(rowind);

    }

    else

    {

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

	{

	    s = 1e-6;

	    diag = -1;

	    for (j = colptr[i]; j < colptr[i + 1]; j++)

	    {

		if (rowind[j] == i) diag = j;

		s += fabs(nzval[j]);

	    }

	    nzval[diag] = s * 3.0;

	}

    }

    Astore->nnz += fill;

    return fill;

}