static char help[] = "Solves a linear system in parallel with KSP.\n\n";

/*T
   Concepts: KSP^solving a Helmholtz equation
   Concepts: complex numbers;
   Concepts: Helmholtz equation
   Processors: n
T*/

/*
   Description: Solves a complex linear system in parallel with KSP.

   The model problem:
      Solve Helmholtz equation on the unit square: (0,1) x (0,1)
          -delta u - sigma1*u + i*sigma2*u = f, 
           where delta = Laplace operator
      Dirichlet b.c.'s on all sides
      Use the 2-D, five-point finite difference stencil.

   Compiling the code:
      This code uses the complex numbers version of PETSc, so configure
      must be run to enable this
*/

/* 
  Include "petscksp.h" so that we can use KSP solvers.  Note that this file
  automatically includes:
     petscsys.h       - base PETSc routines   petscvec.h - vectors
     petscmat.h - matrices
     petscis.h     - index sets            petscksp.h - Krylov subspace methods
     petscviewer.h - viewers               petscpc.h  - preconditioners
*/
#include <petscksp.h>


#undef __FUNCT__
#define __FUNCT__ "ModifySubMatrices"
PetscErrorCode ModifySubMatrices(PC pc,PetscInt nsub,const IS *row,const IS *col,Mat *submat,void *ctx)
{

  return 0;
}

#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **args)
{
  Vec            x,b,u,udd;      /* approx solution, RHS, exact solution */
  Mat            A;            /* linear system matrix */
  KSP            ksp, kspdd;         /* linear solver context */
  PetscReal      norm;         /* norm of solution error */
  PetscInt       dim,i,j,Ii,J,Istart,Iend,n = 6,its, nsub;
  PetscErrorCode ierr;
  PetscScalar    v,none = -1.0,sigma2,pfive = 0.5,*xa;
  PetscReal      h2,sigma1 = 100.0;
  PetscBool      flg = PETSC_FALSE;
  PetscScalar    a=1.0+PETSC_i;
  PC             pc, pcdd;
  IS             is_local[2], is[2];
  PetscMPIInt    rank, size;

  PetscInitialize(&argc,&args,(char *)0,help);
#if !defined(PETSC_USE_COMPLEX)
  SETERRQ(PETSC_COMM_WORLD,1,"This example requires complex numbers");
#endif

  a=1.0+PETSC_i;
  printf("%G+%Gi\n",PetscRealPart(a),PetscImaginaryPart(a));

  ierr = PetscOptionsGetReal(PETSC_NULL,"-sigma1",&sigma1,PETSC_NULL);CHKERRQ(ierr);
  ierr = PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);CHKERRQ(ierr);
  dim = n*n;

  MPI_Comm_size(PETSC_COMM_WORLD,&size);
  if (size>2) {
    PetscPrintf(PETSC_COMM_WORLD,"Please use 1 or 2 processes! Quit\n");
    return 0;
  }

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
         Compute the matrix and right-hand-side vector that define
         the linear system, Ax = b.
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  /* 
     Create parallel matrix, specifying only its global dimensions.
     When using MatCreate(), the matrix format can be specified at
     runtime. Also, the parallel partitioning of the matrix is
     determined by PETSc at runtime.
  */
  ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
  ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,dim,dim);CHKERRQ(ierr);
  ierr = MatSetFromOptions(A);CHKERRQ(ierr);

  /* 
     Currently, all PETSc parallel matrix formats are partitioned by
     contiguous chunks of rows across the processors.  Determine which
     rows of the matrix are locally owned. 
  */
  ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr);

  /* 
     Set matrix elements in parallel.
      - Each processor needs to insert only elements that it owns
        locally (but any non-local elements will be sent to the
        appropriate processor during matrix assembly). 
      - Always specify global rows and columns of matrix entries.
  */

  sigma2 = 10.0*PETSC_i;
  h2 = 1.0/((n+1)*(n+1));
  for (Ii=Istart; Ii<Iend; Ii++) { 
    v = -1.0; i = Ii/n; j = Ii - i*n;  
    if (i>0) {
      J = Ii-n; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);}
    if (i<n-1) {
      J = Ii+n; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);}
    if (j>0) {
      J = Ii-1; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);}
    if (j<n-1) {
      J = Ii+1; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);}
    v = 4.0 - sigma1*h2 + sigma2*h2;
    ierr = MatSetValues(A,1,&Ii,1,&Ii,&v,ADD_VALUES);CHKERRQ(ierr);
  }


  /* 
     Assemble matrix, using the 2-step process:
       MatAssemblyBegin(), MatAssemblyEnd()
     Computations can be done while messages are in transition
     by placing code between these two statements.
  */
  ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

  /* 
     Create parallel vectors.
      - When using VecCreate(), VecSetSizes() and VecSetFromOptions(),
      we specify only the vector's global
        dimension; the parallel partitioning is determined at runtime. 
      - Note: We form 1 vector from scratch and then duplicate as needed.
  */
  ierr = VecCreate(PETSC_COMM_WORLD,&u);CHKERRQ(ierr);
  ierr = VecSetSizes(u,PETSC_DECIDE,dim);CHKERRQ(ierr);
  ierr = VecSetFromOptions(u);CHKERRQ(ierr);
  ierr = VecDuplicate(u,&b);CHKERRQ(ierr); 
  ierr = VecDuplicate(b,&x);CHKERRQ(ierr);
  ierr = ierr = VecDuplicate(b,&udd);CHKERRQ(ierr);
 
  /* 
     Set exact solution; then compute right-hand-side vector.
  */
  
  ierr = VecSet(u,pfive);CHKERRQ(ierr);
  ierr = MatMult(A,u,b);CHKERRQ(ierr);

   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
                     Create the linear solver using PCGASM 
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  
  /*  use prefix sub to set the same ksp options for the subdomain problems -sub_pc_type <pc> -sub_ksp_type <ksp> -sub_ksp_rtol 1.e-4 */
  ierr = KSPCreate(PETSC_COMM_WORLD,&kspdd);CHKERRQ(ierr);
  ierr = KSPSetOperators(kspdd,A,A,SAME_PRECONDITIONER);CHKERRQ(ierr);
  ierr = KSPGetPC(kspdd,&pcdd);CHKERRQ(ierr);
  ierr = PCSetType(pcdd,PCASM);CHKERRQ(ierr);
 
  /*
     Set the two subdomains for PCASM
  */

  /* non-overlapping subdomains */
  MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
  if (size==2) {
    nsub= 1;
    ISCreateStride(PETSC_COMM_SELF,n*n/2,rank*n*n/2,1,&(is_local[0]));
  } 
  else{
    nsub= 2;
    ISCreateStride(PETSC_COMM_SELF,n*n/2,0,1,&(is_local[0]));
    ISCreateStride(PETSC_COMM_SELF,n*n/2,n*n/2,1,&(is_local[1]));
  } 
  /* printf("[%i]: ",rank); */
  /* ISView(is_local[i],PETSC_VIEWER_STDOUT_SELF); */
  
 
  /* overlapping subdomains */
  PetscInt expand= 1;
  if (size==2) {
    if (rank==0) {
      ISCreateStride(PETSC_COMM_SELF,n*(n/2+expand),0,1,&(is[0]));
    }
    else
      ISCreateStride(PETSC_COMM_SELF,n*(n/2+expand),n*(n/2-expand),1,&(is[0]));
  } 
  else{  /* size==1 */
    ISCreateStride(PETSC_COMM_SELF,n*(n/2+expand),0,1,&(is[0]));
    ISCreateStride(PETSC_COMM_SELF,n*(n/2+expand),n*(n/2-expand),1,&(is[1]));
  }
  /* printf("[%i]: ",rank); */
  /* ISView(is[i],PETSC_VIEWER_STDOUT_SELF); */

  /* set PCASM with the subdomains */
  ierr= PCASMSetLocalSubdomains(pcdd,nsub,is,is_local);CHKERRQ(ierr);
  ierr= PCView(pcdd, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);

  /*
     Modify submatrices
  */
  ierr= PCSetModifySubMatrices(pcdd,ModifySubMatrices,PETSC_NULL);CHKERRQ(ierr);

  /* 
     Set options
   */
  ierr= KSPSetFromOptions(kspdd);CHKERRQ(ierr);


  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
                Create the linear solver and set various options
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

  /* 
     Create linear solver context
  */
  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);

  /* 
     Set operators. Here the matrix that defines the linear system
     also serves as the preconditioning matrix.
  */
  ierr = KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);

  /*
     Set type of KSP. 
  */  
  ierr= KSPSetType(ksp,KSPPREONLY);CHKERRQ(ierr);
  ierr= KSPGetPC(ksp,&pc);CHKERRQ(ierr);
  ierr= PCSetType(pc,PCLU);CHKERRQ(ierr);
  if (size==1) {
    ierr= PCFactorSetMatSolverPackage(pc,MATSOLVERSUPERLU);CHKERRQ(ierr);
  }
  else
    ierr= PCFactorSetMatSolverPackage(pc,MATSOLVERSUPERLU_DIST);CHKERRQ(ierr);
  
  /* 
    Set runtime options, e.g.,
        -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
  */
  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
                      Solve the linear system
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

  ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
  ierr = KSPSolve(kspdd,b,udd);CHKERRQ(ierr);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
                      Check solution and clean up
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

  /*
      Print the first 3 entries of x; this demonstrates extraction of the
      real and imaginary components of the complex vector, x.
  */
  flg  = PETSC_FALSE;
  ierr = PetscOptionsGetBool(PETSC_NULL,"-print_x3",&flg,PETSC_NULL);CHKERRQ(ierr);
  if (flg) {
    ierr = VecGetArray(x,&xa);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"The first three entries of x are:\n");CHKERRQ(ierr);
    for (i=0; i<3; i++){
      ierr = PetscPrintf(PETSC_COMM_WORLD,"x[%D] = %G + %G i\n",i,PetscRealPart(xa[i]),PetscImaginaryPart(xa[i]));CHKERRQ(ierr);
  }
    ierr = VecRestoreArray(x,&xa);CHKERRQ(ierr);
  }

  /* 
     Check the error
  */
  ierr = VecAXPY(x,none,u);CHKERRQ(ierr);
  ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr);
  ierr = KSPGetIterationNumber(kspdd,&its);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %A iterations %D\n",norm,its);CHKERRQ(ierr);

  /* 
     Free work space.  All PETSc objects should be destroyed when they
     are no longer needed.
  */
  ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
  ierr = VecDestroy(&u);CHKERRQ(ierr); ierr = VecDestroy(&x);CHKERRQ(ierr);
  ierr = VecDestroy(&b);CHKERRQ(ierr); ierr = MatDestroy(&A);CHKERRQ(ierr);
  ierr = PetscFinalize();
  return 0;
}