Direct LU solver
Amit.Itagi at seagate.com
Amit.Itagi at seagate.com
Tue Mar 4 11:36:36 CST 2008
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