/*
 * Copyright (C) 2008-2012 Delft University of Technology - All rights reserved.
 *
 * Author: Reijer Idema
 */

/*
 * PetscMatSeqAij source.
 */
#include "../src/mat/impls/aij/seq/aij.h"	// PETSc MATSEQAIJ declarations
#include "PetscMatSeqAij.h"

PETSCUTILS_NAMESPACE_BEGIN

#undef __FUNCT__
#define __FUNCT__ "PetscUtils::MatSeqAijGetRowCount"
PetscInt MatSeqAijGetRowCount(const _p_Mat* mat)
{
	if (mat == PETSC_NULL) return 0;

#if !PETSC_UTILS_IMPL_TRUECONST
	PetscInt rows, cols;
	MatGetSize(const_cast<Mat>(mat), &rows, &cols);
	return rows;
#else // PETSC_UTILS_IMPL_TRUECONST
	return mat->rmap->N;	// :NOTE: Use rmap->n for local size
#endif

} // MatSeqAijGetRowCount


#undef __FUNCT__
#define __FUNCT__ "PetscUtils::MatSeqAijGetColCount"
PetscInt MatSeqAijGetColCount(const _p_Mat* mat)
{
	if (mat == PETSC_NULL) return 0;

#if !PETSC_UTILS_IMPL_TRUECONST
	PetscInt rows, cols;
	MatGetSize(const_cast<Mat>(mat), &rows, &cols);
	return cols;
#else // PETSC_UTILS_IMPL_TRUECONST
	return mat->cmap->N;	// :NOTE: Use cmap->n for local size
#endif

} // MatSeqAijGetColCount


#undef __FUNCT__
#define __FUNCT__ "PetscUtils::MatSeqAijGetValCount"
PetscInt MatSeqAijGetValCount(const _p_Mat* mat)
{
	if (mat == PETSC_NULL) return 0;

#if !PETSC_UTILS_IMPL_TRUECONST
	MatInfo info;
	MatGetInfo(const_cast<Mat>(mat), MAT_LOCAL, &info);
	return = int(info.nz_used);
#else // PETSC_UTILS_IMPL_TRUECONST
	Mat_SeqAIJ *data = (Mat_SeqAIJ*)mat->data;
	return (data != NULL ? data->nz : 0);
#endif

} // MatSeqAijGetValCount


#undef __FUNCT__
#define __FUNCT__ "PetscUtils::MatSeqAijGetRow"
PetscErrorCode MatSeqAijGetRow(const _p_Mat* mat, PetscInt row, PetscInt* ncols, const PetscInt *cols[], const PetscScalar *vals[])
{

	// Validate input matrix
	if (mat == PETSC_NULL) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"Matrix is NULL");
	if (strcmp(((PetscObject)mat)->type_name, MATSEQAIJ) != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Only for MATSEQAIJ matrices");

#if !PETSC_UTILS_IMPL_TRUECONST
	PetscErrorCode err = MatGetRow(const_cast<Mat>(mat), row, ncols, cols, vals); CHKERRQ(err);
#else // PETSC_UTILS_IMPL_TRUECONST
//	if (!mat->preallocated) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unpreallocated matrix");
	if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
	if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");

	// Set pointers to the matrix data
	Mat_SeqAIJ *a = (Mat_SeqAIJ*)mat->data;
	if (row < 0 || row >= mat->rmap->N) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D out of range", row);
	PetscInt nnz = a->i[row+1] - a->i[row];				// Number of non-zeros in row
	if (ncols) *ncols = nnz;							// Pointer to number of non-zeros in row
	if (cols) *cols = (nnz==0 ? 0 : a->j + a->i[row]);	// Pointer to column indices of non-zeros
	if (vals) *vals = a->a + a->i[row];					// Pointer to values of non-zeros
#endif

	return 0;
} // MatSeqAijGetRow


#undef __FUNCT__
#define __FUNCT__ "PetscUtils::MatSeqAijRestoreRow"
PetscErrorCode MatSeqAijRestoreRow(ConstMat mat, PetscInt row, PetscInt* ncols, const PetscInt *cols[], const PetscScalar *vals[])
{
	// Validate input matrix
	if (mat == PETSC_NULL) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"Matrix is NULL");
	if (strcmp(((PetscObject)mat)->type_name, MATSEQAIJ) != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Only for MATSEQAIJ matrices");

#if !PETSC_UTILS_IMPL_TRUECONST
	PetscErrorCode err = MatRestoreRow(const_cast<Mat>(mat), row, ncols, cols, vals); CHKERRQ(err);
#endif

	return 0;
} // MatSeqAijRestoreRow


#undef __FUNCT__
#define __FUNCT__ "PetscUtils::MatSeqAijView"
PetscErrorCode MatSeqAijView(const _p_Mat* mat)
{
	// Validate input matrix
	if (mat == PETSC_NULL) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"Matrix is NULL");
	if (strcmp(((PetscObject)mat)->type_name, MATSEQAIJ) != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Only for MATSEQAIJ matrices");

#if !PETSC_UTILS_IMPL_TRUECONST
	PetscErrorCode err = MatView(const_cast<Mat>(mat), PETSC_VIEWER_STDOUT_SELF); CHKERRQ(err);
#else // PETSC_UTILS_IMPL_TRUECONST
//	if (!mat->preallocated) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unpreallocated matrix");
	if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
	if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");

	// Print matrix to standard out
	PetscErrorCode err;
	PetscInt rowCount = MatSeqAijGetRowCount(mat);
	PetscInt colCount;
	const PetscInt *pCols;
	const PetscScalar *pVals;
	for (PetscInt i=0; i<rowCount; ++i)
	{
		err = MatSeqAijGetRow(mat, i, &colCount, &pCols, &pVals); CHKERRQ(err);
		PetscPrintf(PETSC_COMM_WORLD, "row %d:", i);
		for (PetscInt j=0; j<colCount; ++j)
		{
			PetscPrintf(PETSC_COMM_WORLD, " (%d, %g) ", pCols[j], pVals[j]);
		}
		PetscPrintf(PETSC_COMM_WORLD, "\n", i);
		err = MatSeqAijRestoreRow(mat, i, &colCount, &pCols, &pVals); CHKERRQ(err);
	}
#endif

	return 0;
} // MatSeqAijView

PETSCUTILS_NAMESPACE_END
