/* Program usage:  mpiexec ex1 [-help] [all PETSc options] */

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

/*T
   Concepts: KSP^solving a system of linear equations
   Processors: 1
T*/

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

  Note:  The corresponding parallel example is ex23.c
*/
#include <petscksp.h>

#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **args)
{
  Vec            x, b, u;      /* approx solution, RHS, exact solution */
  Mat            A;            /* linear system matrix */
  KSP            ksp;         /* linear solver context */
  PC             pc;           /* preconditioner context */
  PetscReal      norm;         /* norm of solution error */
  PetscErrorCode ierr;
  PetscInt       i,n = 18,col[18],its;
  PetscMPIInt    size;
  PetscScalar    neg_one = -1.0,one = 1.0;
  PetscBool      nonzeroguess = PETSC_FALSE;
  PetscScalar   v0[18],v1[18],v2[18],v3[18],v4[18],v5[18],v6[18],v7[18],R[18];
  PetscScalar v8[18],v9[18],v10[18],v11[18],v12[18],v13[18],v14[18],v15[18],v16[18],v17[18];
 
  //v1 = {1.094137824042913e+04,1.611011381962054e+02,0,0,

  v0={ 10941.378240, 161.101138, 0.000000, 0.000000, 3412.821923, 189518489.348116, -10689.706404, -159.588750, 0.000000, 0.000000, -3382.980017, -202964269.020172, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, };
  v1={ -1235365.658038, 4523848.320728, 0.000000, 0.000000, -629738.943484, -43108874616.217880, 1228.903145, -19546.862832, 0.000000, 0.000000, -28079.890056, -4736882654.909527, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, };
  v2={ -5384820.139963, 169197.651949, 5487897.953266, 0.000000, -447641.135471, 170445063752.133514, 5492.740983, -2738.325628, 0.000000, 0.000000, -41924.276374, -7294123900.730478, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, };
  v3={ 0.000000, 169760.674947, 0.000000, 4113844.014604, -460777.543918, -43806085092.764305, 0.000000, -2132.110663, 0.000000, 0.000000, -28787.867928, -4868027472.556285, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, };
  v4={ -199828752.126023, 172269.840696, 0.000000, 0.000000, 4985880.267183, 7425409382139.539062, 27396507.746542, 5512.078058, 0.000000, 0.000000, -53741.687760, -914450738791.886230, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, };
  v5={ -9499.125441, -127.976567, 0.000000, 0.000000, -2713.813160, 28868833452.505474, 5.656186, 0.527407, 0.000000, 0.000000, 11.371353, -74719210.450024, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, };
  v6={ -10689.696942, -159.588060, 0.000000, 0.000000, -3382.965917, -202962262.885559, 11007.672933, 161.196949, 0.000000, 0.000000, 3413.834628, 185274309.165660, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, };
  v7={ 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, -1130479.394956, 4774054.077672, 0.000000, 0.000000, -452332.493610, 55532905777.686256, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, };
  v8={ 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, -7950763.247911, 177728.522230, 7464458.469369, 0.000000, -339884.891533, 296793011460.141052, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, };
  v9={ 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 179311.149314, 0.000000, 4263268.045563, -339884.891533, 1782066591.634391, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, };
  v10={ 286011.958324, 83.493934, 0.000000, 0.000000, 0.000000, -9454349903.383270, -179266693.994530, 187686.144566, 0.000000, 0.000000, 5413376.728119, 8070741526111.384766, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, };
  v11={ 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, -53200000.000000, -11884.635249, -127.639110, 0.000000, 0.000000, -2672.466150, 26358195654.878288, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, };
  v12={ 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, };
  v13={ 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1.000000, 0.000000, 0.000000, 0.000000, 0.000000, };
  v14={ 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1.000000, 0.000000, 0.000000, 0.000000, };
  v15={ 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1.000000, 0.000000, 0.000000, };
  v16={ 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1.000000, 0.000000, };
  v17={ 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1.000000, };

 
  R = {-2.0830395435112612e+04, -7.7131523863715199e+02,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       -1.6545026895435876e+04,
       2.8322037514216010e+11,
       2.0830488048237599e+04,
       7.7128820369681557e+02,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       1.6545026895435876e+04,
       8.3672094169114275e+09,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       0.0000000000000000e+00,
       0.0000000000000000e+00};
  PetscInitialize(&argc,&args,(char *)0,help);
  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
  if (size != 1) SETERRQ(PETSC_COMM_WORLD,1,"This is a uniprocessor example only!");
  ierr = PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);CHKERRQ(ierr);
  ierr = PetscOptionsGetBool(PETSC_NULL,"-nonzero_guess",&nonzeroguess,PETSC_NULL);CHKERRQ(ierr);


  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
         Compute the matrix and right-hand-side vector that define
         the linear system, Ax = b.
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

  /* 
     Create vectors.  Note that we form 1 vector from scratch and
     then duplicate as needed.
  */
  ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject) x, "Solution");CHKERRQ(ierr);
  ierr = VecSetSizes(x,PETSC_DECIDE,n);CHKERRQ(ierr);
  ierr = VecSetFromOptions(x);CHKERRQ(ierr);
  ierr = VecDuplicate(x,&b);CHKERRQ(ierr);
  ierr = VecDuplicate(x,&u);CHKERRQ(ierr);

  /* 
     Create matrix.  When using MatCreate(), the matrix format can
     be specified at runtime.

     Performance tuning note:  For problems of substantial size,
     preallocation of matrix memory is crucial for attaining good 
     performance. See the matrix chapter of the users manual for details.
  */
  ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
  ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr);
  //ierr = MatSetOption(A,MAT_ROW_ORIENTED,PETSC_FALSE);CHKERRQ(ierr);
  ierr = MatSetFromOptions(A);CHKERRQ(ierr);

  /* 
     Assemble matrix
  */
  /* 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(A,1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr); */
  /* } */
  /* i = n - 1; col[0] = n - 2; col[1] = n - 1; */
  /* ierr = MatSetValues(A,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(A,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr); */
  for (i=0;i<n;i++) col[i] = i;
  i=0;
  ierr = MatSetValues(A,18,col,1,&i,v0,INSERT_VALUES);CHKERRQ(ierr);
  i++;
  ierr = MatSetValues(A,18,col,1,&i,v1,INSERT_VALUES);CHKERRQ(ierr);
  i++;
  ierr = MatSetValues(A,18,col,1,&i,v2,INSERT_VALUES);CHKERRQ(ierr);
  i++;
  ierr = MatSetValues(A,18,col,1,&i,v3,INSERT_VALUES);CHKERRQ(ierr);i++;
  ierr = MatSetValues(A,18,col,1,&i,v4,INSERT_VALUES);CHKERRQ(ierr);i++;
  ierr = MatSetValues(A,18,col,1,&i,v5,INSERT_VALUES);CHKERRQ(ierr);i++;
  ierr = MatSetValues(A,18,col,1,&i,v6,INSERT_VALUES);CHKERRQ(ierr);i++;
  ierr = MatSetValues(A,18,col,1,&i,v7,INSERT_VALUES);CHKERRQ(ierr);i++;
  ierr = MatSetValues(A,18,col,1,&i,v8,INSERT_VALUES);CHKERRQ(ierr);i++;
  ierr = MatSetValues(A,18,col,1,&i,v9,INSERT_VALUES);CHKERRQ(ierr);i++;
  ierr = MatSetValues(A,18,col,1,&i,v10,INSERT_VALUES);CHKERRQ(ierr);i++;
  ierr = MatSetValues(A,18,col,1,&i,v11,INSERT_VALUES);CHKERRQ(ierr);i++;
  ierr = MatSetValues(A,18,col,1,&i,v12,INSERT_VALUES);CHKERRQ(ierr);i++;
  ierr = MatSetValues(A,18,col,1,&i,v13,INSERT_VALUES);CHKERRQ(ierr);i++;
  ierr = MatSetValues(A,18,col,1,&i,v14,INSERT_VALUES);CHKERRQ(ierr);i++;
  ierr = MatSetValues(A,18,col,1,&i,v15,INSERT_VALUES);CHKERRQ(ierr);i++;
  ierr = MatSetValues(A,18,col,1,&i,v16,INSERT_VALUES);CHKERRQ(ierr);i++;
  ierr = MatSetValues(A,18,col,1,&i,v17,INSERT_VALUES);CHKERRQ(ierr);i++;

  ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

  /* 
     Set exact solution; then compute right-hand-side vector.
  */
  ierr = VecSet(u,one);CHKERRQ(ierr);
  ierr = VecSetValues(b,18,col,R,INSERT_VALUES);CHKERRQ(ierr);
  ierr = VecAssemblyBegin(b);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(b);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 linear solver defaults for this problem (optional).
     - By extracting the KSP and PC contexts from the KSP context,
       we can then directly call any KSP and PC routines to set
       various options.
     - The following four statements are optional; all of these
       parameters could alternatively be specified at runtime via
       KSPSetFromOptions();
  */
  ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
  ierr = PCSetType(pc,PCJACOBI);CHKERRQ(ierr);
  ierr = KSPSetTolerances(ksp,1.e-9,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr);

  /* 
    Set runtime options, e.g.,
        -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
    These options will override those specified above as long as
    KSPSetFromOptions() is called _after_ any other customization
    routines.
  */
  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);

  if (nonzeroguess) {
    PetscScalar p = .5;
    ierr = VecSet(x,p);CHKERRQ(ierr);
    ierr = KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);CHKERRQ(ierr);
  }
 
  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
                      Solve the linear system
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  /* 
     Solve linear system
  */
  ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr); 

  /* 
     View solver info; we could instead use the option -ksp_view to
     print this info to the screen at the conclusion of KSPSolve().
  */
  ierr = KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
                      Check solution and clean up
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  /* 
     Check the error
  */
  ierr = VecAXPY(x,neg_one,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 %A, Iterations %D\n",
                     norm,its);CHKERRQ(ierr);
  MatView(A,0);
  VecView(b,0);
  VecView(x,0);
  /* 
     Free work space.  All PETSc objects should be destroyed when they
     are no longer needed.
  */
  ierr = VecDestroy(&x);CHKERRQ(ierr); ierr = VecDestroy(&u);CHKERRQ(ierr);
  ierr = VecDestroy(&b);CHKERRQ(ierr); ierr = MatDestroy(&A);CHKERRQ(ierr);
  ierr = KSPDestroy(&ksp);CHKERRQ(ierr);

  /*
     Always call PetscFinalize() before exiting a program.  This routine
       - finalizes the PETSc libraries as well as MPI
       - provides summary and diagnostic information if certain runtime
         options are chosen (e.g., -log_summary).
  */
  ierr = PetscFinalize();
  return 0;
}