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

[Yaoyu Hu]

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. It is kind of badly organized and the gmail
web page ignores my 'Tab's, forgive me for that.



Thanks ahead!



============Code Begins=============

/*
 * MI.cpp
 *
 *  Created on: Nov 22, 2011
 *      Author: huyaoyu
 */

static char help[] = "Give the inverse of a matrix.\n\n";

#include <slepceps.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;
 Mat BB,XX;
 IS is_row;
 IS is_col;
 MatFactorInfo mfinfo;
 PetscMPIInt n;
 PetscScalar* array_scalar = new PetscScalar[DIMENSION*DIMENSION];

 for(int i=0;i<DIMENSION*DIMENSION;i++)
 {
  array_scalar[i] = 0.0;
 }
 for(int i=0;i<DIMENSION;i++)
 {
  array_scalar[i*DIMENSION + i] = 1.0;
 }

 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!");
 ierr = PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);CHKERRQ(ierr);

 // B & X
 ierr =
MatCreateSeqDense(PETSC_COMM_WORLD,DIMENSION,DIMENSION,PETSC_NULL,&BB);
CHKERRQ(ierr);
 ierr =
MatCreateSeqDense(PETSC_COMM_WORLD,DIMENSION,DIMENSION,PETSC_NULL,&XX);
CHKERRQ(ierr);

 ierr =
MatSetValues(BB,DIMENSION,idxm,DIMENSION,idxm,array_scalar,INSERT_VALUES);
CHKERRQ(ierr);
 ierr = MatAssemblyBegin(BB,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
 ierr = MatAssemblyEnd(BB,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);

 ierr = ISCreateStride(MPI_COMM_WORLD,DIMENSION,0,1,&is_row); CHKERRQ(ierr);
 ierr = ISCreateStride(MPI_COMM_WORLD,DIMENSION,0,1,&is_col); CHKERRQ(ierr);
 ierr = MatFactorInfoInitialize(&mfinfo); CHKERRQ(ierr);

 // matrix A
 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();

 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);

 ierr = MatDuplicate(A,MAT_COPY_VALUES,&CA); CHKERRQ(ierr);

 ierr = MatLUFactor(A,is_row,is_col,&mfinfo); CHKERRQ(ierr);
 ierr = MatMatSolve(A,BB,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 = MatMatMult(CA,XX,MAT_REUSE_MATRIX,PETSC_DEFAULT,&BB);
 ierr = MatView(BB,PETSC_VIEWER_STDOUT_SELF); CHKERRQ(ierr);

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

 ierr = PetscFinalize();

 delete[] array_scalar;

 return 0;
}

============Code Ends==============
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20111122/39961058/attachment-0001.htm>