[petsc-users] Only one column is in the matrix after solving the inverse matrix
Hong Zhang
hzhang at mcs.anl.gov
Sun Nov 27 20:17:51 CST 2011
huyaoyu :
This is a bug in petsc library. I've fixed it in petsc-3.2 and patched
petsc-dev.
You may edit ~petsc-3.2/src/mat/impls/dense/seq/dense.c by following
http://petsc.cs.iit.edu/petsc/petsc-dev/rev/c2fe40bf559c
or get an updated petsc-3.2, then rebuild your petsc-3.2.
Thanks for reporting the problem!
Hong
> 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
>
>
>
>
More information about the petsc-users
mailing list