[petsc-users] Only one column is in the matrix after solving the inverse matrix

[huyaoyu]

Hong,

Thank you for your help!

I tried the example code of petsc-3.2/src/mat/examples/tests/ex1.c.
But the result is the same. I attach the modified code and the
results here. I don't know what is the problem exactly. By the way I am
using PETSc on a 64bit ubuntu 11.04. I complied MPICH2 myself, made
symbolic links in the /usr/bin for mpicc, mpiexec etc, and configured
PETSc use the command line:
./configure --with-cc=mpicc --with-fc=mpif90 --download-f-blas-lapack=1

I break some lines of the code to avoid line wrapping in Evolution mail
program.

> You can take a look at examples
> petsc-3.2/src/mat/examples/tests/ex1.c, ex125.c and ex125.c.
> 
> Hong
> 
> On Tue, Nov 22, 2011 at 9:23 AM, Yaoyu Hu <huyaoyu1986 at gmail.com> wrote:
> > Hi, everyone,
> >
> > I am new to PETSc, and I have just begun to use it together with slepc. It
> > is really fantastic and I like it!
> >
> >
> >
> > It?was not from me but one of my colleges who wanted to solve an inverse of
> > a matrix, which had 400 rows and columns.
> >
> >
> >
> > I know that it is not good to design a algorithm that has a process for
> > solving the inverse of a matrix. I just wanted to?give it a try. However,
> > things turned out wired. I followed the instructions on the FAQ web page of
> > PETSc, the one using MatLUFractor() and MatMatSolve(). After I finished the
> > coding, I tried the program. The result I got?was a matrix which only has
> > its first column but nothing else. I did not know what's happened. The
> > following is the codes I used. 
> >
=============Modified code begins=================
static char help[] = "Give the inverse of a matrix.\n\n";

#include <petscmat.h>

#include <iostream>
#include <fstream>

#undef __FUNCT__
#define __FUNCT__ "main"

#define DIMENSION 2

int main(int argc,char **args)
{
	PetscErrorCode	ierr;
	PetscMPIInt	size;
	Mat		A,CA;	// CA is the copy of A
	Mat		RHS,XX;	// XX is the inverse result
	MatFactorInfo	mfinfo;
	PetscScalar*	array_scalar = 
		new PetscScalar[DIMENSION*DIMENSION];
	PetscScalar*	array_for_RHS;

	// clean the values of array_scalar
	for(int i=0;i<DIMENSION*DIMENSION;i++)
	{
		array_scalar[i] = 0.0;
	}

	// set up the indices
	int idxm[DIMENSION];
	for(int i=0;i<DIMENSION;i++)
	{
		idxm[i] = i;
	}

	PetscInitialize(&argc,&args,(char *)0,help);
	ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
	if (size != 1) SETERRQ(PETSC_COMM_WORLD,
		1,"This is a uniprocessor example only!");

	// RHS & XX
	ierr = MatCreate(PETSC_COMM_WORLD,&RHS);CHKERRQ(ierr);
	ierr =MatSetSizes(RHS,
		PETSC_DECIDE,PETSC_DECIDE,
		DIMENSION,DIMENSION);CHKERRQ(ierr);
	ierr = MatSetType(RHS,MATDENSE);CHKERRQ(ierr);
	ierr = MatSetFromOptions(RHS);CHKERRQ(ierr);
	ierr = MatGetArray(RHS,&array_for_RHS); CHKERRQ(ierr);
	for(int i=0;i<DIMENSION;i++)
	{
			array_for_RHS[i*DIMENSION + i] = 1.0;
	}
	ierr = MatRestoreArray(RHS,&array_for_RHS); CHKERRQ(ierr);

	ierr = MatDuplicate(RHS,
		MAT_DO_NOT_COPY_VALUES,&XX);CHKERRQ(ierr);

	/* matrix A */
	/* read in the file A.txt. It is a single process program so I think it
is OK to read a file like this.*/
	std::fstream in_file;
	double temp_scalar;
	in_file.open("./A.txt",std::ifstream::in);
	if(!in_file.good())
	{
		ierr = PetscPrintf(PETSC_COMM_SELF,
			"File open failed!\n"); CHKERRQ(ierr);
		return 1;
	}

	for(int i=0;i<DIMENSION;i++)
	{
		for(int j=0;j<DIMENSION;j++)
		{
			in_file>>temp_scalar;

			array_scalar[i*DIMENSION + j] = temp_scalar;
		}
	}
	in_file.close();

	// matrices creation and initialization
	ierr = MatCreateSeqDense(PETSC_COMM_WORLD,
		DIMENSION,DIMENSION,PETSC_NULL,&A); CHKERRQ(ierr);
	ierr = MatCreateSeqDense(PETSC_COMM_WORLD,
		DIMENSION,DIMENSION,PETSC_NULL,&CA); CHKERRQ(ierr);
	ierr = MatSetValues(A,DIMENSION,idxm,DIMENSION,idxm,
			array_scalar,INSERT_VALUES); CHKERRQ(ierr);

	ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
	ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
	// create CA
	ierr = MatDuplicate(A,MAT_COPY_VALUES,&CA); CHKERRQ(ierr);

	ierr = PetscPrintf(PETSC_COMM_SELF,
		"The A matrix is:\n"); CHKERRQ(ierr);
	ierr = MatView(A,PETSC_VIEWER_STDOUT_SELF); CHKERRQ(ierr);

	// in-place LUFactor
	ierr = MatLUFactor(A,0,0,0); CHKERRQ(ierr);
	// solve for the inverse matrix XX
	ierr = MatMatSolve(A,RHS,XX); CHKERRQ(ierr);

	ierr = PetscPrintf(PETSC_COMM_SELF,
		"The inverse of A matrix is:\n"); CHKERRQ(ierr);
	ierr = MatView(XX,PETSC_VIEWER_STDOUT_SELF); CHKERRQ(ierr);

	ierr = PetscPrintf(PETSC_COMM_SELF,
		"The multiplied result is:\n"); CHKERRQ(ierr);
	ierr = MatMatMult(CA,XX,MAT_REUSE_MATRIX,PETSC_DEFAULT,&RHS);
	ierr = MatView(RHS,PETSC_VIEWER_STDOUT_SELF); CHKERRQ(ierr);

	// destroy
	ierr = MatDestroy(&A); CHKERRQ(ierr);
	ierr = MatDestroy(&RHS); CHKERRQ(ierr);
	ierr = MatDestroy(&XX); CHKERRQ(ierr);
	ierr = MatDestroy(&CA); CHKERRQ(ierr);

	ierr = PetscFinalize();

	delete[] array_scalar;

	return 0;
}
===========Modified code ends==============
The input file(DIMENSION = 2) and the results are
A.txt:
2.0 3.0
-20.0 55.0
Results:
The A matrix is:
Matrix Object: 1 MPI processes
  type: seqdense
2.0000000000000000e+00 3.0000000000000000e+00 
-2.0000000000000000e+01 5.5000000000000000e+01 
The inverse of A matrix is:
Matrix Object: 1 MPI processes
  type: seqdense
3.2352941176470590e-01 -0.0000000000000000e+00 
1.1764705882352941e-01 0.0000000000000000e+00 
The multiplied result is:
Matrix Object: 1 MPI processes
  type: seqdense
1.0000000000000000e+00 0.0000000000000000e+00 
0.0000000000000000e+00 0.0000000000000000e+00