[petsc-users] use Superlu as ilu preconditioner

[Ping Rong]

Dear Hong,

I have found the problem. I don't know whether this is a bug.

	ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat(),
         		               this->same_preconditioner ? SAME_PRECONDITIONER : DIFFERENT_NONZERO_PATTERN);


I used my system matrix as the precondition matrix, where matrix->mat() 
and precond->mat() are in deed the same matrix. I have left the 
KSPSetOperators(ksp) out before the KSPSolve(..) function, so 
KSPSetOperators() is called only once while the petsc_linear_solver 
object is initialized, which eliminated the petsc error ("object is in 
wrong state").
Also when I used petsc default ilu preconditioner, everything is fine, 
also the solutions are correct. However, as long as I apply the superlu 
as my solver package, it overwrites the system matrix. I compared the 
system matrix before and after solve(). With petsc ilu, it remains 
identical, but with superlu ilu, the system matrix is changed to the 
preconditioner matrix itself. Now I simply duplicate the system matrix, 
and use this copy as precondition matrix. Everything works perfectly.

The code I used is pretty standard, just as in ex2. It took me quite a 
while to realize the problem, because when the system matrix is replaced 
by the preconditioner matrix, you still get a "correct" solution, which 
is however incorrectly scaled.

Ping


Am 28.07.2011 18:21, schrieb Hong Zhang:
> Ping :
> I tested similar calling procedure using
> petsc-dev/src/ksp/ksp/examples/tutorials/ex5.c
> successfully.
>
> Try following:
> 1) switch to petsc-dev. See
> http://www.mcs.anl.gov/petsc/petsc-as/developers/index.html on how to
> get it.
> 2) send us a short and stand-alone code that reproduce the error. We
> can investigate it.
>
> Hong
>
>> Dear Hong,
>>
>> thank you very much for your time. Can you tell me when should
>>   KSPSetFromOptions(ksp) be called?
>> well, libMesh has a class called "petsc_linear_solver", in its original code
>> KSPSetFromOptions(ksp) is called while the object is created. I couldn't get
>> any output, when I run the program with the option
>>
>>    -pc_type ilu -pc_factor_mat_solver_package superlu
>> -mat_superlu_ilu_droptol 0.1 -ksp_view
>>
>> The actual solve() function of the "petsc_linear_solver" contains the
>> similar code as in the ex2.c
>> ...
>>     ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat(),
>>                                              this->same_preconditioner ?
>> SAME_PRECONDITIONER : DIFFERENT_NONZERO_PATTERN);
>>     CHKERRABORT(libMesh::COMM_WORLD,ierr);
>>
>>   // Set the tolerances for the iterative solver.  Use the user-supplied
>>   // tolerance for the relative residual&  leave the others at default
>> values.
>>   ierr = KSPSetTolerances (_ksp, tol, PETSC_DEFAULT, PETSC_DEFAULT, max_its);
>>   CHKERRABORT(libMesh::COMM_WORLD,ierr);
>>
>>   ierr = KSPSetFromOptions (_ksp);<--------------  I added these two lines,
>>   CHKERRABORT(libMesh::COMM_WORLD,ierr);<--------------  then I got proper
>> output.
>>
>>   // Solve the linear system
>>   ierr = KSPSolve (_ksp, rhs->vec(), solution->vec());
>>   CHKERRABORT(libMesh::COMM_WORLD,ierr);
>> ....
>>
>> The problem is that I have to solve a multiple frequency system, and for
>> each frequency I have several RHS. The same ksp context is used to solve the
>> systems repeatedly. I always call KSPSetOperators() with the option
>> DIFFERENT_NONZERO_PATTERN for the first RHS and SAME_PRECONDITIONER for the
>> rest, since the system matrix remains the same.  For the first frequency,
>> everything is fine, both DIFFERENT_NONZERO_PATTERN and SAME_PRECONDITIONER
>> options. I got following output from the -ksp_view.
>> ...
>> PC Object:
>>   type: ilu
>>     ILU: out-of-place factorization
>>     0 levels of fill
>>     tolerance for zero pivot 1e-12
>>     using diagonal shift to prevent zero pivot
>>     matrix ordering: natural
>>     factor fill ratio given 0, needed 0
>>       Factored matrix follows:
>>         Matrix Object:
>>           type=seqaij, rows=2883, cols=2883
>>           package used to perform factorization: superlu<----------------
>> superlu is properly set
>> ....
>>
>> but when the second frequency comes up, KSPSetOperators() is called again
>> with DIFFERENT_NONZERO_PATTERN option for the first RHS. the command line
>> option doesn't seem to work anymore.
>> ...
>> PC Object:
>>   type: ilu
>>     ILU: out-of-place factorization
>>     0 levels of fill
>>     tolerance for zero pivot 1e-12
>>     using diagonal shift to prevent zero pivot
>>     matrix ordering: natural
>>     factor fill ratio given 1, needed 1
>>       Factored matrix follows:
>>         Matrix Object:
>>           type=seqaij, rows=2883, cols=2883
>>           package used to perform factorization: petsc<------------ not
>> superlu anymore
>> ...
>>
>> and then for the second RHS, when it calls KSPSetOperators(ksp) in the
>> solve() routine, petsc throws the following error message.
>>
>> [0]PETSC ERROR: --------------------- Error Message
>> ------------------------------------
>> [0]PETSC ERROR: Object is in wrong state!
>> [0]PETSC ERROR: Cannot change solver matrix package after PC has been setup
>> or used!
>> [0]PETSC ERROR:
>> ------------------------------------------------------------------------
>> [0]PETSC ERROR: Petsc Release Version 3.1.0, Patch 8, Thu Mar 17 13:37:48
>> CDT 2011
>> [0]PETSC ERROR: See docs/changes/index.html for recent updates.
>> [0]PETSC ERROR: See docs/faq.html for hints about trouble shooting.
>> [0]PETSC ERROR: See docs/index.html for manual pages.
>> [0]PETSC ERROR:
>> ------------------------------------------------------------------------
>> [0]PETSC ERROR: ./PLTMainTest.opt on a linux-gnu named abe-vm by mubpr Wed
>> Jul 27 19:41:22 2011
>> [0]PETSC ERROR: Libraries linked from
>> /home/mubpr/libs/petsc/linux-gnu-acml-shared-opt/lib
>> [0]PETSC ERROR: Configure run at Tue Jul 26 22:09:39 2011
>> [0]PETSC ERROR: Configure options
>> --with-blas-lapack-lib=/home/mubpr/libs/acml/gfortran64_mp/lib/libacml_mp.so
>> --with-blas-lapack-include=/home/mubpr/libs/acml/gfortran64_mp/include
>> --with-scalar-type=complex --with-clanguage=c++
>> --with-mpi-dir=/home/mubpr/libs/mpich2 --download-superlu=yes
>> --with-superlu=1 --download-parmetis=yes --with-parmetis=1
>> --download-superlu_dist=yes --with-superlu_dist=1 --download-mumps=yes
>> --with-mumps=1 --download-blacs=yes --with-blacs=1 --download-scalapack=yes
>> --with-scalapack=1 --download-spooles=yes --with-spooles=1
>> --download-umfpack=yes --with-umfpack=1 --with-debugging=no --with-shared=1
>> [0]PETSC ERROR:
>> ------------------------------------------------------------------------
>> [0]PETSC ERROR: PCFactorSetMatSolverPackage_Factor() line 183 in
>> src/ksp/pc/impls/factor/factimpl.c
>> [0]PETSC ERROR: PCFactorSetMatSolverPackage() line 276 in
>> src/ksp/pc/impls/factor/factor.c
>> [0]PETSC ERROR: PCSetFromOptions_Factor() line 275 in
>> src/ksp/pc/impls/factor/factimpl.c
>> [0]PETSC ERROR: PCSetFromOptions_ILU() line 112 in
>> src/ksp/pc/impls/factor/ilu/ilu.c
>> [0]PETSC ERROR: PCSetFromOptions() line 185 in src/ksp/pc/interface/pcset.c
>> [0]PETSC ERROR: KSPSetFromOptions() line 323 in src/ksp/ksp/interface/itcl.c
>> [0]PETSC ERROR: User provided function() line 696 in
>> "unknowndirectory/"src/solvers/petsc_linear_solver.C
>>
>> Any idea what the problem may be? I would be very grateful, if you can give
>> me any hint. thanks alot
>>
>> best regards
>> Ping
>>
>>
>>
>> Am 27.07.2011 04:13, schrieb Hong Zhang:
>>> Did you call KSPSetFromOptions(ksp)?
>>> Run your code with '-options_table' to dump list of options inputted
>>> or '-options_left' to dump list of unused options.
>>>
>>> I tested with petsc-3.1/src/ksp/ksp/examples/tutorials/ex2.c:
>>> ./ex2
>>> ...
>>>              SuperLU run parameters:
>>>              ...
>>>                ILU_DropTol: 0.1
>>>                ILU_FillTol: 0.01
>>>                ILU_FillFactor: 10
>>>                ILU_DropRule: 9
>>>                ILU_Norm: 2
>>>                ILU_MILU: 2
>>>
>>> Hong
>>>
>>> On Tue, Jul 26, 2011 at 3:53 PM, Ping Rong<ping.rong at tuhh.de>    wrote:
>>>> Dear developers,
>>>>
>>>> I have compiled petsc-3.1-p8 for a while. Now I would like to use superlu
>>>> as
>>>> an ilu preconditioner, since it offers the drop tolerance option. I have
>>>> read in a thread
>>>>
>>>> (https://lists.mcs.anl.gov/mailman/htdig/petsc-users/2010-December/007439.html)
>>>> that one can run the code with option
>>>>
>>>> -pc_type ilu -pc_factor_mat_solver_package superlu
>>>> -mat_superlu_ilu_droptol
>>>> <>
>>>>
>>>> to setup the ilu preconditioner. I also use the option "-help |grep
>>>> superlu"
>>>> to check the settings. However, no matter how I change the value of
>>>> -mat_superlu_ilu_droptol, I always got the same result
>>>> ...
>>>>   -mat_superlu_lwork<0>: size of work array in bytes used by factorization
>>>> (None)
>>>>   -mat_superlu_ilu_droptol<0.0001>: ILU_DropTol (None)
>>>>   -mat_superlu_ilu_filltol<0.01>: ILU_FillTol (None)
>>>> ...
>>>> I dont know whether I understand it correctly. But it seems to me the
>>>> value
>>>> of the droptol has never been changed. In that maillist thread, it was
>>>> also
>>>> mentioned that the dev-version is recommended, because of some bugs in
>>>> the
>>>> superlu interface. I have then compiled the current dev-version. but the
>>>> problem is my program is based on another library (libMesh) which uses
>>>> petsc
>>>> as a solver package. Since some of the interfaces have been changed in
>>>> the
>>>> dev-version, I would not be able to compile the libMesh with petsc
>>>> anymore.
>>>> Is there anyway I can correct the superlu interface in the 3.1-p8 version
>>>> directly? Any help will be appreciated!!
>>>>
>>>> Best regards
>>>>
>>>> --
>>>> Ping Rong, M.Sc.
>>>> Hamburg University of Technology
>>>> Institut of modelling and computation
>>>> Denickestraße 17 (Room 3031)
>>>> 21073 Hamburg
>>>>
>>>> Tel.: ++49 - (0)40 42878 2749
>>>> Fax:  ++49 - (0)40 42878 43533
>>>> Email: ping.rong at tuhh.de
>>>>
>>>>
>>
>> --
>> Ping Rong, M.Sc.
>> Hamburg University of Technology
>> Institut of modelling and computation
>> Denickestraße 17 (Room 3031)
>> 21073 Hamburg
>>
>> Tel.: ++49 - (0)40 42878 2749
>> Fax:  ++49 - (0)40 42878 43533
>> Email: ping.rong at tuhh.de
>>
>>


-- 
Ping Rong, M.Sc.
Hamburg University of Technology
Institut of modelling and computation
Denickestraße 17 (Room 3031)
21073 Hamburg

Tel.: ++49 - (0)40 42878 2749
Fax:  ++49 - (0)40 42878 43533
Email: ping.rong at tuhh.de