Direct LU solver
Hong Zhang
hzhang at mcs.anl.gov
Tue Mar 4 11:51:31 CST 2008
Amit,
I'll test it and let you know the result.
Hong
On Tue, 4 Mar 2008, Amit.Itagi at seagate.com wrote:
>
>
>
> Here is an example. I am creating the matrices, the vectors and the solver.
> I can do an LU solve on two processes just fine. However, when I destroy
> the entities and repeat the process (without any change), I get a crash.
>
> #include <iostream>
> #include <complex>
> #include <stdlib.h>
> #include "petsc.h"
> #include "petscmat.h"
> #include "petscvec.h"
> #include "petscksp.h"
>
> using namespace std;
>
> int main( int argc, char *argv[] ) {
>
> int rank, size;
> Mat A;
> PetscErrorCode ierr;
> PetscInt loc;
> PetscScalar val;
> Vec x, y;
> KSP solver;
> PC prec;
> MPI_Comm comm;
>
> // Number of non-zeros in each row
>
> ierr=PetscInitialize(&argc,&argv,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
> // Initialization (including MPI)
>
> comm=PETSC_COMM_WORLD;
>
> ierr=MPI_Comm_size(PETSC_COMM_WORLD,&size); CHKERRQ(ierr);
> ierr=MPI_Comm_rank(PETSC_COMM_WORLD,&rank); CHKERRQ(ierr);
>
>
> //================== LU : First Time =============================
>
> // Assemble matrix A
>
> PetscInt d_nnz=1, o_nnz=1;
> ierr=MatCreate(comm,&A); CHKERRQ(ierr);
> ierr=MatSetSizes(A,1,1,2,2); CHKERRQ(ierr);
> ierr=MatMPIAIJSetPreallocation(A,0,&d_nnz,0,&o_nnz); CHKERRQ(ierr);
> ierr=MatSetType(A,MATSUPERLU_DIST); CHKERRQ(ierr);
> ierr=MatSetFromOptions(A); CHKERRQ(ierr);
>
>
> if(rank==0) {
> val=complex<double>(1.0,0.0);
> ierr=MatSetValue(A,0,0,val,ADD_VALUES);CHKERRQ(ierr);
> val=complex<double>(0.0,1.0);
> ierr=MatSetValue(A,0,1,val,ADD_VALUES);CHKERRQ(ierr);
> }
> else {
> val=complex<double>(1.0,1.0);
> ierr=MatSetValue(A,1,1,val,ADD_VALUES);CHKERRQ(ierr);
> val=complex<double>(0.0,1.0);
> ierr=MatSetValue(A,1,0,val,ADD_VALUES);CHKERRQ(ierr);
> }
>
> ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
> ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
>
> cout << "============ Mat A ==================" << endl;
> ierr=MatView(A,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
> cout << "======================================" << endl;
>
> // Direct LU solver
> ierr=KSPCreate(comm,&solver); CHKERRQ(ierr);
> ierr=KSPSetOperators(solver,A,A,SAME_PRECONDITIONER); CHKERRQ(ierr);
> ierr=KSPSetType(solver,KSPPREONLY); CHKERRQ(ierr);
> ierr=KSPGetPC(solver,&prec); CHKERRQ(ierr);
> ierr=PCSetType(prec,PCLU); CHKERRQ(ierr);
> ierr=PCFactorSetShiftNonzero(prec,PETSC_DECIDE);
> ierr=KSPSetFromOptions(solver); CHKERRQ(ierr);
>
> //============ Vector assembly ========================
>
> if(rank==0) {
> ierr=VecCreateMPI(PETSC_COMM_WORLD,1,2,&x); CHKERRQ(ierr);
> val=complex<double>(1.0,0.0);
> loc=0;
> ierr=VecSetValues(x,1,&loc,&val,ADD_VALUES);CHKERRQ(ierr);
> }
> else {
> ierr=VecCreateMPI(PETSC_COMM_WORLD,1,2,&x); CHKERRQ(ierr);
> val=complex<double>(-1.0,0.0);
> loc=1;
> ierr=VecSetValues(x,1,&loc,&val,ADD_VALUES);CHKERRQ(ierr);
> }
>
> ierr=VecAssemblyBegin(x); CHKERRQ(ierr);
> ierr=VecAssemblyEnd(x); CHKERRQ(ierr);
>
> cout << "============== Vec x ==================" << endl;
> ierr=VecView(x,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
> cout << "======================================" << endl;
>
> VecDuplicate(x,&y); // Duplicate the matrix storage
>
> // Solve the matrix equation
> ierr=KSPSolve(solver,x,y); CHKERRQ(ierr);
>
> cout << "============== Vec y =================" << endl;
> ierr=VecView(y,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
> cout << "======================================" << endl;
>
>
> // Destructors
> ierr=KSPDestroy(solver); CHKERRQ(ierr);
> ierr=VecDestroy(x); CHKERRQ(ierr);
> ierr=VecDestroy(y); CHKERRQ(ierr);
> ierr=MatDestroy(A); CHKERRQ(ierr);
>
>
> //============== LU Second time ==============================
> //======= Everything is exactly the same =====================
>
> ierr=MatCreate(comm,&A); CHKERRQ(ierr);
> ierr=MatSetSizes(A,1,1,2,2); CHKERRQ(ierr);
> ierr=MatMPIAIJSetPreallocation(A,0,&d_nnz,0,&o_nnz); CHKERRQ(ierr);
> ierr=MatSetType(A,MATSUPERLU_DIST); CHKERRQ(ierr);
> ierr=MatSetFromOptions(A); CHKERRQ(ierr);
>
>
> if(rank==0) {
> val=complex<double>(1.0,0.0);
> ierr=MatSetValue(A,0,0,val,ADD_VALUES);CHKERRQ(ierr);
> val=complex<double>(0.0,1.0);
> ierr=MatSetValue(A,0,1,val,ADD_VALUES);CHKERRQ(ierr);
> }
> else {
> val=complex<double>(1.0,1.0);
> ierr=MatSetValue(A,1,1,val,ADD_VALUES);CHKERRQ(ierr);
> val=complex<double>(0.0,1.0);
> ierr=MatSetValue(A,1,0,val,ADD_VALUES);CHKERRQ(ierr);
> }
>
> ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
> ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
>
> cout << "============ Mat A ==================" << endl;
> ierr=MatView(A,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
> cout << "======================================" << endl;
>
> // Direct LU solver
> ierr=KSPCreate(comm,&solver); CHKERRQ(ierr);
> ierr=KSPSetOperators(solver,A,A,SAME_PRECONDITIONER); CHKERRQ(ierr);
> ierr=KSPSetType(solver,KSPPREONLY); CHKERRQ(ierr);
> ierr=KSPGetPC(solver,&prec); CHKERRQ(ierr);
> ierr=PCSetType(prec,PCLU); CHKERRQ(ierr);
> ierr=PCFactorSetShiftNonzero(prec,PETSC_DECIDE);
> ierr=KSPSetFromOptions(solver); CHKERRQ(ierr);
>
> //============ Vector assembly ========================
>
> if(rank==0) {
> ierr=VecCreateMPI(PETSC_COMM_WORLD,1,2,&x); CHKERRQ(ierr);
> val=complex<double>(1.0,0.0);
> loc=0;
> ierr=VecSetValues(x,1,&loc,&val,ADD_VALUES);CHKERRQ(ierr);
> }
> else {
> ierr=VecCreateMPI(PETSC_COMM_WORLD,1,2,&x); CHKERRQ(ierr);
> val=complex<double>(-1.0,0.0);
> loc=1;
> ierr=VecSetValues(x,1,&loc,&val,ADD_VALUES);CHKERRQ(ierr);
> }
>
> ierr=VecAssemblyBegin(x); CHKERRQ(ierr);
> ierr=VecAssemblyEnd(x); CHKERRQ(ierr);
>
> cout << "============== Vec x ==================" << endl;
> ierr=VecView(x,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
> cout << "======================================" << endl;
>
> VecDuplicate(x,&y); // Duplicate the matrix storage
>
> // Solve the matrix equation
> ierr=KSPSolve(solver,x,y); CHKERRQ(ierr);
>
> cout << "============== Vec y =================" << endl;
> ierr=VecView(y,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
> cout << "======================================" << endl;
>
>
> // Destructors
> ierr=KSPDestroy(solver); CHKERRQ(ierr);
> ierr=VecDestroy(x); CHKERRQ(ierr);
> ierr=VecDestroy(y); CHKERRQ(ierr);
> ierr=MatDestroy(A); CHKERRQ(ierr);
>
>
> // Finalize
> ierr=PetscFinalize(); CHKERRQ(ierr);
>
>
> return 0;
>
> }
>
> Thanks
>
> Rgds,
> Amit
>
> Hong Zhang <hzhang at mcs.anl.gov> wrote on 03/03/2008 09:55:43 PM:
>
>>
>> Are you using petsc-dev or a released version?
>> If petsc-dev with the newly released superlu_dist_2.2,
>> there is a bug in reusing parallel LU symbolic factorization
>> which I've reported to the developer of superlu_dist_2.2.
>>
>>>
>>> Your comments were helpful. Now I can do the LU solves multiple times.
>>> However within the program, if I destroy the matrix and the KSP and
> create
>>> them again (for a different matrix), I get an error in superlu.
>>
>> Otherwise,
>> can you send me a simplified code that display the crash?
>> Then I may run it and see if petsc superlu_dist interface has bug.
>> You may also test other parallel direct LU solvers, e.g.,
>> mumps or spooles.
>>
>> Hong
>>
>>>
>>> Here is the log from a successful run.
>>>
>>> [1] 0.249931 Event begin: KSPSetup
>>> [1] 0.24994 Event end: KSPSetup
>>> [1] 0.249948 Event begin: PCSetUp
>>> [1] 0.249963 Event begin: MatGetOrdering
>>> [0] 0.151691 Event begin: KSPSetup
>>> [0] 0.151699 Event end: KSPSetup
>>> [0] PCSetUp(): Setting up new PC
>>> [0] 0.151713 Event begin: PCSetUp
>>> [0] 0.151722 Event begin: MatGetOrdering
>>> [1] PetscCommDuplicate(): returning tag 2147483493
>>> [0] PetscCommDuplicate(): returning tag 2147483493
>>> [1] PetscCommDuplicate(): returning tag 2147483492
>>> [0] PetscCommDuplicate(): returning tag 2147483492
>>> [1] 0.250169 Event end: MatGetOrdering
>>> [1] 0.250182 Event begin: MatLUFactorSym
>>> [0] 0.151885 Event end: MatGetOrdering
>>> [0] 0.151896 Event begin: MatLUFactorSym
>>> [1] PetscCommDuplicate(): returning tag 2147483491
>>> [0] PetscCommDuplicate(): returning tag 2147483491
>>> [1] MatConvert_AIJ_SuperLU_DIST(): Using SuperLU_DIST for SeqAIJ
>> LU factorization and solves.
>>> [0] MatConvert_AIJ_SuperLU_DIST(): Using SuperLU_DIST for SeqAIJ
>> LU factorization and solves.
>>> [1] PetscCommDuplicate(): Using internal PETSc communicator
>> 1083412768 143685848
>>> [1] PetscCommDuplicate(): returning tag 2147483627
>>> [1] PetscCommDuplicate(): Using internal PETSc communicator
>> 1083412768 143685848
>>> [1] PetscCommDuplicate(): returning tag 2147483626
>>> [0] PetscCommDuplicate(): Using internal PETSc communicator
>> 1083412768 144562608
>>> [0] PetscCommDuplicate(): returning tag 2147483627
>>> [0] PetscCommDuplicate(): Using internal PETSc communicator
>> 1083412768 144562608
>>> [0] PetscCommDuplicate(): returning tag 2147483626
>>> [1] 0.252072 Event end: MatLUFactorSym
>>> [1] 0.252098 Event begin: MatLUFactorNum
>>> [0] 0.153821 Event end: MatLUFactorSym
>>> [0] 0.153835 Event begin: MatLUFactorNum
>>> Nonzeros in L 171877
>>> Nonzeros in U 171877
>>> nonzeros in L+U-I 342195
>>> nonzeros in LSUB 52137
>>> Mat conversion(PETSc->SuperLU_DIST) time (max/min/avg):
>>> 0.000296116 / 0.00028491 / 0.000290513
>>> [1] 0.623389 Event end: MatLUFactorNum
>>> [1] 0.623425 Event end: PCSetUp
>>> EQUIL time 0.00
>>> ROWPERM time 0.01
>>> COLPERM time 0.01
>>> SYMBFACT time 0.01
>>> DISTRIBUTE time 0.02
>>> FACTOR time 0.29
>>> Factor flops 2.652935e+08 Mflops 924.82
>>> SOLVE time 0.00
>>> [0] 0.525149 Event end: MatLUFactorNum
>>> [0] 0.525167 Event end: PCSetUp
>>>
>>> Here is the log from a crash.
>>>
>>> [0] 1.41665 Event begin: KSPSetup
>>> [0] 1.41666 Event end: KSPSetup
>>> [0] PCSetUp(): Setting up new PC
>>> [0] 1.41668 Event begin: PCSetUp
>>> [0] 1.41669 Event begin: MatGetOrdering
>>> [1] 1.51505 Event begin: KSPSetup
>>> [1] 1.51506 Event end: KSPSetup
>>> [1] 1.51507 Event begin: PCSetUp
>>> [1] 1.51507 Event begin: MatGetOrdering
>>> [0] PetscCommDuplicate(): returning tag 2147483226
>>> [1] PetscCommDuplicate(): returning tag 2147483226
>>> [0] PetscCommDuplicate(): returning tag 2147483225
>>> [1] PetscCommDuplicate(): returning tag 2147483225
>>> [0] 1.4169 Event end: MatGetOrdering
>>> [0] 1.41692 Event begin: MatLUFactorSym
>>> [1] 1.51526 Event end: MatGetOrdering
>>> [1] 1.51526 Event begin: MatLUFactorSym
>>> [0] PetscCommDuplicate(): returning tag 2147483224
>>> [1] PetscCommDuplicate(): returning tag 2147483224
>>> [0] MatConvert_AIJ_SuperLU_DIST(): Using SuperLU_DIST for SeqAIJ
>> LU factorization and solves.
>>> [1] MatConvert_AIJ_SuperLU_DIST(): Using SuperLU_DIST for SeqAIJ
>> LU factorization and solves.
>>> [0] PetscCommDuplicate(): Using internal PETSc communicator
>> 1083412768 144562608
>>> [0] PetscCommDuplicate(): returning tag 2147483627
>>> [0] PetscCommDuplicate(): Using internal PETSc communicator
>> 1083412768 144562608
>>> [0] PetscCommDuplicate(): returning tag 2147483626
>>> [1] PetscCommDuplicate(): Using internal PETSc communicator
>> 1083412768 143685848
>>> [1] PetscCommDuplicate(): returning tag 2147483627
>>> [1] PetscCommDuplicate(): Using internal PETSc communicator
>> 1083412768 143685848
>>> [1] PetscCommDuplicate(): returning tag 2147483626
>>> [0] 1.4187 Event end: MatLUFactorSym
>>> [0] 1.41872 Event begin: MatLUFactorNum
>>> [1] 1.51706 Event end: MatLUFactorSym
>>> [1] 1.51707 Event begin: MatLUFactorNum
>>> MPI_Alltoallv: invalid datatype argument: Invalid argument (rank 0,
> comm 5)
>>> Rank (0, MPI_COMM_WORLD): Call stack within LAM:
>>> Rank (0, MPI_COMM_WORLD): - MPI_Alltoallv()
>>> Rank (0, MPI_COMM_WORLD): - main()
>>>
>>> If I trace the place of the error in the debugger, the trace gives
>>>
>>> KSPSolve() -> KSPSetUp() -> PCSetUp() -> PCSetUp_LU() ->
>> MatLuFactorNumeric() -> MatLuFactorNumeric_SuperLU_Dist()
>>> -> pzgssvx() -> pzCompRow_loc_to_CompCol_global() -> MPI_Alltoallv
>>>
>>>
>>> I could not make much head-way by looking at the errors. Could you
>> give me some tips on what might be causing this error, and what to look
> for ?
>>>
>>>
>>> Thanks
>>>
>>> Rgds,
>>> Amit
>>>
>>> owner-petsc-users at mcs.anl.gov wrote on 03/03/2008 10:55:21 AM:
>>>
>>>>
>>>> Amit,
>>>>
>>>>> Thanks for your response. A couple of follow-up questions -
>>>>>
>>>>> I go over the steps
>>>>>
>>>>> MatCreate()
>>>>> MatSetType(MPIAIJ)
>>>>> MatMPIAIJSetPreallocation()
>>>>> MatSetType(MATSUPERLU_DIST).
>>>>>
>>>>> Even though I want to do the direct LU solves repeatedly (with the
> same
>>>>> matrix), I don't want the program to do LU factorization repeatedly.
> I
>>> hope
>>>>> I can get that functionality by using the "SAME_PRECONDITIONER"
> option
>>>>> (along with -ksppreonly), while defining the KSP. When the program
> does
>>> the
>>>>> factorization, does it do it in-place or does it allocate new storage
> ?
>>>>
>>>> allocate new storage unless you call
>>>> MatLUFactor() explicitly.
>>>>
>>>> It is done by
>>>> KSPSetUp()
>>>> which calls PCSetUp() -> PCSetUp_LU()
>>>> -> MatLUFactorSymbolic() & MatLUFactorNumeric(
>>>>
>>>> You can do following for reusing factored L and U:
>>>> KSPCreate()
>>>> KSPSetOperators()
>>>> KSPSetFromOptions()
>>>> KSPSetUp()
>>>> while ( num_rhs-- ) {
>>>> KSPSolve(ksp,b,x);
>>>> }
>>>>
>>>> See src/ksp/ksp/examples/tutorials/ex10.c
>>>>
>>>>> After doing the LU factorization, are the other operations such as
>>>>> MatMult() preserved ?
>>>>
>>>> Yes.
>>>>
>>>> Hong
>>>>
>>>>>
>>>>> Thanks
>>>>>
>>>>> Rgds,
>>>>> Amit
>>>>>
>>>>>
>>>
>>>
>
>
More information about the petsc-users
mailing list