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

[Hong Zhang]

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. 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==============