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            u, x, xdd, b;      /* exact solution, approx solutions, RHS */
  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      dd_gasm = PETSC_FALSE, dd_asm = PETSC_FALSE, flg = PETSC_FALSE, printx;
  char           dd_type[11];
  PetscScalar    a=1.0+PETSC_i;
  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,PETSC_ERR_SUP,"This example requires complex scalars.");
#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);
  ierr = MatSetUp(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 = VecDuplicate(b,&xdd);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 PC(G)ASM.
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  
  /* 
     Use the -dd_ prefix to set the options for the DD solver: -dd_ksp_type, etc. 
     Use prefix sub to set the subdomain solver options:
      -dd_sub_pc_type <pc_type> -dd_sub_ksp_type <ksp_type> -dd_sub_ksp_rtol 1.e-4 
  */
  ierr = KSPCreate(PETSC_COMM_WORLD,&kspdd);CHKERRQ(ierr);
  ierr = PetscObjectSetOptionsPrefix((PetscObject)kspdd, "dd_"); CHKERRQ(ierr);
  ierr = KSPSetOperators(kspdd,A,A,SAME_PRECONDITIONER);CHKERRQ(ierr);
  ierr = KSPGetPC(kspdd,&pcdd);CHKERRQ(ierr);
  ierr = PetscObjectSetOptionsPrefix((PetscObject)pcdd, "dd_"); CHKERRQ(ierr);
  dd_type[0] = 0;
  ierr = PetscOptionsGetString(PETSC_NULL, "-dd_type", dd_type, 10, &flg); CHKERRQ(ierr);
  ierr = PetscStrcmp(dd_type, "gasm", &dd_gasm); CHKERRQ(ierr);
  if(dd_gasm) {
    ierr = PCSetType(pcdd,PCGASM);CHKERRQ(ierr);
  }
  else {
    ierr = PetscStrcmp(dd_type, "asm", &dd_asm); CHKERRQ(ierr);
    if(dd_asm) {
      ierr = PCSetType(pcdd,PCASM);CHKERRQ(ierr);
    }
    else {
      SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Unsupported -dd_type: %s", dd_type);
    }
  }
 
  /*
     Set the two subdomains for PC(G)ASM
  */

  /* 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 PC(G)ASM with the subdomains */
  if(dd_gasm) {
    ierr= PCGASMSetSubdomains(pcdd,nsub,is_local, is);CHKERRQ(ierr);
  }
  else {
    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 the DD solver options; use the -dd_ prefix: -dd_ksp_type, etc.  
   */
  ierr= KSPSetFromOptions(kspdd);CHKERRQ(ierr);
     


  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
                Create the benchmark linear solver.
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

  /* 
     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);


  /* 
     Use command-line options to set the benchmark KSP and PC options.
  */
#if 0
  /*
     Set type of KSP. 
  */  
  ierr= KSPSetType(ksp,KSPPREONLY);CHKERRQ(ierr);
  {
    PC pc;
    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);
    }
  }
#endif
  
  /* 
    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 with two different solvers
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
                      Check solution x
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

  /*
      Print the first 3 entries of x; this demonstrates extraction of the
      real and imaginary components of the complex vector, x.
  */
  printx  = PETSC_FALSE;
  ierr = PetscOptionsGetBool(PETSC_NULL,"-print_x",&printx,PETSC_NULL);CHKERRQ(ierr);
  if (printx) {
    PetscInt nentries = 3, xlsize;
    ierr = PetscOptionsGetInt(PETSC_NULL, "-print_x_entries", &nentries,PETSC_NULL); CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"The first %D entries of x on each processor are:\n", nentries);CHKERRQ(ierr);
    ierr = VecGetLocalSize(x,&xlsize); CHKERRQ(ierr);
    nentries = PetscMin(xlsize, nentries);
    ierr = VecGetArray(x,&xa);CHKERRQ(ierr);
    for (i=0; i<nentries; i++){
      ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%D:%D]: x[%D] = %G + %G i\n",rank, size, i,PetscRealPart(xa[i]),PetscImaginaryPart(xa[i]));CHKERRQ(ierr);
  }
    ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD); 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(ksp,&its);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error in x after %D iterations: %G\n",its, norm);CHKERRQ(ierr);


  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
                      Check solution xdd
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

  /*
      Print the first 3 entries of xdd; this demonstrates extraction of the
      real and imaginary components of the complex vector, x.
  */
  printx  = PETSC_FALSE;
  ierr = PetscOptionsGetBool(PETSC_NULL,"-dd_print_x",&printx,PETSC_NULL);CHKERRQ(ierr);
  if (printx) {
    PetscInt nentries = 3, xlsize;
    ierr = PetscOptionsGetInt(PETSC_NULL, "-dd_print_x_entries", &nentries,PETSC_NULL); CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"The first %D entries of xdd on each processor are:\n", nentries);CHKERRQ(ierr);
    ierr = VecGetLocalSize(xdd,&xlsize); CHKERRQ(ierr);
    nentries = PetscMin(xlsize, nentries);
    ierr = VecGetArray(xdd,&xa);CHKERRQ(ierr);
    for (i=0; i<nentries; i++){
      ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%D:%D]: xdd[%D] = %G + %G i\n",rank, size, i,PetscRealPart(xa[i]),PetscImaginaryPart(xa[i]));CHKERRQ(ierr);
    }
    ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD); CHKERRQ(ierr);
    ierr = VecRestoreArray(xdd,&xa);CHKERRQ(ierr);
  }
  /* Check the error. */
  ierr = VecAXPY(xdd,none,u);CHKERRQ(ierr);
  ierr = VecNorm(xdd,NORM_2,&norm);CHKERRQ(ierr);
  ierr = KSPGetIterationNumber(kspdd,&its);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error in xdd after %D iterations: %G\n",its, norm);CHKERRQ(ierr);


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