[petsc-users] Solving complex linear sparse matrix in parallel + external library

[Manuel Jung]

Hi,

I can comment on 2):
Petsc does not directly allow to set the minimum iterations, but you can
define your own Convergence Routine (See the manual).
You can define the routine similar to this (Fortran) code:

SUBROUTINE MyConvergenceTest(ksp, it, rnorm, reason, ctx,
ierr)
  use
PetscInterfaces
  use
Cool_petsc
  use
Cool_data
  implicit
none
  type(tKSP) ::
ksp
  PetscInt
it
  PetscFortranAddr
ctx
  PetscErrorCode
ierr
  PetscReal
rnorm
  KSPConvergedReason :: reason

  call
KSPConvergedDefault(ksp,it,rnorm,reason,ctx,ierr)
  IF(it.lt.cool_minits)
THEN
    ierr =
0
    reason =
0
  END
IF
END SUBROUTINE
MyConvergenceTest

You might need to translate that to C in your case and don't forget tell to
tell petsc about your custom Convergence test.
Cheers, Manuel

Am Fr., 16. Nov. 2018 um 10:23 Uhr schrieb Sal Am via petsc-users <
petsc-users at mcs.anl.gov>:

> Hi,
>
> I have a few questions:
>
> 1. The following issue/misunderstanding:
> My code reads in two files one PETSc vector and one PETSc matrix (b and A
> from Ax=b, size ~65000x65000).
> and then calls KSP solver to solve it by running the following in the
> terminal:
>
> mpiexec - n 2 ./SolveSys -ksp_type preonly -pc_type lu
> -pc_factor_mat_solver mumps
>
> Now mumps is supposed to work in parallel and complex, but the code is not
> solved in parallel it seems. It just prints the result twice. Adding
> -log_view gives me
>
> "./SolveSys on a linux-opt named F8434 with 1 processor..." printed twice.
>
> 2. Using iterative solvers, I am having difficulty getting convergence. I
> found that there is a way to set the maximum number of iterations, but is
> there a minimum I can increase?
>
> 3. The residual is not computed when using direct external solvers. What
> is proper PETSc way of doing this?
> -----------------------------------------------------The
> code---------------------------------------
> #include <petscksp.h>
> #include <iostream>
> int main(int argc,char **args)
> {
>   Vec            x,b;      /* approx solution, RHS */
>   Mat            A;            /* linear system matrix */
>   KSP            ksp;         /* linear solver context */
>   PetscReal      norm;         /* norm of solution error */
>   PC             pc;
>   PetscMPIInt    rank, size;
>   PetscViewer    viewer;
>   PetscInt       its, i;
>   PetscErrorCode ierr;
>   PetscScalar    *xa;
>   PetscBool      flg = PETSC_FALSE;
>
>   ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
>   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
>   MPI_Comm_size(PETSC_COMM_WORLD,&size);
>
>   #if !defined(PETSC_USE_COMPLEX)
>     SETERRQ(PETSC_COMM_WORLD,1,"This example requires complex numbers");
>   #endif
>   /*
>    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
>          Compute the matrix and right-hand-side vector that define
>          the linear system, Ax = b.
>      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
>   */
>   ierr = PetscPrintf(PETSC_COMM_WORLD,"reading vector in binary from
> Vector_b.dat ...\n");CHKERRQ(ierr);
>   ierr =
> PetscViewerBinaryOpen(PETSC_COMM_WORLD,"../../python/petscpy/Vector_b.dat",FILE_MODE_READ,&viewer);CHKERRQ(ierr);
>   ierr = VecCreate(PETSC_COMM_WORLD, &b);CHKERRQ(ierr);
>   ierr = VecLoad(b,viewer); CHKERRQ(ierr);
>
>   ierr = PetscPrintf(PETSC_COMM_WORLD,"reading matrix in binary from
> Matrix_A.dat ...\n");CHKERRQ(ierr);
>   ierr =
> PetscViewerBinaryOpen(PETSC_COMM_WORLD,"../../python/petscpy/Matrix_A.dat",FILE_MODE_READ,&viewer);CHKERRQ(ierr);
>   ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
>   ierr = MatLoad(A,viewer);CHKERRQ(ierr);
>   ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
>
>   ierr = VecDuplicate(b,&x);CHKERRQ(ierr);
>
>   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
>                 Create the linear solver and set various options
>      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
>   */
>   PetscPrintf(PETSC_COMM_WORLD, "Creating KSP\n");
>   ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
>
>
>   PetscPrintf(PETSC_COMM_WORLD, "KSP Operators\n");
>   ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
>
>   ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr);
>   ierr = PCSetType(pc, PCLU);CHKERRQ(ierr);
>   ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
>
>   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
>                       Solve the linear system
>      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
>   ierr = KSPSetUp(ksp);CHKERRQ(ierr);
>   ierr = KSPSetUpOnBlocks(ksp);CHKERRQ(ierr);
>   ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
>   PetscPrintf(PETSC_COMM_WORLD, "Solved");
>
>   /*
>      Free work space.  All PETSc objects should be destroyed when they
>      are no longer needed.
>   */
>   ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
>   ierr = VecDestroy(&x);CHKERRQ(ierr);
>   ierr = VecDestroy(&b);CHKERRQ(ierr);
>   ierr = MatDestroy(&A);CHKERRQ(ierr);
>   ierr = PetscFinalize();
>   return ierr;
> }
>
>
> Kind regards,
> Sal
>


-- 
Manuel Jung, Dr. rer. nat.
Hamburger Sternwarte
Universität Hamburg
Gojenbergsweg 112
21029 Hamburg
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20181116/7fc37221/attachment-0001.html>