static char help[] = "Reads a PETSc matrix and vector from a socket connection,  solves a linear system and sends the result back.\n";

/*T
   Concepts: KSP^solving a linear system
   Processors: n
T*/

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

/* Define context for user-provided preconditioner */
typedef struct {
  Mat SD,Jsp,Jss;
} GS2PC;
typedef struct {
 /*  Mat SD,Jss; */
 Mat SD,Jss;
} BJ2PC;
typedef struct {
  Mat SD,Jsp,Jss,Jps;
  Vec diag;
} DP2PC;

/* Declare routines for user-provided preconditioner */
extern PetscErrorCode GS2PCCreate(GS2PC**);
extern PetscErrorCode GS2PCSetUp(GS2PC*,Mat,Vec);
extern PetscErrorCode GS2PCApply(void*,Vec x,Vec y);
extern PetscErrorCode GS2PCDestroy(GS2PC*);
extern PetscErrorCode BJ2PCCreate(BJ2PC**);
extern PetscErrorCode BJ2PCSetUp(BJ2PC*,Mat,Vec);
extern PetscErrorCode BJ2PCApply(void*,Vec x,Vec y);
extern PetscErrorCode BJ2PCDestroy(BJ2PC*);
extern PetscErrorCode DP2PCCreate(DP2PC**);
extern PetscErrorCode DP2PCSetUp(DP2PC*,Mat,Vec);
extern PetscErrorCode DP2PCApply(void*,Vec x,Vec y);
extern PetscErrorCode DP2PCDestroy(DP2PC*);
/* extern FILE *f,*f1,*f2; */
/* PetscTruth test,flg4; */
/* char  fname[PETSC_MAX_PATH_LEN]; */



#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **args)
{
  KSP            ksp;             /* linear solver context */
  Mat            A;            /* matrices */
  Vec            x,b,contr;          /* approx solution, RHS, exact solution */
  PetscViewer    fd;               /* viewer */
  PC             pc;
  FILE           *f;
  char           fname[PETSC_MAX_PATH_LEN];
  char    fnamepre[PETSC_MAX_PATH_LEN];
  GS2PC  *shell;    /* user-defined preconditioner context */
  BJ2PC  *shell1;
  DP2PC  *shell2;
  PetscErrorCode ierr;
  PetscReal      max,min; /*Maximun and minimun Singular Values*/
  PetscInt        its,N,i;
  PetscTruth     user_defined_pc,flg2,flg3,flg5;
  PetscTruth     flg4;
 /*  PetscInt       its2; */

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

  ierr = VecLoad(fd,VECSEQ,&b);CHKERRQ(ierr);
  ierr = MatLoad(fd,MATAIJ,&A);CHKERRQ(ierr);
   ierr = VecDuplicate(b,&x);CHKERRQ(ierr);
  ierr = VecGetSize(b,&N);CHKERRQ(ierr);
  ierr = VecCreate(PETSC_COMM_SELF,&contr);CHKERRQ(ierr);
  ierr = VecSetSizes(contr,N,N);CHKERRQ(ierr);
  ierr = VecSetType(contr,VECSEQ);CHKERRQ(ierr);
  ierr = VecSet(contr,1);CHKERRQ(ierr);
  ierr = VecEqual(b,contr,&flg5);CHKERRQ(ierr);
  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);

 


   ierr = KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
  ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
  ierr = PetscOptionsHasName(PETSC_NULL,"-2SGS_pc",&user_defined_pc);CHKERRQ(ierr);
   ierr = PetscOptionsHasName(PETSC_NULL,"-2SBJ_pc",&flg2);CHKERRQ(ierr);
  ierr = PetscOptionsHasName(PETSC_NULL,"-2SDP_pc",&flg3);CHKERRQ(ierr);
  if (user_defined_pc) {

    ierr = PCSetType(pc,PCSHELL);CHKERRQ(ierr);

    ierr = GS2PCCreate(&shell);CHKERRQ(ierr);
    ierr = PCShellSetApply(pc,GS2PCApply);CHKERRQ(ierr);
    ierr = PCShellSetContext(pc,shell);CHKERRQ(ierr);
    ierr = PCShellSetName(pc,"Two Stage Gauss-Seidel");CHKERRQ(ierr);
    ierr = GS2PCSetUp(shell,A,x);CHKERRQ(ierr);

  } 

  else if  (flg2){

    ierr = PCSetType(pc,PCSHELL);CHKERRQ(ierr);
    ierr = BJ2PCCreate(&shell1);CHKERRQ(ierr);
    ierr = PCShellSetApply(pc,BJ2PCApply);CHKERRQ(ierr);
    ierr = PCShellSetContext(pc,shell1);CHKERRQ(ierr);
    ierr = PCShellSetName(pc,"Two Stage Block Jacobi");CHKERRQ(ierr);
    ierr = BJ2PCSetUp(shell1,A,x);CHKERRQ(ierr);


  } 
  else if (flg3){
    ierr = PCSetType(pc,PCSHELL);CHKERRQ(ierr);
    ierr = DP2PCCreate(&shell2);CHKERRQ(ierr);
    ierr = PCShellSetApply(pc,DP2PCApply);CHKERRQ(ierr);
    ierr = PCShellSetContext(pc,shell2);CHKERRQ(ierr);
    ierr = PCShellSetName(pc,"Two Stage Discrete Projection");CHKERRQ(ierr);
    ierr = DP2PCSetUp(shell2,A,x);CHKERRQ(ierr);

  }  else {
    ierr = PCSetType(pc,PCBJACOBI);CHKERRQ(ierr);
  }


  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
  ierr = KSPSetComputeSingularValues(ksp,1);CHKERRQ(ierr);
  ierr = KSPSetUp(ksp);CHKERRQ(ierr);
   ierr = KSPGMRESSetRestart(ksp,30);CHKERRQ(ierr);
  ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
                        PetscOptionsGetString(PETSC_NULL,"-f",fname,PETSC_MAX_PATH_LEN-1,&flg4);
                        if (!flg4) SETERRQ(1,"Must indicate binary file with the -f option");
                        ierr = KSPComputeExtremeSingularValues(ksp,&max,&min);CHKERRQ(ierr);
                        ierr = KSPGetIterationNumber(ksp,&its);
			sprintf(fnamepre, "Out_it_");
			strcat(fnamepre,fname);
			f = fopen(fnamepre, "a+");
			ierr = PetscFPrintf(PETSC_COMM_WORLD,f,"%D %g %g \n",its,max,min);CHKERRQ(ierr);
			fclose(f);
            


  ierr = VecView(x,fd);CHKERRQ(ierr);
  ierr = MatDestroy(A);CHKERRQ(ierr); 
  ierr = VecDestroy(b);CHKERRQ(ierr);
  ierr = VecDestroy(x);CHKERRQ(ierr);


  ierr = KSPDestroy(ksp);CHKERRQ(ierr); 

  if (user_defined_pc) {
    ierr = GS2PCDestroy(shell);CHKERRQ(ierr);
  }
  if (flg2) {
    ierr = BJ2PCDestroy(shell1);CHKERRQ(ierr);
  }
  if (flg3) {
    ierr = DP2PCDestroy(shell2);CHKERRQ(ierr);
  }

  for (i=0;!flg5;i++){
                     ierr = VecLoad(fd,VECSEQ,&b);CHKERRQ(ierr);
		     ierr = VecEqual(b,contr,&flg5);
		       if (!flg5) {
  ierr = MatLoad(fd,MATAIJ,&A);CHKERRQ(ierr);
  ierr = VecDuplicate(b,&x);CHKERRQ(ierr);
  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);

   ierr = KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
  ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);

  ierr = PetscOptionsHasName(PETSC_NULL,"-2SGS_pc",&user_defined_pc);CHKERRQ(ierr);
   ierr = PetscOptionsHasName(PETSC_NULL,"-2SBJ_pc",&flg2);CHKERRQ(ierr);
  ierr = PetscOptionsHasName(PETSC_NULL,"-2SDP_pc",&flg3);CHKERRQ(ierr);
  if (user_defined_pc) {

    ierr = PCSetType(pc,PCSHELL);CHKERRQ(ierr);

 
    ierr = GS2PCCreate(&shell);CHKERRQ(ierr);


    ierr = PCShellSetApply(pc,GS2PCApply);CHKERRQ(ierr);
    ierr = PCShellSetContext(pc,shell);CHKERRQ(ierr);


    ierr = PCShellSetName(pc,"Two Stage Gauss-Seidel");CHKERRQ(ierr);

 
    ierr = GS2PCSetUp(shell,A,x);CHKERRQ(ierr);

  } 

  else if  (flg2){

    ierr = PCSetType(pc,PCSHELL);CHKERRQ(ierr);
    ierr = BJ2PCCreate(&shell1);CHKERRQ(ierr);
    ierr = PCShellSetApply(pc,BJ2PCApply);CHKERRQ(ierr);
    ierr = PCShellSetContext(pc,shell1);CHKERRQ(ierr);
    ierr = PCShellSetName(pc,"Two Stage Block Jacobi");CHKERRQ(ierr);
    ierr = BJ2PCSetUp(shell1,A,x);CHKERRQ(ierr);


  } 
  else if (flg3){
    ierr = PCSetType(pc,PCSHELL);CHKERRQ(ierr);
    ierr = DP2PCCreate(&shell2);CHKERRQ(ierr);
    ierr = PCShellSetApply(pc,DP2PCApply);CHKERRQ(ierr);
    ierr = PCShellSetContext(pc,shell2);CHKERRQ(ierr);
    ierr = PCShellSetName(pc,"Two Stage Discrete Projection");CHKERRQ(ierr);
    ierr = DP2PCSetUp(shell2,A,x);CHKERRQ(ierr);

  }  else {
    ierr = PCSetType(pc,PCBJACOBI);CHKERRQ(ierr);
  }


  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
  ierr = KSPSetComputeSingularValues(ksp,1);CHKERRQ(ierr);
  ierr = KSPSetUp(ksp);CHKERRQ(ierr);
   ierr = KSPGMRESSetRestart(ksp,30);CHKERRQ(ierr);
  ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
                        PetscOptionsGetString(PETSC_NULL,"-f",fname,PETSC_MAX_PATH_LEN-1,&flg4);
                        if (!flg4) SETERRQ(1,"Must indicate binary file with the -f option");
                        ierr = KSPComputeExtremeSingularValues(ksp,&max,&min);CHKERRQ(ierr);
                        ierr = KSPGetIterationNumber(ksp,&its);
			sprintf(fnamepre, "Out_it_");
			strcat(fnamepre,fname);
			f = fopen(fnamepre, "a+");
			ierr = PetscFPrintf(PETSC_COMM_WORLD,f,"%D %g %g \n",its,max,min);CHKERRQ(ierr);
			fclose(f);
            



  ierr = VecView(x,fd);CHKERRQ(ierr);
  ierr = MatDestroy(A);CHKERRQ(ierr); 
  ierr = VecDestroy(b);CHKERRQ(ierr);
  ierr = VecDestroy(x);CHKERRQ(ierr);


  ierr = KSPDestroy(ksp);CHKERRQ(ierr); 

  if (user_defined_pc) {
    ierr = GS2PCDestroy(shell);CHKERRQ(ierr);
  }
  if (flg2) {
    ierr = BJ2PCDestroy(shell1);CHKERRQ(ierr);
  }
  if (flg3) {
    ierr = DP2PCDestroy(shell2);CHKERRQ(ierr);
  }



         }


}



  ierr = VecDestroy(contr);CHKERRQ(ierr);


  ierr = PetscFinalize();CHKERRQ(ierr);
  return 0;
}


/***********************************************************************/
/*          Routines for a user-defined shell preconditioner           */
/***********************************************************************/

#undef __FUNCT__  
#define __FUNCT__ "GS2PCCreate"

PetscErrorCode GS2PCCreate(GS2PC **shell)
{
  GS2PC  *newctx;
  PetscErrorCode ierr;

  ierr         = PetscNew(GS2PC,&newctx);CHKERRQ(ierr);
  *shell       = newctx;
  return 0;
}
/* ------------------------------------------------------------------- */
#undef __FUNCT__  
#define __FUNCT__ "GS2PCSetUp"

PetscErrorCode GS2PCSetUp(GS2PC *shell,Mat pmat,Vec x)
{
 
  Mat            *submat;
  PetscErrorCode ierr;
  PetscInt       N;
  IS             set1[4],set2[4]; 
  


  ierr = VecGetSize(x,&N);CHKERRQ(ierr); 
  ierr = ISCreateStride(PETSC_COMM_SELF,N/2,0,1,&set1[0]);CHKERRQ(ierr);
  ierr = ISCreateStride(PETSC_COMM_SELF,N/2,0,1,&set2[0]);CHKERRQ(ierr);
  ierr = ISCreateStride(PETSC_COMM_SELF,N/2,N/2,1,&set1[1]);CHKERRQ(ierr);
  ierr = ISCreateStride(PETSC_COMM_SELF,N/2,0,1,&set2[1]);CHKERRQ(ierr);
  ierr = ISCreateStride(PETSC_COMM_SELF,N/2,N/2,1,&set1[2]);CHKERRQ(ierr);
  ierr = ISCreateStride(PETSC_COMM_SELF,N/2,N/2,1,&set2[2]);CHKERRQ(ierr);
  ierr = MatGetSubMatrices(pmat,3,set1,set2,MAT_INITIAL_MATRIX,&submat);CHKERRQ(ierr);
  shell->SD =  submat[0];
  shell->Jsp = submat[1];
  shell->Jss = submat[2];

  ierr = ISDestroy(set1[0]);CHKERRQ(ierr);
  ierr = ISDestroy(set2[0]);CHKERRQ(ierr);
  ierr = ISDestroy(set1[1]);CHKERRQ(ierr);
  ierr = ISDestroy(set2[1]);CHKERRQ(ierr);
  ierr = ISDestroy(set1[2]);CHKERRQ(ierr);
  ierr = ISDestroy(set2[2]);CHKERRQ(ierr);

  return 0;
}
/* ------------------------------------------------------------------- */
#undef __FUNCT__  
#define __FUNCT__ "GS2PCApply"

PetscErrorCode GS2PCApply(void *ctx,Vec x,Vec y)
{
  GS2PC   *shell = (GS2PC*)ctx;
  PetscErrorCode  ierr;
  KSP             ksp2,ksp3;
  PC              pc2,pc3;
  PetscInt       N,i,pos,posR;
  Vec             Rn,Rw,p,s,y1; 
  PetscReal      *xx,value1,value2;
  PetscTruth   flg4;
  FILE *f2,*f3;
  char fnamepre[PETSC_MAX_PATH_LEN];
  char fname[PETSC_MAX_PATH_LEN];
  PetscInt its;
  PetscReal   max,min;
  ierr = VecGetSize(x,&N);CHKERRQ(ierr);
  ierr = VecCreate(PETSC_COMM_WORLD,&Rn);CHKERRQ(ierr);
  ierr = VecSetSizes(Rn,PETSC_DECIDE,N/2);CHKERRQ(ierr);
  ierr = VecSetFromOptions(Rn);CHKERRQ(ierr);
  ierr = VecDuplicate(Rn,&Rw);CHKERRQ(ierr);
  ierr = VecDuplicate(Rn,&p);CHKERRQ(ierr);
  ierr = VecDuplicate(Rn,&s);CHKERRQ(ierr);
  ierr = VecDuplicate(Rn,&y1);CHKERRQ(ierr);
  ierr = VecGetArray(x,&xx);CHKERRQ(ierr);
  for (i=0;i<N/2;i++){
    pos=N-i-1;
    posR=pos-N/2;
    value1=xx[i];
    value2=xx[pos];
    VecSetValues(Rn,1,&i,&value1,INSERT_VALUES);
    VecSetValues(Rw,1,&posR,&value2,INSERT_VALUES);
  }
  ierr = VecAssemblyBegin(Rn);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(Rn);CHKERRQ(ierr);
  ierr = VecAssemblyBegin(Rw);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(Rw);CHKERRQ(ierr);
  ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr);




 
 ierr = KSPCreate(PETSC_COMM_WORLD,&ksp2);CHKERRQ(ierr);
  ierr = KSPSetOperators(ksp2,shell->SD,shell->SD,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
  ierr = KSPGetPC(ksp2,&pc2);CHKERRQ(ierr);
  ierr = PCSetType(pc2,PCILU);CHKERRQ(ierr);
  ierr = KSPSetFromOptions(ksp2);CHKERRQ(ierr);
  ierr = KSPSetComputeSingularValues(ksp2,1);CHKERRQ(ierr);
  ierr = KSPSetUp(ksp2);CHKERRQ(ierr);
  ierr = KSPSolve(ksp2,Rn,p);CHKERRQ(ierr);

                        ierr=PetscOptionsGetString(PETSC_NULL,"-f",fname,PETSC_MAX_PATH_LEN-1,&flg4);CHKERRQ(ierr);
                        if (!flg4) SETERRQ(1,"Must indicate binary file with the -f option");
                        ierr = KSPComputeExtremeSingularValues(ksp2,&max,&min);CHKERRQ(ierr);
                        ierr = KSPGetIterationNumber(ksp2,&its);
		
			sprintf(fnamepre, "Inn_it_pri_");
			strcat(fnamepre,fname);
			f2 = fopen(fnamepre, "a+");
			ierr = PetscFPrintf(PETSC_COMM_WORLD,f2,"%D %g %g \n",its,max,min);CHKERRQ(ierr);
			fclose(f2);
                      





  /********************Step 4:y=Rw-Jcp*p **********************************/
  ierr = VecScale(p,-1);CHKERRQ(ierr);
  ierr = MatMultAdd(shell->Jsp,p,Rw,y1);CHKERRQ(ierr);

  /******* Step 5: First sublinear system solution Jcc*q=Rw *************/
  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp3);CHKERRQ(ierr);
  ierr = KSPSetOperators(ksp3,shell->Jss,shell->Jss,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
  ierr = KSPGetPC(ksp3,&pc3);CHKERRQ(ierr);
  ierr = PCSetType(pc3,PCBJACOBI);CHKERRQ(ierr);
  ierr = KSPSetFromOptions(ksp3);CHKERRQ(ierr);
  ierr = KSPSetComputeSingularValues(ksp3,1);CHKERRQ(ierr);
  ierr = KSPSetUp(ksp3);CHKERRQ(ierr);
  ierr = KSPSolve(ksp3,y1,s);CHKERRQ(ierr); 

                        ierr = KSPComputeExtremeSingularValues(ksp3,&max,&min);CHKERRQ(ierr);
                        ierr = KSPGetIterationNumber(ksp3,&its);
	
			sprintf(fnamepre, "Inn_it_sec_");
			strcat(fnamepre,fname);
			f3 = fopen(fnamepre, "a+");
			ierr = PetscFPrintf(PETSC_COMM_WORLD,f3,"%D %g %g \n",its,max,min);CHKERRQ(ierr);
			fclose(f3);
         


 
  ierr = VecGetArray(p,&xx);CHKERRQ(ierr);
  for (i=0;i<N/2;i++){
    value1=xx[i];
    VecSetValues(y,1,&i,&value1,INSERT_VALUES);
  }
 ierr = VecRestoreArray(p,&xx);CHKERRQ(ierr);

 ierr = VecGetArray(s,&xx);CHKERRQ(ierr);
 /*Ver VecSetValuesBlocked*/
  for (i=0;i<N/2;i++){
    value1=xx[i];
    pos=i+N/2-1;
    VecSetValues(y,1,&pos,&value1,INSERT_VALUES);
  }
 ierr = VecRestoreArray(s,&xx);CHKERRQ(ierr);



  ierr = VecAssemblyBegin(y);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(y);CHKERRQ(ierr);





  ierr = VecDestroy(Rn);CHKERRQ(ierr);
  ierr = VecDestroy(Rw);CHKERRQ(ierr);
  ierr = VecDestroy(p);CHKERRQ(ierr);
  ierr = VecDestroy(s);CHKERRQ(ierr);
  ierr = VecDestroy(y1);CHKERRQ(ierr);
  ierr = KSPDestroy(ksp3);CHKERRQ(ierr);
  ierr = KSPDestroy(ksp2);CHKERRQ(ierr);

  return 0;
}
/* ------------------------------------------------------------------- */
#undef __FUNCT__  
#define __FUNCT__ "GS2PCDestroy"

PetscErrorCode GS2PCDestroy(GS2PC *shell)
{
  PetscErrorCode ierr;

  ierr = MatDestroy(shell->SD);CHKERRQ(ierr);
  ierr = MatDestroy(shell->Jsp);CHKERRQ(ierr);
  ierr = MatDestroy(shell->Jss);CHKERRQ(ierr);
  ierr = PetscFree(shell);CHKERRQ(ierr);
 
  return 0;
}
/*********************************************/
/*                                           */ 
/*           2SBJ Functions                  */
/*                                           */
/*********************************************/


#undef __FUNCT__  
#define __FUNCT__ "BJ2PCCreate"

PetscErrorCode BJ2PCCreate(BJ2PC **shell)
{
  BJ2PC  *newctx;
  PetscErrorCode ierr;

  ierr         = PetscNew(BJ2PC,&newctx);CHKERRQ(ierr);
  *shell       = newctx;
  return 0;
}
/* ------------------------------------------------------------------- */
#undef __FUNCT__  
#define __FUNCT__ "BJ2PCSetUp"

PetscErrorCode BJ2PCSetUp(BJ2PC *shell,Mat pmat,Vec x)
{

  Mat            *submat;
  PetscErrorCode ierr;
  PetscInt       N;
  IS             set1[2],set2[2]; 
  


  ierr = VecGetSize(x,&N);CHKERRQ(ierr); 
  ierr = ISCreateStride(PETSC_COMM_SELF,N/2,0,1,&set1[0]);CHKERRQ(ierr);
  ierr = ISCreateStride(PETSC_COMM_SELF,N/2,0,1,&set2[0]);CHKERRQ(ierr);
  ierr = ISCreateStride(PETSC_COMM_SELF,N/2,N/2,1,&set1[1]);CHKERRQ(ierr);
  ierr = ISCreateStride(PETSC_COMM_SELF,N/2,N/2,1,&set2[1]);CHKERRQ(ierr);

  ierr = MatGetSubMatrices(pmat,2,set1,set2,MAT_INITIAL_MATRIX,&submat);CHKERRQ(ierr);
  shell->SD =  submat[0];
  shell->Jss = submat[1];

  ierr = ISDestroy(set1[0]);CHKERRQ(ierr);
  ierr = ISDestroy(set2[0]);CHKERRQ(ierr);
  ierr = ISDestroy(set1[1]);CHKERRQ(ierr);
  ierr = ISDestroy(set2[1]);CHKERRQ(ierr);
  return 0;
}
/* ------------------------------------------------------------------- */
#undef __FUNCT__  
#define __FUNCT__ "BJ2PCApply"

PetscErrorCode BJ2PCApply(void *ctx,Vec x,Vec y)
{
  BJ2PC   *shell = (BJ2PC*)ctx;
  PetscErrorCode  ierr;
  KSP             ksp2,ksp3;
  PC              pc2,pc3;
  PetscInt       N,i,pos,posR;
  Vec             Rn,Rw,p,s; 
  PetscReal      *xx,value1,value2;
  PetscTruth   flg4;
  FILE *f2,*f3;
  char fnamepre[PETSC_MAX_PATH_LEN];
  char fname[PETSC_MAX_PATH_LEN];
  PetscInt its;
  PetscReal   max,min;


ierr = VecGetSize(x,&N);CHKERRQ(ierr);

   ierr = VecCreate(PETSC_COMM_WORLD,&Rn);CHKERRQ(ierr);
  ierr = VecSetSizes(Rn,PETSC_DECIDE,N/2);CHKERRQ(ierr);
  ierr = VecSetFromOptions(Rn);CHKERRQ(ierr);
  ierr = VecDuplicate(Rn,&Rw);CHKERRQ(ierr);
  ierr = VecDuplicate(Rn,&p);CHKERRQ(ierr);
  ierr = VecDuplicate(Rn,&s);CHKERRQ(ierr);

  ierr = VecGetArray(x,&xx);CHKERRQ(ierr);
  for (i=0;i<N/2;i++){
    pos=N-i-1;
    posR=pos-N/2;
    value1=xx[i];
    value2=xx[pos];
    VecSetValues(Rn,1,&i,&value1,INSERT_VALUES);
    VecSetValues(Rw,1,&posR,&value2,INSERT_VALUES);
  }

  ierr = VecAssemblyBegin(Rn);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(Rn);CHKERRQ(ierr);

  ierr = VecAssemblyBegin(Rw);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(Rw);CHKERRQ(ierr);
  ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr);
 




  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp2);CHKERRQ(ierr);
  ierr = KSPSetOperators(ksp2,shell->SD,shell->SD,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
  ierr = KSPGetPC(ksp2,&pc2);CHKERRQ(ierr);
  ierr = PCSetType(pc2,PCILU);CHKERRQ(ierr);
  ierr = KSPSetFromOptions(ksp2);CHKERRQ(ierr);
  ierr = KSPSetComputeSingularValues(ksp2,1);CHKERRQ(ierr);
  ierr = KSPSetUp(ksp2);CHKERRQ(ierr);
  ierr = KSPSolve(ksp2,Rn,p);CHKERRQ(ierr);

                        ierr=PetscOptionsGetString(PETSC_NULL,"-f",fname,PETSC_MAX_PATH_LEN-1,&flg4);CHKERRQ(ierr);
                        if (!flg4) SETERRQ(1,"Must indicate binary file with the -f option");
                        ierr = KSPComputeExtremeSingularValues(ksp2,&max,&min);CHKERRQ(ierr);
                        ierr = KSPGetIterationNumber(ksp2,&its);
			sprintf(fnamepre, "Inn_it_pri_");
			strcat(fnamepre,fname);
			f2 = fopen(fnamepre, "a+");
			ierr = PetscFPrintf(PETSC_COMM_WORLD,f2,"%D %g %g \n",its,max,min);CHKERRQ(ierr);
			fclose(f2);


  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp3);CHKERRQ(ierr);
  ierr = KSPSetOperators(ksp3,shell->Jss,shell->Jss,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
  ierr = KSPGetPC(ksp3,&pc3);CHKERRQ(ierr);
  ierr = PCSetType(pc3,PCBJACOBI);CHKERRQ(ierr);
  ierr = KSPSetFromOptions(ksp3);CHKERRQ(ierr);
  ierr = KSPSetComputeSingularValues(ksp3,1);CHKERRQ(ierr);
  ierr = KSPSetUp(ksp3);CHKERRQ(ierr);
  ierr = KSPSolve(ksp3,Rw,s);CHKERRQ(ierr); 

                        ierr = KSPComputeExtremeSingularValues(ksp3,&max,&min);CHKERRQ(ierr);
                        ierr = KSPGetIterationNumber(ksp3,&its);
			sprintf(fnamepre, "Inn_it_sec_");
			strcat(fnamepre,fname);
			f3 = fopen(fnamepre, "a+");
			ierr = PetscFPrintf(PETSC_COMM_WORLD,f3,"%D %g %g \n",its,max,min);CHKERRQ(ierr);
			fclose(f3);

  
 
  ierr = VecGetArray(p,&xx);CHKERRQ(ierr);
  for (i=0;i<N/2;i++){
    value1=xx[i];
    VecSetValues(y,1,&i,&value1,INSERT_VALUES);
  }
 ierr = VecRestoreArray(p,&xx);CHKERRQ(ierr);

 ierr = VecGetArray(s,&xx);CHKERRQ(ierr);
  for (i=0;i<N/2;i++){
    value1=xx[i];
    pos=i+N/2-1;
    VecSetValues(y,1,&pos,&value1,INSERT_VALUES);
  }
 ierr = VecRestoreArray(s,&xx);CHKERRQ(ierr);

ierr = VecAssemblyBegin(y);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(y);CHKERRQ(ierr);

  ierr = VecDestroy(Rn);CHKERRQ(ierr);
  ierr = VecDestroy(Rw);CHKERRQ(ierr);
  ierr = VecDestroy(p);CHKERRQ(ierr);
  ierr = VecDestroy(s);CHKERRQ(ierr);
  ierr = KSPDestroy(ksp3);CHKERRQ(ierr);
  ierr = KSPDestroy(ksp2);CHKERRQ(ierr);

  return 0;
}
/* ------------------------------------------------------------------- */
#undef __FUNCT__  
#define __FUNCT__ "BJ2PCDestroy"

PetscErrorCode BJ2PCDestroy(BJ2PC *shell)
{
  PetscErrorCode ierr;

  ierr = MatDestroy(shell->SD);CHKERRQ(ierr);
  ierr = MatDestroy(shell->Jss);CHKERRQ(ierr);
  ierr = PetscFree(shell);CHKERRQ(ierr);
 
  return 0;
}





/*********************************************/
/*                                           */ 
/*           2SDP Functions                  */
/*                                           */
/*********************************************/






#undef __FUNCT__  
#define __FUNCT__ "DP2PCCreate"

PetscErrorCode DP2PCCreate(DP2PC **shell)
{
  DP2PC  *newctx;
  PetscErrorCode ierr;

  ierr         = PetscNew(DP2PC,&newctx);CHKERRQ(ierr);
  *shell       = newctx;
  return 0;
}
/* ------------------------------------------------------------------- */
#undef __FUNCT__  
#define __FUNCT__ "DP2PCSetUp"

PetscErrorCode DP2PCSetUp(DP2PC *shell,Mat pmat,Vec x)
{
  Vec             diag;
  Mat            *submat,A;
  PetscErrorCode ierr;
  PetscInt       N;
  IS             set1[4],set2[4]; 
  


  ierr = VecGetSize(x,&N);CHKERRQ(ierr); 
  ierr = VecCreate(PETSC_COMM_SELF,&diag);CHKERRQ(ierr);
  ierr = VecSetSizes(diag,PETSC_DECIDE,N/2);CHKERRQ(ierr);
  ierr = VecSetFromOptions(diag);CHKERRQ(ierr);
  ierr = ISCreateStride(PETSC_COMM_SELF,N/2,0,1,&set1[0]);CHKERRQ(ierr);
  ierr = ISCreateStride(PETSC_COMM_SELF,N/2,0,1,&set2[0]);CHKERRQ(ierr);
  ierr = ISCreateStride(PETSC_COMM_SELF,N/2,0,1,&set1[1]);CHKERRQ(ierr);
  ierr = ISCreateStride(PETSC_COMM_SELF,N/2,N/2,1,&set2[1]);CHKERRQ(ierr);
  ierr = ISCreateStride(PETSC_COMM_SELF,N/2,N/2,1,&set1[2]);CHKERRQ(ierr);
  ierr = ISCreateStride(PETSC_COMM_SELF,N/2,0,1,&set2[2]);CHKERRQ(ierr);
  ierr = ISCreateStride(PETSC_COMM_SELF,N/2,N/2,1,&set1[3]);CHKERRQ(ierr);
  ierr = ISCreateStride(PETSC_COMM_SELF,N/2,N/2,1,&set2[3]);CHKERRQ(ierr);
  ierr = MatGetSubMatrices(pmat,4,set1,set2,MAT_INITIAL_MATRIX,&submat);CHKERRQ(ierr);
 /*  ierr = MatDuplicate(submat[0],MAT_DO_NOT_COPY_VALUES,&A); */
  ierr = MatGetDiagonal(submat[3],diag);CHKERRQ(ierr);
  ierr = VecReciprocal(diag);CHKERRQ(ierr);
  shell->SD =  submat[3];
  ierr = MatDiagonalScale(shell->SD,diag,PETSC_NULL);CHKERRQ(ierr);
  ierr = MatMatMult(submat[2],shell->SD,MAT_INITIAL_MATRIX,2,&A);CHKERRQ(ierr);
  ierr = MatScale(A,-1);CHKERRQ(ierr);
  ierr = MatAXPY(A,1,submat[0],DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
  shell->SD = A;
  shell->Jps = submat[1];
  shell->Jsp = submat[2];
  shell->Jss = submat[3];
  shell->diag= diag;

  ierr = ISDestroy(set1[0]);CHKERRQ(ierr);
  ierr = ISDestroy(set2[0]);CHKERRQ(ierr);
  ierr = ISDestroy(set1[1]);CHKERRQ(ierr);
  ierr = ISDestroy(set2[1]);CHKERRQ(ierr);
  ierr = ISDestroy(set1[2]);CHKERRQ(ierr);
  ierr = ISDestroy(set2[2]);CHKERRQ(ierr);
  ierr = ISDestroy(set1[3]);CHKERRQ(ierr);
  ierr = ISDestroy(set2[3]);CHKERRQ(ierr);
  ierr = MatDestroy(submat[0]);CHKERRQ(ierr);

  return 0;
}
/* ------------------------------------------------------------------- */
#undef __FUNCT__  
#define __FUNCT__ "DP2PCApply"

PetscErrorCode DP2PCApply(void *ctx,Vec x,Vec y)
{
  DP2PC   *shell = (DP2PC*)ctx;
  PetscErrorCode  ierr;
  KSP             ksp2,ksp3;
  PC              pc2,pc3;
  PetscInt       N,i,pos,posR;
  Vec             Rn,Rw,p,s,y1,w,q; 
  PetscReal      *xx,value1,value2;
  PetscTruth   flg4;
  FILE *f2,*f3;
  char fnamepre[PETSC_MAX_PATH_LEN];
  char fname[PETSC_MAX_PATH_LEN];
  PetscInt its;
  PetscReal   max,min;




ierr = VecGetSize(x,&N);CHKERRQ(ierr);

   ierr = VecCreate(PETSC_COMM_WORLD,&Rn);CHKERRQ(ierr);
  ierr = VecSetSizes(Rn,PETSC_DECIDE,N/2);CHKERRQ(ierr);
  ierr = VecSetFromOptions(Rn);CHKERRQ(ierr);
  ierr = VecDuplicate(Rn,&Rw);CHKERRQ(ierr);
  ierr = VecDuplicate(Rn,&p);CHKERRQ(ierr);
   ierr = VecDuplicate(Rn,&q);CHKERRQ(ierr);
 ierr = VecDuplicate(Rn,&w);CHKERRQ(ierr);
  ierr = VecDuplicate(Rn,&s);CHKERRQ(ierr);
  ierr = VecDuplicate(Rn,&y1);CHKERRQ(ierr);

  ierr = VecGetArray(x,&xx);CHKERRQ(ierr);
  for (i=0;i<N/2;i++){
    pos=N-i-1;
    posR=pos-N/2;
    value1=xx[i];
    value2=xx[pos];
    VecSetValues(Rn,1,&i,&value1,INSERT_VALUES);
    VecSetValues(Rw,1,&posR,&value2,INSERT_VALUES);
  }
  /*Ensamblo los Vectores*/
  ierr = VecAssemblyBegin(Rn);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(Rn);CHKERRQ(ierr);

  ierr = VecAssemblyBegin(Rw);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(Rw);CHKERRQ(ierr);
  /*Devuelvo el arreglo*/
  ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr);
 




/******* Step 1: First sublinear system solution Jcc*q=Rw *************/
 ierr = VecPointwiseMult(q,Rw,shell->diag);CHKERRQ(ierr);
  /********************Step 2:w=Rn-Jpc*q **********************************/
ierr = VecScale(q,-1);CHKERRQ(ierr);
  ierr = MatMultAdd(shell->Jps,q,Rn,w);CHKERRQ(ierr);
  /******* Step 3: First sublinear system solution SD*p=w *************/ 
  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp2);CHKERRQ(ierr);
  ierr = KSPSetOperators(ksp2,shell->SD,shell->SD,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
  ierr = KSPGetPC(ksp2,&pc2);CHKERRQ(ierr);
  ierr = PCSetType(pc2,PCILU);CHKERRQ(ierr);
  ierr = KSPSetFromOptions(ksp2);CHKERRQ(ierr);
  ierr = KSPSetComputeSingularValues(ksp2,1);CHKERRQ(ierr);
  ierr = KSPSetUp(ksp2);CHKERRQ(ierr);
  ierr = KSPSolve(ksp2,w,p);CHKERRQ(ierr);

                        ierr=PetscOptionsGetString(PETSC_NULL,"-f",fname,PETSC_MAX_PATH_LEN-1,&flg4);CHKERRQ(ierr);
                        if (!flg4) SETERRQ(1,"Must indicate binary file with the -f option");
                        ierr = KSPComputeExtremeSingularValues(ksp2,&max,&min);CHKERRQ(ierr);
                        ierr = KSPGetIterationNumber(ksp2,&its);
			/* fnamepre="Out_It_"; */
			sprintf(fnamepre, "Inn_it_pri_");
			strcat(fnamepre,fname);
			f2 = fopen(fnamepre, "a+");
			ierr = PetscFPrintf(PETSC_COMM_WORLD,f2,"%D %g %g \n",its,max,min);CHKERRQ(ierr);
			fclose(f2);


  /********************Step 4:y=Rw-Jcp*p **********************************/
  ierr = VecScale(p,-1);CHKERRQ(ierr);
  ierr = MatMultAdd(shell->Jsp,p,Rw,y1);CHKERRQ(ierr);

  /******* Step 5: First sublinear system solution Jcc*q=Rw *************/
  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp3);CHKERRQ(ierr);
  ierr = KSPSetOperators(ksp3,shell->Jss,shell->Jss,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
  ierr = KSPGetPC(ksp3,&pc3);CHKERRQ(ierr);
  ierr = PCSetType(pc3,PCBJACOBI);CHKERRQ(ierr);
  ierr = KSPSetFromOptions(ksp3);CHKERRQ(ierr);
  ierr = KSPSetComputeSingularValues(ksp3,1);CHKERRQ(ierr);
  ierr = KSPSetUp(ksp3);CHKERRQ(ierr);
  ierr = KSPSolve(ksp3,y1,s);CHKERRQ(ierr); 
                        ierr = KSPComputeExtremeSingularValues(ksp3,&max,&min);CHKERRQ(ierr);
                        ierr = KSPGetIterationNumber(ksp3,&its);
			/* fnamepre="Out_It_"; */
			sprintf(fnamepre, "Inn_it_sec_");
			strcat(fnamepre,fname);
			f3 = fopen(fnamepre, "a+");
			ierr = PetscFPrintf(PETSC_COMM_WORLD,f3,"%D %g %g \n",its,max,min);CHKERRQ(ierr);
			fclose(f3);
  


  /*****************Step 6: Return (p,s)**********************************/
  
 
  ierr = VecGetArray(p,&xx);CHKERRQ(ierr);
  for (i=0;i<N/2;i++){
    value1=xx[i];
    VecSetValues(y,1,&i,&value1,INSERT_VALUES);
  }
 ierr = VecRestoreArray(p,&xx);CHKERRQ(ierr);

 ierr = VecGetArray(s,&xx);CHKERRQ(ierr);
 /*Improve with VecSetValuesBlocked*/
  for (i=0;i<N/2;i++){
    value1=xx[i];
    pos=i+N/2-1;
    VecSetValues(y,1,&pos,&value1,INSERT_VALUES);
  }
 ierr = VecRestoreArray(s,&xx);CHKERRQ(ierr);



  ierr = VecDestroy(Rn);CHKERRQ(ierr);
  ierr = VecDestroy(Rw);CHKERRQ(ierr);
  ierr = VecDestroy(p);CHKERRQ(ierr);
 ierr = VecDestroy(w);CHKERRQ(ierr);
 ierr = VecDestroy(y1);CHKERRQ(ierr);
 ierr = VecDestroy(q);CHKERRQ(ierr);
  ierr = VecDestroy(s);CHKERRQ(ierr);
  ierr = KSPDestroy(ksp3);CHKERRQ(ierr);
  ierr = KSPDestroy(ksp2);CHKERRQ(ierr);
  return 0;
}
/* ------------------------------------------------------------------- */
#undef __FUNCT__  
#define __FUNCT__ "DP2PCDestroy"

PetscErrorCode DP2PCDestroy(DP2PC *shell)
{
  PetscErrorCode ierr;
  ierr = MatDestroy(shell->SD);CHKERRQ(ierr);
  ierr = MatDestroy(shell->Jsp);CHKERRQ(ierr);
  ierr = MatDestroy(shell->Jss);CHKERRQ(ierr);
  ierr = MatDestroy(shell->Jps);CHKERRQ(ierr);
  ierr = VecDestroy(shell->diag);CHKERRQ(ierr);
  ierr = PetscFree(shell);CHKERRQ(ierr);
 
  return 0;
}