kusano 7d535a
kusano 7d535a
/*! @file sdiagonal.c
kusano 7d535a
 * \brief Auxiliary routines to work with diagonal elements
kusano 7d535a
 *
kusano 7d535a
 * 
kusano 7d535a
 * -- SuperLU routine (version 4.0) --
kusano 7d535a
 * Lawrence Berkeley National Laboratory
kusano 7d535a
 * June 30, 2009
kusano 7d535a
 * 
kusano 7d535a
 */
kusano 7d535a
kusano 7d535a
#include "slu_sdefs.h"
kusano 7d535a
kusano 7d535a
int sfill_diag(int n, NCformat *Astore)
kusano 7d535a
/* fill explicit zeros on the diagonal entries, so that the matrix is not
kusano 7d535a
   structurally singular. */
kusano 7d535a
{
kusano 7d535a
    float *nzval = (float *)Astore->nzval;
kusano 7d535a
    int *rowind = Astore->rowind;
kusano 7d535a
    int *colptr = Astore->colptr;
kusano 7d535a
    int nnz = colptr[n];
kusano 7d535a
    int fill = 0;
kusano 7d535a
    float *nzval_new;
kusano 7d535a
    float zero = 0.0;
kusano 7d535a
    int *rowind_new;
kusano 7d535a
    int i, j, diag;
kusano 7d535a
kusano 7d535a
    for (i = 0; i < n; i++)
kusano 7d535a
    {
kusano 7d535a
	diag = -1;
kusano 7d535a
	for (j = colptr[i]; j < colptr[i + 1]; j++)
kusano 7d535a
	    if (rowind[j] == i) diag = j;
kusano 7d535a
	if (diag < 0) fill++;
kusano 7d535a
    }
kusano 7d535a
    if (fill)
kusano 7d535a
    {
kusano 7d535a
	nzval_new = floatMalloc(nnz + fill);
kusano 7d535a
	rowind_new = intMalloc(nnz + fill);
kusano 7d535a
	fill = 0;
kusano 7d535a
	for (i = 0; i < n; i++)
kusano 7d535a
	{
kusano 7d535a
	    diag = -1;
kusano 7d535a
	    for (j = colptr[i] - fill; j < colptr[i + 1]; j++)
kusano 7d535a
	    {
kusano 7d535a
		if ((rowind_new[j + fill] = rowind[j]) == i) diag = j;
kusano 7d535a
		nzval_new[j + fill] = nzval[j];
kusano 7d535a
	    }
kusano 7d535a
	    if (diag < 0)
kusano 7d535a
	    {
kusano 7d535a
		rowind_new[colptr[i + 1] + fill] = i;
kusano 7d535a
		nzval_new[colptr[i + 1] + fill] = zero;
kusano 7d535a
		fill++;
kusano 7d535a
	    }
kusano 7d535a
	    colptr[i + 1] += fill;
kusano 7d535a
	}
kusano 7d535a
	Astore->nzval = nzval_new;
kusano 7d535a
	Astore->rowind = rowind_new;
kusano 7d535a
	SUPERLU_FREE(nzval);
kusano 7d535a
	SUPERLU_FREE(rowind);
kusano 7d535a
    }
kusano 7d535a
    Astore->nnz += fill;
kusano 7d535a
    return fill;
kusano 7d535a
}
kusano 7d535a
kusano 7d535a
int sdominate(int n, NCformat *Astore)
kusano 7d535a
/* make the matrix diagonally dominant */
kusano 7d535a
{
kusano 7d535a
    float *nzval = (float *)Astore->nzval;
kusano 7d535a
    int *rowind = Astore->rowind;
kusano 7d535a
    int *colptr = Astore->colptr;
kusano 7d535a
    int nnz = colptr[n];
kusano 7d535a
    int fill = 0;
kusano 7d535a
    float *nzval_new;
kusano 7d535a
    int *rowind_new;
kusano 7d535a
    int i, j, diag;
kusano 7d535a
    double s;
kusano 7d535a
kusano 7d535a
    for (i = 0; i < n; i++)
kusano 7d535a
    {
kusano 7d535a
	diag = -1;
kusano 7d535a
	for (j = colptr[i]; j < colptr[i + 1]; j++)
kusano 7d535a
	    if (rowind[j] == i) diag = j;
kusano 7d535a
	if (diag < 0) fill++;
kusano 7d535a
    }
kusano 7d535a
    if (fill)
kusano 7d535a
    {
kusano 7d535a
	nzval_new = floatMalloc(nnz + fill);
kusano 7d535a
	rowind_new = intMalloc(nnz+ fill);
kusano 7d535a
	fill = 0;
kusano 7d535a
	for (i = 0; i < n; i++)
kusano 7d535a
	{
kusano 7d535a
	    s = 1e-6;
kusano 7d535a
	    diag = -1;
kusano 7d535a
	    for (j = colptr[i] - fill; j < colptr[i + 1]; j++)
kusano 7d535a
	    {
kusano 7d535a
		if ((rowind_new[j + fill] = rowind[j]) == i) diag = j;
kusano 7d535a
		s += fabs(nzval_new[j + fill] = nzval[j]);
kusano 7d535a
	    }
kusano 7d535a
	    if (diag >= 0) {
kusano 7d535a
		nzval_new[diag+fill] = s * 3.0;
kusano 7d535a
	    } else {
kusano 7d535a
		rowind_new[colptr[i + 1] + fill] = i;
kusano 7d535a
		nzval_new[colptr[i + 1] + fill] = s * 3.0;
kusano 7d535a
		fill++;
kusano 7d535a
	    }
kusano 7d535a
	    colptr[i + 1] += fill;
kusano 7d535a
	}
kusano 7d535a
	Astore->nzval = nzval_new;
kusano 7d535a
	Astore->rowind = rowind_new;
kusano 7d535a
	SUPERLU_FREE(nzval);
kusano 7d535a
	SUPERLU_FREE(rowind);
kusano 7d535a
    }
kusano 7d535a
    else
kusano 7d535a
    {
kusano 7d535a
	for (i = 0; i < n; i++)
kusano 7d535a
	{
kusano 7d535a
	    s = 1e-6;
kusano 7d535a
	    diag = -1;
kusano 7d535a
	    for (j = colptr[i]; j < colptr[i + 1]; j++)
kusano 7d535a
	    {
kusano 7d535a
		if (rowind[j] == i) diag = j;
kusano 7d535a
		s += fabs(nzval[j]);
kusano 7d535a
	    }
kusano 7d535a
	    nzval[diag] = s * 3.0;
kusano 7d535a
	}
kusano 7d535a
    }
kusano 7d535a
    Astore->nnz += fill;
kusano 7d535a
    return fill;
kusano 7d535a
}