[petsc-users] 답장: Re: The inverse of seqaij matrix
Joon Hee Choi
choi240 at purdue.edu
Fri Sep 27 08:50:52 CDT 2013
Dear Brry,
Thank you for your fast reply. I generated a SEQAIJ matrix in MPI rank 1 from a MPIAIJ matrix, because I needed the inverse matrix of the MPI matrix.
Actually, the furmula is very simple, but my code is very complex.
What I want to get is the inverse of (temp1 * temp2). temp1 and temp2 are mpiaij and * is hadamard product. Also, I want to get mpiaij matrix as the result. I think the ordering used to factor the matrix was not good because seqM2 is not singular. From my code, could you let me know why the ordering is wrong? I am sorry, but I am attaching my very long code. Thank you again.
Joon
MPI_Comm_size(PETSC_COMM_WORLD, &size);
MPI_Comm_rank(MPI_COMM_SELF, &myrank);
ierr = MatCreate(PETSC_COMM_WORLD, &invM2); CHKERRQ(ierr);
ierr = MatSetType(invM2, MATDENSE); CHKERRQ(ierr);
ierr = MatSetSizes(invM2, PETSC_DECIDE, PETSC_DECIDE, R, R); CHKERRQ(ierr);
ierr = MatSeqDenseSetPreallocation(invM2, NULL); CHKERRQ(ierr);
ierr = MatMPIDenseSetPreallocation(invM2, NULL); CHKERRQ(ierr);
ierr = MatGetOwnershipRange(temp1, &rstart, &rend); CHKERRQ(ierr);
ierr = PetscMalloc((rend-rstart)*R*sizeof(PetscScalar), &cell); CHKERRQ(ierr);
for (i=rstart; i<rend; i++) {
for (j=0; j<R; j++) {
ierr = MatGetValues(temp1, 1, &i, 1, &j, &cell1); CHKERRQ(ierr);
ierr = MatGetValues(temp2, 1, &i, 1, &j, &cell2); CHKERRQ(ierr);
cell[(i-rstart)*R+j] = cell1 * cell2;
}
}
rbuf = (PetscScalar *)malloc(R*R*sizeof(PetscScalar));
PetscInt *idx;
idx = (PetscInt *)malloc(R*sizeof(PetscInt));
for (i=0; i<R; i++) idx[i]=i;
MPI_Gather(cell, (rend-rstart)*R, MPIU_SCALAR, rbuf, (rend-rstart)*R, MPIU_SCALAR, 0, MPI_COMM_WORLD);
if (myrank == 0) {
ierr = MatCreateSeqAIJ(PETSC_COMM_SELF, R, R, R, NULL, &invM2); CHKERRQ(ierr);
ierr = MatSetValues(invM2, R, idx, R, idx, rbuf, INSERT_VALUES); CHKERRQ(ierr);
ierr = MatAssemblyBegin(invM2, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
ierr = MatAssemblyEnd(invM2, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
ierr = MatTranspose(invM2, MAT_INITIAL_MATRIX, &invM2T); CHKERRQ(ierr);
ierr = MatCreateSeqDense(PETSC_COMM_SELF, R, R, NULL, &I); CHKERRQ(ierr);
ierr = MatZeroEntries(I); CHKERRQ(ierr);
ierr = MatShift(I, 1.0); CHKERRQ(ierr);
ierr = MatCreateSeqDense(PETSC_COMM_SELF, R, R, NULL, &seqM2T); CHKERRQ(ierr);
ierr = MatGetOrdering(invM2T, MATORDERINGNATURAL, &isr, &isc); CHKERRQ(ierr);
ierr = MatFactorInfoInitialize(&info); CHKERRQ(ierr);
ierr = MatLUFactor(invM2T, isr, isc, &info); CHKERRQ(ierr);
PetscPrintf(PETSC_COMM_WORLD, "invM2T \n");
MatView(invM2T, PETSC_VIEWER_STDOUT_WORLD);
ierr = MatMatSolve(invM2T, I, seqM2T); CHKERRQ(ierr);
PetscPrintf(PETSC_COMM_WORLD, "seqM2T \n");
MatView(seqM2T, PETSC_VIEWER_STDOUT_WORLD);
ierr = MatDenseGetArray(seqM2T, &rbuf); CHKERRQ(ierr);;
}
MPI_Scatter(rbuf, (rend-rstart)*R, MPIU_SCALAR, cell, (rend-rstart)*R, MPIU_SCALAR, 0, MPI_COMM_WORLD);
M = R;
clocm = PETSC_DECIDE;
ierr = PetscSplitOwnership(PETSC_COMM_WORLD, &clocm, &M); CHKERRQ(ierr);
ierr = MPI_Scan(&clocm, &cend, size, MPI_INT, MPI_SUM, PETSC_COMM_WORLD); CHKERRQ(ierr);
cstart = cend - clocm;
ierr = MatCreate(PETSC_COMM_WORLD, M2T); CHKERRQ(ierr);
ierr = MatSetType(*M2T, MATAIJ); CHKERRQ(ierr);
ierr = MatSetSizes(*M2T, PETSC_DECIDE, PETSC_DECIDE, R, R); CHKERRQ(ierr);
ierr = MatSeqAIJSetPreallocation(*M2T, R, NULL); CHKERRQ(ierr);
ierr = MatMPIAIJSetPreallocation(*M2T, clocm, NULL, R-clocm, NULL); CHKERRQ(ierr);
for (i=rstart; i<rend; i++) {
for (j=0; j<R; j++) {
ierr = MatSetValue(*M2T, i, j, cell[(i-rstart)*R+j], INSERT_VALUES); CHKERRQ(ierr);
}
}
ierr = MatAssemblyBegin(*M2T, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
ierr = MatAssemblyEnd(*M2T, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
}
----- 원본 메시지 -----
보낸 사람: Barry Smith <bsmith at mcs.anl.gov>
받는 사람: Joon Hee Choi <choi240 at purdue.edu>
참조: petsc-users at mcs.anl.gov
보냄: Fri, 27 Sep 2013 09:05:58 -0400 (EDT)
제목: Re: [petsc-users] The inverse of seqaij matrix
> Based on the output below it appears that the code block you sent was completely successfully since it printed all three matrices you tried to proint. So it is something else in your code that is generating the zero pivot.
> BTW: it should not be
>MatView(invM2T, PETSC_VIEWER_STDOUT_WORLD);
>it should be
>MatView(seqM2T, PETSC_VIEWER_STDOUT_WORLD);
>to display the inverse
>And how is it possible that the error occurred on MPI rank 1 but the matrices are sequential matrices viewed on PETSC_VIEWER_STDOUT_WORLD?
>A zero privet means that either the matrix is singular or the ordering used to factor the matrix is not a good ordering. Perhaps the matrix was generated incorrectly.
>The code crashed because likely your code returned an error code returned from PETSc and tried to keep running. For example if the factorization generated an error code but you called the MatSolveSolve on the result anyway, you cannot do that.
On Sep 27, 2013, at 7:21 AM, Joon Hee Choi <choi240 at purdue.edu> wrote:
>> Hi.
>> I am trying to get the inverse matrix of seqaij. I followed FAQ of petsc manual, but I got error.
>> Could someone tell me what I am missing? I am attaching my code, results, and error message.
>> Thank you.
>>
>> Joon
>>
>>
>>
>> PetscPrintf(PETSC_COMM_WORLD, "I \n");
>> MatView(I, PETSC_VIEWER_STDOUT_WORLD);
>> PetscPrintf(PETSC_COMM_WORLD, "invM2T \n");
>> MatView(invM2T, PETSC_VIEWER_STDOUT_WORLD);
>>
>> ierr = MatCreateSeqDense(PETSC_COMM_SELF, R, R, NULL, &seqM2T); CHKERRQ(ierr);
>> ierr = MatGetOrdering(invM2T, MATORDERINGNATURAL, &isr, &isc); CHKERRQ(ierr);
>> ierr = MatFactorInfoInitialize(&info); CHKERRQ(ierr);
>> ierr = MatLUFactor(invM2T, isr, isc, &info); CHKERRQ(ierr);
>>
>> ierr = MatMatSolve(invM2T, I, seqM2T); CHKERRQ(ierr);
>> PetscPrintf(PETSC_COMM_WORLD, "invM2T \n");
>> MatView(invM2T, PETSC_VIEWER_STDOUT_WORLD);
>>
>>
>> Matrix Object: 1 MPI processes
>> type: seqdense
>> 1.0000000000000000e+00 0.0000000000000000e+00
>> 0.0000000000000000e+00 1.0000000000000000e+00
>> invM2T
>> Matrix Object: 1 MPI processes
>> type: seqaij
>> row 0: (0, 3.22689) (1, 2.47603)
>> row 1: (0, 2.47603) (1, 2.25324)
>> invM2T
>> Matrix Object: 1 MPI processes
>> type: seqaij
>> row 0: (0, 3.22689) (1, 2.47603)
>> row 1: (0, 0.76731) (1, 0.353361)
>>
>>
>> [1]PETSC ERROR: --------------------- Error Message ------------------------------------
>> [1]PETSC ERROR: Detected zero pivot in LU factorization:
>> see http://www.mcs.anl.gov/petsc/documentation/faq.html#ZeroPivot!
>> [1]PETSC ERROR: Zero pivot row 0 value 0 tolerance 0!
>> [1]PETSC ERROR: ------------------------------------------------------------------------
>> [1]PETSC ERROR: Petsc Release Version 3.4.1, Jun, 10, 2013
>> [1]PETSC ERROR: See docs/changes/index.html for recent updates.
>> [1]PETSC ERROR: See docs/faq.html for hints about trouble shooting.
>> [1]PETSC ERROR: See docs/index.html for manual pages.
>> [1]PETSC ERROR: ------------------------------------------------------------------------
>> [1]PETSC ERROR: ./tenfac on a linux named rossmann-b000.rcac.purdue.edu by choi240 Fri Sep 27 08:01:31 2013
>> [1]PETSC ERROR: Libraries linked from /apps/rhel5/petsc-3.4.1/64/impi-4.1.0.030_intel-13.0.1.117_ind64/linux/lib
>> [1]PETSC ERROR: Configure run at Sun Jun 23 10:20:51 2013
>> [1]PETSC ERROR: Configure options --with-cc=mpiicc --with-cxx=mpiicpc --with-fc=mpiifort --download-sowing --with-scalar-type=real --with-shared-libraries=1 --with-pic=1 --with-clanguage=C++ --with-fortran --with-fortran-kernels=1 --with-64-bit-indices=1 --with-debugging=0 --with-blas-lapack-dir=/opt/intel/composer_xe_2013.1.117/mkl/lib/intel64 --COPTFLAGS=-O3 --CXXOPTFLAGS=-O3 --FOPTFLAGS=-O3 --download-hdf5=no --download-metis=no --download-parmetis=no --download-superlu_dist=no --download-mumps=no --download-scalapack=yes --download-blacs=no --download-hypre=no --download-spooles=no
>> [1]PETSC ERROR: ------------------------------------------------------------------------
>> [1]PETSC ERROR: MatPivotCheck_none() line 589 in /apps/rhel5/petsc-3.4.1/64/impi-4.1.0.030_intel-13.0.1.117_ind64/include/petsc-private/matimpl.h
>> [1]PETSC ERROR: MatPivotCheck() line 608 in /apps/rhel5/petsc-3.4.1/64/impi-4.1.0.030_intel-13.0.1.117_ind64/include/petsc-private/matimpl.h
>> [1]PETSC ERROR: MatLUFactorNumeric_SeqAIJ_Inode() line 1422 in /apps/rhel5/petsc-3.4.1/64/impi-4.1.0.030_intel-13.0.1.117_ind64/src/mat/impls/aij/seq/inode.c
>> [1]PETSC ERROR: MatLUFactorNumeric() line 2889 in /apps/rhel5/petsc-3.4.1/64/impi-4.1.0.030_intel-13.0.1.117_ind64/src/mat/interface/matrix.c
>> [1]PETSC ERROR: MatLUFactor_SeqAIJ() line 971 in /apps/rhel5/petsc-3.4.1/64/impi-4.1.0.030_intel-13.0.1.117_ind64/src/mat/impls/aij/seq/aijfact.c
>> [1]PETSC ERROR: MatLUFactor() line 2716 in /apps/rhel5/petsc-3.4.1/64/impi-4.1.0.030_intel-13.0.1.117_ind64/src/mat/interface/matrix.c
>> [1]PETSC ERROR: cal_M2() line 83 in tenfac.cpp
>> [1]PETSC ERROR: ------------------------------------------------------------------------
>> [1]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range
>> [1]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
>> [1]PETSC ERROR: or see http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind[1]PETSC ERROR: or try http://valgrind.org on GNU/linux and Apple Mac OS X to find memory corruption errors
>> [1]PETSC ERROR: configure using --with-debugging=yes, recompile, link, and run
>> [1]PETSC ERROR: to get more information on the crash.
>> [1]PETSC ERROR: --------------------- Error Message ------------------------------------
>> [1]PETSC ERROR: Signal received!
>> [1]PETSC ERROR: ------------------------------------------------------------------------
>> [1]PETSC ERROR: Petsc Release Version 3.4.1, Jun, 10, 2013
>> [1]PETSC ERROR: See docs/changes/index.html for recent updates.
>> [1]PETSC ERROR: See docs/faq.html for hints about trouble shooting.
>> [1]PETSC ERROR: See docs/index.html for manual pages.
>> [1]PETSC ERROR: ------------------------------------------------------------------------
>> [1]PETSC ERROR: ./tenfac on a linux named rossmann-b000.rcac.purdue.edu by choi240 Fri Sep 27 08:01:31 2013
>> [1]PETSC ERROR: Libraries linked from /apps/rhel5/petsc-3.4.1/64/impi-4.1.0.030_intel-13.0.1.117_ind64/linux/lib
>> [1]PETSC ERROR: Configure run at Sun Jun 23 10:20:51 2013
>> [1]PETSC ERROR: Configure options --with-cc=mpiicc --with-cxx=mpiicpc --with-fc=mpiifort --download-sowing --with-scalar-type=real --with-shared-libraries=1 --with-pic=1 --with-clanguage=C++ --with-fortran --with-fortran-kernels=1 --with-64-bit-indices=1 --with-debugging=0 --with-blas-lapack-dir=/opt/intel/composer_xe_2013.1.117/mkl/lib/intel64 --COPTFLAGS=-O3 --CXXOPTFLAGS=-O3 --FOPTFLAGS=-O3 --download-hdf5=no --download-metis=no --download-parmetis=no --download-superlu_dist=no --download-mumps=no --download-scalapack=yes --download-blacs=no --download-hypre=no --download-spooles=no
>> [1]PETSC ERROR: ------------------------------------------------------------------------
>> [1]PETSC ERROR: User provided function() line 0 in unknown directory unknown file
>> application called MPI_Abort(MPI_COMM_WORLD, 59) - process 1
More information about the petsc-users
mailing list