#include "petsc.h"
#include <petscksp.h>
#include <Box.H>
using namespace std;

int load( Mat& a_A, Vec& a_b )
{
#ifdef CH_MPI
  MPI_Comm a_communicator = Chombo_MPI::comm;
#else
  MPI_Comm a_communicator = PETSC_COMM_SELF;
#endif


  double b[28];
  double a[28][28];

  for( int i=0; i<28; i++ ) {
    b[i] = 0;
    for( int j=0; j<28; j++ ) a[i][j] = 0;
  }

b[0]=-0.01547629572451115;
a[0][1]=65.58024691358024;
a[0][0]=-35.16049382716049;
a[0][2]=-1.975308641975309;
a[0][6]=-21.72839506172839;
a[0][8]=1.975308641975309;
a[0][7]=-8.691358024691358;
b[1]=-0.002485675737261772;
a[1][0]=18.21673525377229;
a[1][1]=-28.53223593964335;
a[1][2]=38.7599451303155;
a[1][7]=-28.35665294924554;
a[1][6]=0.7462277091906734;
a[1][8]=-0.8340192043895763;
b[2]=0.0001520376652479477;
a[2][1]=22.29728395061728;
a[2][2]=-28.47604938271605;
a[2][3]=34.62320987654321;
a[2][8]=-28.41283950617284;
a[2][7]=0.45827160493827;
a[2][9]=-0.4898765432098747;
b[3]=-0.0009238988693271262;
a[3][2]=24.050390526581;
a[3][3]=-28.46056941295037;
a[3][4]=32.8546233308138;
a[3][9]=-28.42831947593852;
a[3][8]=0.330561854371377;
a[3][10]=-0.3466868228772978;
b[4]=-0.005844825878739357;
a[4][3]=25.02545343697607;
a[4][4]=-28.4541990550221;
a[4][5]=31.87319006249047;
a[4][10]=-28.43468983386678;
a[4][9]=0.258497180307879;
a[4][11]=-0.2682517908855347;
b[6]=-0.03303125314414501;
a[6][7]=74.27160493827159;
a[6][6]=-13.4320987654321;
a[6][8]=-3.950617283950617;
a[6][0]=-21.72839506172839;
a[6][2]=1.975308641975309;
a[6][1]=-8.691358024691358;
a[6][12]=-21.72839506172839;
a[6][14]=1.975308641975309;
a[6][13]=-8.691358024691358;
b[7]=-0.003860184922814369;
a[7][6]=17.47050754458161;
a[7][7]=-0.1755829903978139;
a[7][8]=39.59396433470508;
a[7][1]=-28.35665294924554;
a[7][0]=0.7462277091906734;
a[7][2]=-0.8340192043895763;
a[7][13]=-28.35665294924554;
a[7][12]=0.7462277091906734;
a[7][14]=-0.8340192043895763;
b[8]=0.001571375504136135;
a[8][7]=21.83901234567902;
a[8][8]=-0.06320987654320476;
a[8][9]=35.11308641975308;
a[8][2]=-28.41283950617284;
a[8][1]=0.45827160493827;
a[8][3]=-0.4898765432098747;
a[8][14]=-28.41283950617284;
a[8][13]=0.45827160493827;
a[8][15]=-0.4898765432098747;
b[9]=-0.003092412171619249;
a[9][8]=23.71982867220962;
a[9][9]=-0.03224993701184786;
a[9][10]=33.2013101536911;
a[9][3]=-28.42831947593852;
a[9][2]=0.330561854371377;
a[9][4]=-0.3466868228772978;
a[9][15]=-28.42831947593852;
a[9][14]=0.330561854371377;
a[9][16]=-0.3466868228772978;
b[12]=-0.01821462996304035;
a[12][13]=74.27160493827159;
a[12][12]=-13.4320987654321;
a[12][14]=-3.950617283950617;
a[12][6]=-21.72839506172839;
a[12][8]=1.975308641975309;
a[12][7]=-8.691358024691358;
a[12][17]=-21.72839506172839;
a[12][19]=1.975308641975309;
a[12][18]=-8.691358024691358;
b[13]=0.00294872559607029;
a[13][12]=17.47050754458161;
a[13][13]=-0.1755829903978139;
a[13][14]=39.59396433470508;
a[13][7]=-28.35665294924554;
a[13][6]=0.7462277091906734;
a[13][8]=-0.8340192043895763;
a[13][18]=-28.35665294924554;
a[13][17]=0.7462277091906734;
a[13][19]=-0.8340192043895763;
b[14]=0.004772236570715904;
a[14][13]=21.83901234567902;
a[14][14]=-0.06320987654320476;
a[14][15]=35.11308641975308;
a[14][8]=-28.41283950617284;
a[14][7]=0.45827160493827;
a[14][9]=-0.4898765432098747;
a[14][19]=-28.41283950617284;
a[14][18]=0.45827160493827;
a[14][20]=-0.4898765432098747;
b[15]=-0.007340403805885992;
a[15][14]=23.71982867220962;
a[15][15]=-0.03224993701184786;
a[15][16]=33.2013101536911;
a[15][9]=-28.42831947593852;
a[15][8]=0.330561854371377;
a[15][10]=-0.3466868228772978;
a[15][20]=-28.42831947593852;
a[15][19]=0.330561854371377;
a[15][21]=-0.3466868228772978;
b[17]=0.02448355592787266;
a[17][18]=74.27160493827159;
a[17][17]=-13.4320987654321;
a[17][19]=-3.950617283950617;
a[17][12]=-21.72839506172839;
a[17][14]=1.975308641975309;
a[17][13]=-8.691358024691358;
a[17][22]=-21.72839506172839;
a[17][24]=1.975308641975309;
a[17][23]=-8.691358024691358;
b[18]=0.01470631174743176;
a[18][17]=17.47050754458161;
a[18][18]=-0.1755829903978139;
a[18][19]=39.59396433470508;
a[18][13]=-28.35665294924554;
a[18][12]=0.7462277091906734;
a[18][14]=-0.8340192043895763;
a[18][23]=-28.35665294924554;
a[18][22]=0.7462277091906734;
a[18][24]=-0.8340192043895763;
b[19]=0.005728088691830557;
a[19][18]=21.83901234567902;
a[19][19]=-0.06320987654320476;
a[19][20]=35.11308641975308;
a[19][14]=-28.41283950617284;
a[19][13]=0.45827160493827;
a[19][15]=-0.4898765432098747;
a[19][24]=-28.41283950617284;
a[19][23]=0.45827160493827;
a[19][25]=-0.4898765432098747;
b[22]=0.05146667920053005;
a[22][23]=74.27160493827159;
a[22][22]=-13.4320987654321;
a[22][24]=-3.950617283950617;
a[22][17]=-21.72839506172839;
a[22][19]=1.975308641975309;
a[22][18]=-8.691358024691358;
a[22][26]=-23.7037037037037;
a[22][27]=-4.74074074074074;
b[5]=-0.008092587680904055;
a[5][5]=-25.10036289008332;
a[5][4]=23.53691538109786;
a[5][11]=-0;
a[5][10]=1.952315031405406;
a[5][3]=-4.688714472653889;
a[5][9]=-1.172951932263014;
b[10]=-0.006743181000547374;
a[10][10]=-101.6200730047915;
a[10][9]=22.49424958870105;
a[10][4]=28.43468983386678;
a[10][3]=-0.4377259680188596;
a[10][11]=23.85083804178547;
a[10][5]=0.2682517908855347;
a[10][8]=-0.2215884662341047;
a[10][2]=-0.2601515562095363;
a[10][15]=-0.1021026077601177;
b[11]=-0.01394322014469498;
a[11][11]=-21.92340011793079;
a[11][10]=21.58320544438007;
a[11][9]=-1.463800868184856;
a[11][3]=-2.048025439714809;
a[11][4]=-0.9244192450806664;
a[11][16]=0.2440298979792365;
b[16]=0.05231242979406764;
a[16][16]=-22.88186020223699;
a[16][15]=22.31324750748786;
a[16][10]=26.68072232566674;
a[16][9]=-1.056397384272231;
a[16][11]=0.2517049276006288;
a[16][14]=-0.9476519776162851;
a[16][8]=-1.192884645996966;
a[16][2]=-1.438117314377646;
a[16][3]=-1.059078031510487;
b[20]=0.1920435374921324;
a[20][20]=-78.24992372250692;
a[20][19]=21.08377569743498;
a[20][15]=27.98026612697332;
a[20][14]=-0.9453979282060349;
a[20][21]=14.70934439343081;
a[20][16]=0.3466868228772977;
a[20][18]=-0;
a[20][13]=-0.7816187987041225;
a[20][7]=-1.02472678972443;
a[20][8]=-0.8579440648549654;
b[21]=-0.00979688810311187;
a[21][21]=-13.14687114362383;
a[21][20]=13.09070698962783;
a[21][19]=-0.4769940673139393;
a[21][14]=-0.7854956726398815;
a[21][8]=-1.093997277965824;
a[21][15]=-0.3646657593219411;
a[21][9]=-0.6731673646478833;
b[23]=0.03415699077578808;
a[23][23]=-87.87303378547674;
a[23][22]=10.85099903453563;
a[23][18]=28.18730931199816;
a[23][17]=-0.8530820455878043;
a[23][24]=25.38439645856634;
a[23][19]=0.602186266291964;
a[23][12]=-0.2405798615121256;
a[23][13]=-0.3030691623623662;
b[24]=0.4407623145955407;
a[24][24]=-57.11137390558506;
a[24][23]=17.11009793263905;
a[24][25]=11.55683152850157;
a[24][19]=27.4758509148833;
a[24][18]=-1.099800918518962;
a[24][20]=0.4898765432098747;
a[24][22]=-0;
a[24][17]=-0.3460700358718459;
a[24][12]=-1.041972599580332;
a[24][13]=-1.337431877289178;
a[24][14]=-1.632891154998025;
b[25]=0.03273250180037153;
a[25][25]=-9.714610583504175;
a[25][24]=9.714610583504175;
a[25][20]=13.78224099415373;
a[25][19]=-0.5664952782728514;
a[25][21]=0.1707855569868072;
a[25][18]=-0.5851317923644941;
a[25][13]=-1.210959706693501;
a[25][14]=-1.029481149893507;
b[26]=-0.002235136553914946;
a[26][27]=55.20640652879609;
a[26][26]=-83.65085097324054;
a[26][22]=21.85703511361476;
a[26][24]=-1.975308641975309;
a[26][23]=7.223338619802056;
a[26][18]=-2.807398757892236;
a[26][17]=-1.210739301116568;
b[27]=0.02343671046291189;
a[27][27]=-3.390831067623068;
a[27][26]=3.390831067623068;
a[27][23]=22.88051584836229;
a[27][22]=-1.717985452632928;
a[27][24]=1.307381248385104;
a[27][17]=-2.4397530031676;
a[27][18]=-1.571325037323317;


  PetscInt *pdnnz = new PetscInt[28]; // diagonals
  PetscInt *ponnz = new PetscInt[28]; // offdiagonals
  PetscInt *snnz = new PetscInt[28];  // sum of both
  PetscInt snnzrow = 0;
  PetscInt pdnnzrow = 0;
  PetscInt ponnzrow = 0;
  for( int i=0; i<28; i++ ) {
    int n=0;
    for( int j=0; j<28; j++ ) {
      if( a[i][j] ) n++;
    }
pout() << "row " << i << " " << n << endl;
    pdnnz[i] = n;
    ponnz[i] = 0;
    snnz[i] = n;
    if( snnzrow < n ) snnzrow = n;
  }

  PetscInt ierr;
  PetscInt localSize=28;
  PetscInt globalSize=28;
  ierr = MatCreate( a_communicator, &a_A ); CHKERRQ(ierr);
  ierr = MatSetSizes(a_A, localSize, localSize, globalSize, globalSize ); CHKERRQ(ierr);
  ierr = MatSetBlockSize(a_A, 1); CHKERRQ(ierr);
  ierr = MatSetType(a_A, MATAIJ); CHKERRQ(ierr);

  ierr = MatSeqAIJSetPreallocation(a_A, snnzrow, snnz); CHKERRQ(ierr);
  ierr = MatMPIAIJSetPreallocation(a_A, pdnnzrow, pdnnz, ponnzrow, ponnz); CHKERRQ(ierr);
  ierr = MatSetOption(a_A, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE); CHKERRQ(ierr);
  ierr = MatSetFromOptions(a_A); CHKERRQ(ierr);
  ierr = MatZeroEntries(a_A); CHKERRQ(ierr);


  delete[] pdnnz;
  delete[] ponnz;
  delete[] snnz;

  ierr = VecCreate( a_communicator, &a_b ); CHKERRQ(ierr);
  ierr = VecSetFromOptions( a_b ); CHKERRQ(ierr);
  ierr = VecSetSizes( a_b, localSize, globalSize ); CHKERRQ(ierr);
  ierr = VecZeroEntries(a_b); CHKERRQ(ierr);


  for( PetscInt i=0; i<28; i++ ) {
    Real val = b[i];
    ierr = VecSetValues(a_b,1,&i,&val,INSERT_VALUES); CHKERRQ(ierr);
    for( PetscInt j=0; j<28; j++ ) {
      Real val = a[i][j];
      if( val ) {
        ierr = MatSetValues(a_A,1,&i,1,&j,&val,INSERT_VALUES); CHKERRQ(ierr);
      }
    }
  }

  ierr = MatAssemblyBegin(a_A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
  ierr = MatAssemblyEnd(a_A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
  ierr = VecAssemblyBegin(a_b); CHKERRQ(ierr);
  ierr = VecAssemblyEnd(a_b); CHKERRQ(ierr);


}

int main(int argc, char *argv[])
{
  #ifdef CH_MPI
  MPI_Comm a_communicator = Chombo_MPI::comm;
#else
  MPI_Comm a_communicator = PETSC_COMM_SELF;
#endif
  PetscLogStage stageProj;
  int ierr=PetscInitialize(&argc,&argv,"./petscrc",PETSC_NULL);
  PetscLogStageRegister("Projection",&stageProj);
  PetscLogStagePush(stageProj);


  Mat A;
  Vec b;
  load(A,b);
  

  KSP ksp;
  ierr = KSPCreate( a_communicator, &ksp ); CHKERRQ(ierr);
  ierr = KSPSetOptionsPrefix(ksp,"pressure_"); CHKERRQ(ierr);
  ierr = KSPSetFromOptions(ksp); CHKERRQ(ierr);
#if (PETSC_VERSION_MAJOR >3 || (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR >= 5))
  ierr = KSPSetOperators(ksp,A,A); CHKERRQ(ierr);
#else
  ierr = KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN); CHKERRQ(ierr);
#endif



  PC pc;
  ierr = KSPGetPC(ksp,&pc); CHKERRQ(ierr);
  ierr = MatView(A,PETSC_VIEWER_STDOUT_SELF);

  ierr = KSPSolve( ksp, b, b ); CHKERRQ(ierr);

  KSPConvergedReason reason;
  ierr = KSPGetConvergedReason(ksp,&reason); CHKERRQ(ierr);
  pout() << " KSP " << ((reason>=0) ? "converged" : "diverged") << ". reason= " << reason << endl;
}