[petsc-users] Inverse matrix empty

Mon Aug 9 14:14:33 CDT 2010

Here is my code,

Thank you,
regards,
Nemanja

#include "math_petsc.h"

#include <petscksp.h>

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

double** getInverseMatrix(
		int num_rows_in,
		int num_columns_in,
		const char *file_name,
		const char *output_file_name)
{
	double **X = NULL;
	Mat Ap;							// input matrix
	Mat Bp,							// identity matrix
		Xp;							// result matrix
	PetscErrorCode	ierr;			// error code
	PetscInt i,	j,					// iterator variable
			 *col,				// column indices
			 *row,				// row indices
			 ap_num_rows,					// number of rows
			 ap_num_columns;					// number of columns
	FILE *file;					// input file
	IS perm_rows, perm_cols;
	MatFactorInfo lu_info, fact_info;
	KSP ksp;
	PC pc;
	IS perm, iperm;
	MatFactorInfo info;

	PetscInt rank, size;
	ap_num_rows = num_rows_in;
	ap_num_columns = num_columns_in;

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

	printf("MPI_Comm_rank = %d\nMPI_Comm_size = %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(ap_num_columns * sizeof(PetscInt));
	row = (PetscInt*) malloc(ap_num_rows * sizeof(PetscInt));

	// create matrices
	ierr = MatCreate(PETSC_COMM_WORLD, &Ap); CHKERRQ(ierr);
	ierr = MatSetSizes(Ap, PETSC_DECIDE, PETSC_DECIDE, ap_num_rows, ap_num_columns); CHKERRQ(ierr);
	ierr = MatSetFromOptions(Ap); CHKERRQ(ierr);

	ierr = MatCreate(PETSC_COMM_WORLD, &Bp); CHKERRQ(ierr);
	ierr = MatSetSizes(Bp, PETSC_DECIDE, PETSC_DECIDE, ap_num_rows, ap_num_columns); CHKERRQ(ierr);
	ierr = MatSetFromOptions(Bp); CHKERRQ(ierr);

	ierr = MatCreate(PETSC_COMM_WORLD, &Xp); CHKERRQ(ierr);
	ierr = MatSetSizes(Xp, PETSC_DECIDE, PETSC_DECIDE, ap_num_rows, ap_num_columns); CHKERRQ(ierr);
	ierr = MatSetFromOptions(Xp); CHKERRQ(ierr);

/*	ierr = MatCreate(PETSC_COMM_WORLD, &fact); CHKERRQ(ierr);
	ierr = MatSetSizes(fact, PETSC_DECIDE, PETSC_DECIDE, m, n); CHKERRQ(ierr);
	ierr = MatSetFromOptions(fact); CHKERRQ(ierr);
*/

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

	// now populate matrix
	double *row_values = (double*) malloc(ap_num_columns * sizeof(double));
	for (i = 0; i < ap_num_rows; ++i)
	{
		for (j = 0; j < ap_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, ap_num_columns, col, row_values, INSERT_VALUES); CHKERRQ(ierr);
	}
	ierr = MatAssemblyBegin(Ap, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
	ierr = MatAssemblyEnd(Ap, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);

	ierr = MatView(Ap, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);

	// factor the matrix
	ierr = MatGetOrdering(Ap, MATORDERING_1WD, &perm, &iperm);CHKERRQ(ierr);
	ierr = MatFactorInfoInitialize(&info);CHKERRQ(ierr);
	ierr = MatLUFactor(Ap, perm, iperm, &info);CHKERRQ(ierr);

	// the identity matrix
	double identity_values[ap_num_columns];
	memset(identity_values, 0, ap_num_columns * sizeof(double));
	for (i = 0; i < ap_num_rows; ++i)
	{
		identity_values[i] = 1;
		if (i != 0)
			identity_values[i-1] = 0;
		ierr = MatSetValues(Bp, 1, &i, ap_num_columns, col, identity_values, INSERT_VALUES); CHKERRQ(ierr);
	}
	ierr = MatAssemblyBegin(Bp, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
	ierr = MatAssemblyEnd(Bp, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);

	// view identity matrix
	printf("This is the identity matrix:\n");
	ierr = MatView(Bp, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);

	// create solver context
	ierr = KSPCreate(PETSC_COMM_WORLD, &ksp); CHKERRQ(ierr);
	ierr = KSPSetOperators(ksp, Ap, Ap, DIFFERENT_NONZERO_PATTERN); CHKERRQ(ierr);
//	ierr = KSPSetOperators(ksp, Bp, Bp, DIFFERENT_NONZERO_PATTERN); CHKERRQ(ierr);
//	ierr = KSPSetOperators(ksp, Xp, Xp, DIFFERENT_NONZERO_PATTERN); CHKERRQ(ierr);
	ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
	ierr = PCSetType(pc,PCJACOBI);CHKERRQ(ierr);
	ierr = KSPSetTolerances(ksp,1.e-7,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr);
	ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);

	// !!!!!		get inverse matrix		!!!!!
	ierr = MatMatSolve(Ap, Bp, Xp); CHKERRQ(ierr);

	// assemble matrix
	ierr = MatAssemblyBegin(Xp, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
	ierr = MatAssemblyEnd(Xp, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);

	// print result
	printf("This is the result:\n");
	ierr = MatView(Xp, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);

	// write result to file

	PetscViewer viewer;
	ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD, output_file_name, &viewer); CHKERRQ(ierr);
	ierr = PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_COMMON); CHKERRQ(ierr);
	ierr = MatView(Xp, viewer); CHKERRQ(ierr);

	// clean-up
	ierr = MatDestroy(Ap); CHKERRQ(ierr);
	ierr = MatDestroy(Bp); CHKERRQ(ierr);
	ierr = MatDestroy(Xp); CHKERRQ(ierr);

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

	return X;
}

On Monday 09 August 2010 20:14:01 you wrote:
> Send me the segment of your code for computing inv(A).
> I'll check it.
> 
> Hong
> 
> On Mon, Aug 9, 2010 at 1:00 PM, Немања Илић (Nemanja Ilic)
> <nemanja.ilic.81 at gmail.com> wrote:
> > Hello,
> >
> > I want to get the inverse of a matrix. I use the following algorithm: I have input matrix (A), identity matrix (B) and I try to solve A * X = B, where X is the inverse matrix. I use MatMatSolve function. I included blas, fblas, and superlu libraries in my project but still I get empty inverse matrix. I factored A with LUFactor before calling MatMatSolve. What am I doing wrong?
> >
> > Thank you in advance,
> > Regards,
> > Nemanja
> >
> 

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20100809/3e3768de/attachment.htm>