/*
 * math_get_matrix_determinant.c
 *
 *  Created on: Jul 10, 2009
 *      Author: krang
 */

#include "math_petsc.h"

#include <petscksp.h>

#include <stdio.h>
#include <stdlib.h>
#include <string.h>


double getMatrixDeterminant(
		int num_rows_in,
		int num_columns_in,
		const char *file_name,
		const char *output_file_name)
{
	double a = 0;
	Mat Ap;							// matrix
	PetscErrorCode	ierr;			// error code
	PetscInt i, j,					// iterator variables
			 *col,					// column indices
			 *row,					// row indices
			 num_rows,				// number of rows of matrix
			 num_columns;			// number of columns of matrix
	FILE *file;						// the input file
	IS perm, iperm;
	MatFactorInfo info;
	Vec diagonal;				// the diagonal of the LU factored matrix

	num_rows = num_rows_in;
	num_columns = num_columns_in;

	PetscInt rank, size;

	ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
	ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size); CHKERRQ(ierr);

	if (size != 1)
	{
		SETERRQ(1, "This works only on single processor."); CHKERRQ(ierr);
	}

	printf("Hello from process No. %d of %d\n", rank, size);

	// open file and read number of rows and columns
	ierr = PetscFOpen(PETSC_COMM_SELF, file_name, "r", &file); CHKERRQ(ierr);

	// alocate column and row indices
	col = (PetscInt*) malloc(num_columns * sizeof(PetscInt));
	row = (PetscInt*) malloc(num_rows * sizeof(PetscInt));

	// create matrix
	ierr = MatCreate(PETSC_COMM_SELF, &Ap); CHKERRQ(ierr);
	ierr = MatSetSizes(Ap, PETSC_DECIDE, PETSC_DECIDE, num_rows, num_columns); CHKERRQ(ierr);
	ierr = MatSetFromOptions(Ap); CHKERRQ(ierr);

	// populate row index vector
	for (i = 0; i < num_rows; ++i)
		row[i] = i;
	// populate column index variable
	for (i = 0; i < num_columns; ++i)
		col[i] = i;


	// we read the data from file
	double *row_values = (double*) malloc(num_columns * sizeof(double));
	for (i = 0; i < num_rows; ++i)
	{
		for (j = 0; j < num_columns; ++j)
		{
			fscanf(file, "%lg", &row_values[j]);
		}
		printf("%f %f %f\n", row_values[j-3], row_values[j-2], row_values[j-1]);
		// add row to matrix
		ierr = MatSetValues(Ap, 1, &i, num_columns, col, row_values, INSERT_VALUES); CHKERRQ(ierr);
	}


	// now populate matrix
	ierr = MatAssemblyBegin(Ap, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
	ierr = MatAssemblyEnd(Ap, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);

	// display the matrix
	printf("This is the matrix:\n");
//	ierr = PetscViewerSetFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_COMMON); CHKERRQ(ierr);
	ierr = MatView(Ap, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);

	ierr = MatGetOrdering(Ap,MATORDERING_1WD,&perm,&iperm);CHKERRQ(ierr);
	ierr = MatFactorInfoInitialize(&info);CHKERRQ(ierr);
	info.fill      = 0.0;
	info.dtcol     = 0.0;
	info.shiftnz   = 0.0;
	info.pivotinblocks = 0.0;
	ierr = MatLUFactor(Ap, perm, iperm, &info);CHKERRQ(ierr);
	ierr = MatSetUnfactored(Ap);
	ierr = MatAssemblyBegin(Ap, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
	ierr = MatAssemblyEnd(Ap, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);

	ierr = VecCreateSeq(PETSC_COMM_SELF, num_rows, &diagonal);CHKERRQ(ierr);
	ierr = MatGetDiagonal(Ap, diagonal);CHKERRQ(ierr);
	ierr = VecAssemblyBegin(diagonal); CHKERRQ(ierr);
	ierr = VecAssemblyEnd(diagonal); CHKERRQ(ierr);

	PetscScalar *values = NULL;
	PetscMalloc(sizeof(PetscScalar) * num_rows, &values);
	for(i = 0; i < num_rows; i++) {
	col[i] = i;
	}

	ierr = VecGetValues(diagonal, num_rows, col, values); CHKERRQ(ierr);
	for(i = 0; i < num_rows; i++)
	{
		a = a * values[i];
	}



	// write output to file

	// you have to use petsc viewer to write to file, you *cannot* use fprintf

	// we use Petsc vector to write one element to file
	Vec det;
	ierr = VecCreate(PETSC_COMM_SELF, &det); CHKERRQ(ierr);
	ierr = VecSetSizes(det, PETSC_DECIDE, 1); CHKERRQ(ierr);
	ierr = VecSetFromOptions(det); CHKERRQ(ierr);

	int index = 0;
	// populate vector
	ierr = VecSetValues(det, 1, &index, &a, INSERT_VALUES); CHKERRQ(ierr);
	ierr = VecAssemblyBegin(det); CHKERRQ(ierr);
	ierr = VecAssemblyEnd(det); CHKERRQ(ierr);

	PetscViewer viewer;
	ierr = PetscViewerASCIIOpen(PETSC_COMM_SELF, output_file_name, &viewer); CHKERRQ(ierr);
	ierr = PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_COMMON); CHKERRQ(ierr);
	ierr = VecView(det, viewer); CHKERRQ(ierr);

	// clean-up
	ierr = MatDestroy(Ap); CHKERRQ(ierr);
	ierr = VecDestroy(det); CHKERRQ(ierr);
	ierr = VecDestroy(diagonal); CHKERRQ(ierr);

	free(col);
	free(row);
	free(row_values);

	return a;
}