/*$Id: ex4.c,v 1.45 2001/08/07 03:03:43 balay Exp $*/

static char help[] = "Demonstrates the use of fast Richardson for SOR. And tests\n\
the MatRelax() routines.\n\n";

#include "petscpc.h"

#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **args)
{
  Mat       mat;
  Vec       b,u;
  PC        pc;
  int       ierr,n = 5,i,col[3];
  PetscScalar    value[3],zero = 0.0;

  PetscInitialize(&argc,&args,(char *)0,help);

  /* Create vectors */
  ierr = VecCreateSeq(PETSC_COMM_SELF,n,&b);CHKERRQ(ierr);
  ierr = VecCreateSeq(PETSC_COMM_SELF,n,&u);CHKERRQ(ierr);

  /* Create and assemble matrix */
  ierr = MatCreateSeqDense(PETSC_COMM_SELF,n,n,PETSC_NULL,&mat);CHKERRQ(ierr);
  value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
  for (i=1; i<n-1; i++) {
    col[0] = i-1; col[1] = i; col[2] = i+1;
    ierr = MatSetValues(mat,1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr);
  }
  i = n - 1; col[0] = n - 2; col[1] = n - 1;
  ierr = MatSetValues(mat,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr);
  i = 0; col[0] = 0; col[1] = 1; value[0] = 2.0; value[1] = -1.0;
  ierr = MatSetValues(mat,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr);
  ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

  /* Create PC context and set up data structures */
  ierr = PCCreate(PETSC_COMM_WORLD,&pc);CHKERRQ(ierr);
  ierr = PCSetType(pc,PCSOR);CHKERRQ(ierr);
  ierr = PCSetFromOptions(pc);CHKERRQ(ierr);
  ierr = PCSetOperators(pc,mat,mat,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
  ierr = PCSetVector(pc,u);CHKERRQ(ierr);
  ierr = PCSetUp(pc);CHKERRQ(ierr);

  value[0] = 1.0;
  for (i=0; i<n; i++) {
    ierr = VecSet(&zero,u);CHKERRQ(ierr);
    ierr = VecSetValues(u,1,&i,value,INSERT_VALUES);CHKERRQ(ierr);
    ierr = PCApply(pc,u,b);CHKERRQ(ierr);
    ierr = VecView(b,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
  }

  /* Free data structures */
  ierr = MatDestroy(mat);CHKERRQ(ierr);
  ierr = PCDestroy(pc);CHKERRQ(ierr);
  ierr = VecDestroy(u);CHKERRQ(ierr);
  ierr = VecDestroy(b);CHKERRQ(ierr); 
  ierr = PetscFinalize();CHKERRQ(ierr);
  return 0;
}