Commits

sarich committed 8042ba1

lnks fails to converge

Comments (0)

Files changed (3)

src/tao/pde_constrained/examples/tutorials/elliptic.c

   ierr = MatShellSetOperation(user->MQ,MATOP_MULT_TRANSPOSE,(void(*)(void))QMatMultTranspose); CHKERRQ(ierr);
 
 
-  ierr = MatCreateShell(PETSC_COMM_WORLD,dsubnlocal*user->ns,PETSC_DECIDE,user->ndata*user->ns,user->nstate,user,&user->MQTQ); CHKERRQ(ierr);
+  ierr = MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,user->nstate,user->nstate,user,&user->MQTQ); CHKERRQ(ierr);
   ierr = MatShellSetOperation(user->MQTQ,MATOP_MULT,(void(*)(void))QTQMatMult); CHKERRQ(ierr);
 
-  ierr = MatCreateShell(PETSC_COMM_WORLD,dsubnlocal*user->ns,PETSC_DECIDE,user->ndata*user->ns,user->nstate,user,&user->LTL); CHKERRQ(ierr);
+  ierr = MatCreateShell(PETSC_COMM_WORLD,ysubnlocal,ysubnlocal,n,n,user,&user->LTL); CHKERRQ(ierr);
   ierr = MatShellSetOperation(user->LTL,MATOP_MULT,(void(*)(void))LTLMatMult); CHKERRQ(ierr);
 
 

src/tao/pde_constrained/impls/lnks/lnks.c

 #include "lnks.h"
 #include "../src/tao/matrix/lmvmmat.h"
 
-#define LNKS_PC_NONE   0
-#define LNKS_PC_P2     1
-#define LNKS_PC_P4     2
-#define LNKS_PC_PETSC  3
-#define LNKS_PC_TYPES  4
-static const char *LNKS_PC[64] = {"none","p2","p4","petsc"};
-static PetscErrorCode LNKSApplyPC(PC,Vec,Vec);
-
-/*
-static PetscErrorCode LNKSComputeLagrangianAndGradient(Tao,Vec,Vec,PetscReal*,Vec,Vec);
-static PetscErrorCode LNKSComputeAugmentedLagrangianAndGradient(Tao,Vec,Vec,PetscReal*,Vec,Vec);
-
-static PetscErrorCode LNKSComputeAugmentedLagrangianAndGradient_LS(TaoLineSearch,Vec,PetscReal*,Vec,void*);
-static PetscErrorCode LNKSScatterLS(TAO_LNKS*,Vec,Vec,Vec);
-static PetscErrorCode LNKSGatherLS(TAO_LNKS*,Vec,Vec,Vec);
- */
+#define LNKS_PC_NONE     0
+#define LNKS_PC_P2       1
+#define LNKS_PC_P4       2
+#define LNKS_PC_P2_LMVM  3
+#define LNKS_PC_P4_LMVM  4
+#define LNKS_PC_PETSC    5
+#define LNKS_PC_TYPES    6
+static const char *LNKS_PC[64] = {"none","p2","p4","p2lmvm","p4lmvm","petsc"};
+static PetscErrorCode LNKSApplyPC_P2(PC,Vec,Vec);
+//static PetscErrorCode LNKSApplyPC_P4(PC,Vec,Vec);
 static PetscErrorCode LNKSScatterBig(TAO_LNKS*,Vec,Vec,Vec,Vec);
 static PetscErrorCode LNKSGatherBig(TAO_LNKS*,Vec,Vec,Vec,Vec);
 static PetscErrorCode LNKSInitializeKKT(Tao);
 static PetscErrorCode LNKSUpdateKKTandRHS(Tao);
 static PetscErrorCode LNKSMatMult(Mat,Vec,Vec);
+static PetscErrorCode LNKSWMatMult(Mat,Vec,Vec);
 
 #undef __FUNCT__
 #define __FUNCT__ "LNKSMatMult"
   ierr = MatShellGetContext(K,&ptr);CHKERRQ(ierr);
   tao = (Tao)ptr;
   lnksP = (TAO_LNKS*)(tao->data);
-  /* TODO: zero solution, check existence of matrices, check for symmetry */
+  ierr = VecSet(lnksP->work_state,0.0);
+  ierr = VecSet(lnksP->work_state2,0.0);
+  ierr = VecSet(lnksP->work_design,0.0);
+  ierr = VecSet(lnksP->work_design2,0.0);
+  ierr = VecSet(lnksP->work_lamda,0.0);
+  ierr = VecSet(lnksP->ys,0.0);
+  ierr = VecSet(lnksP->yd,0.0);
+  ierr = VecSet(lnksP->yl,0.0);
+
   ierr = LNKSScatterBig(lnksP,x,lnksP->rhs_state,lnksP->rhs_design,lnksP->rhs_lamda);CHKERRQ(ierr);
-  ierr = MatMult(tao->Hss,lnksP->rhs_state,lnksP->work_state);CHKERRQ(ierr);
-  ierr = MatMult(tao->Hsd,lnksP->rhs_design,lnksP->work_state2);CHKERRQ(ierr);
-  ierr = MatMultTranspose(tao->jacobian_equality_state,lnksP->rhs_lamda,lnksP->ys);CHKERRQ(ierr);
+  if (tao->Hss) {
+    ierr = MatMult(tao->Hss,lnksP->rhs_state,lnksP->work_state);CHKERRQ(ierr);
+  }
+  if (tao->Hsd) {
+    ierr = MatMult(tao->Hsd,lnksP->rhs_design,lnksP->work_state2);CHKERRQ(ierr);
+    ierr = MatMultTranspose(tao->Hsd,lnksP->rhs_state,lnksP->work_design);CHKERRQ(ierr);
+  }
+  if (tao->Hdd) {
+    ierr = MatMult(tao->Hdd,lnksP->rhs_design,lnksP->work_design2);CHKERRQ(ierr);
+  }
+
+
+  if (tao->jacobian_equality_state) {
+    ierr = MatMultTranspose(tao->jacobian_equality_state,lnksP->rhs_lamda,lnksP->ys);CHKERRQ(ierr);
+    ierr = MatMult(tao->jacobian_equality_state,lnksP->rhs_state,lnksP->work_lamda);CHKERRQ(ierr);
+  }
+  if (tao->jacobian_equality_design) {
+    ierr = MatMultTranspose(tao->jacobian_equality_design,lnksP->rhs_lamda,lnksP->yd);CHKERRQ(ierr);
+    ierr = MatMult(tao->jacobian_equality_design,lnksP->rhs_design,lnksP->yl);CHKERRQ(ierr);
+  }
+
   ierr = VecAXPY(lnksP->ys,1.0,lnksP->work_state);CHKERRQ(ierr);
   ierr = VecAXPY(lnksP->ys,1.0,lnksP->work_state2);CHKERRQ(ierr);
 
-  ierr = MatMultTranspose(tao->Hsd,lnksP->rhs_state,lnksP->work_design);CHKERRQ(ierr);
-  ierr = MatMult(tao->Hdd,lnksP->rhs_design,lnksP->work_design2);CHKERRQ(ierr);
-  ierr = MatMultTranspose(tao->jacobian_equality_design,lnksP->rhs_lamda,lnksP->yd);CHKERRQ(ierr);
   ierr = VecAXPY(lnksP->yd,1.0,lnksP->work_design);CHKERRQ(ierr);
   ierr = VecAXPY(lnksP->yd,1.0,lnksP->work_design2);CHKERRQ(ierr);
 
-  ierr = MatMult(tao->jacobian_equality_state,lnksP->rhs_state,lnksP->work_lamda);CHKERRQ(ierr);
-  ierr = MatMult(tao->jacobian_equality_design,lnksP->rhs_design,lnksP->yl);CHKERRQ(ierr);
   ierr = VecAXPY(lnksP->yl,1.0,lnksP->work_lamda);CHKERRQ(ierr);
 
   ierr = LNKSGatherBig(lnksP,lnksP->bigstep,lnksP->ys,lnksP->yd,lnksP->yl);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 
+#undef __FUNCT__
+#define __FUNCT__ "LNKSWMatMult"
+/* apply reduced space y=  Wz *x_d */
+static PetscErrorCode LNKSWMatMult(Mat W, Vec x, Vec y) {
+  void *ptr;
+  Tao tao;
+  TAO_LNKS *lnksP;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  /* Wz = Ad^T As^-T Hss As^-1 Ad - Ad^T As^-T Hsd - Hsd^T As^-1 Ad + Hdd */
+  /* This used inside of LNKSPC, so don't use work_design* vectors */
+  ierr = MatShellGetContext(W,&ptr);CHKERRQ(ierr);
+  tao = (Tao)ptr;
+  lnksP = (TAO_LNKS*)(tao->data);
+  if (tao->Hdd) {
+    ierr = MatMult(tao->Hdd,x,lnksP->wd1);CHKERRQ(ierr);
+  } else {
+    ierr = VecSet(lnksP->wd1,0.0);CHKERRQ(ierr);
+  }
+
+  if (tao->Hsd) {
+    /*TODO */
+  }
+  if (tao->Hss && tao->jacobian_equality_state && tao->jacobian_equality_design) {
+    ierr = MatMult(tao->jacobian_equality_design,x,lnksP->wl1);CHKERRQ(ierr);
+    if (tao->jacobian_equality_state_inv) {
+      ierr = MatMult(tao->jacobian_equality_state_inv,lnksP->wl1,lnksP->ws1);CHKERRQ(ierr);
+      ierr = MatMult(tao->Hss,lnksP->ws1,lnksP->ws2);CHKERRQ(ierr);
+      ierr = MatMultTranspose(tao->jacobian_equality_state_inv,lnksP->ws2,lnksP->wl1);CHKERRQ(ierr);
+
+    } else {
+      ierr = KSPSolve(lnksP->as_ksp,lnksP->wl1,lnksP->ws1);CHKERRQ(ierr);
+      ierr = MatMult(tao->Hss,lnksP->ws1,lnksP->ws2);CHKERRQ(ierr);
+      ierr = KSPSolveTranspose(lnksP->as_ksp,lnksP->ws2,lnksP->wl1);CHKERRQ(ierr);
+    }
+
+    ierr = MatMultTranspose(tao->jacobian_equality_design,lnksP->wl1,y);CHKERRQ(ierr);
+  } else {
+    ierr = VecSet(y,0.0);CHKERRQ(ierr);
+  }
+
+  ierr = VecAXPY(y,1.0,lnksP->wd1);CHKERRQ(ierr);
+  /*TODO add middle terms */
+
+  PetscFunctionReturn(0);
+}
 
 #undef __FUNCT__
 #define __FUNCT__ "TaoDestroy_LNKS"
     ierr = VecDestroy(&lnksP->work_state);CHKERRQ(ierr);
     ierr = VecDestroy(&lnksP->work_state2);CHKERRQ(ierr);
     ierr = VecDestroy(&lnksP->rhs_state);CHKERRQ(ierr);
+    ierr = VecDestroy(&lnksP->ws1);CHKERRQ(ierr);
+    ierr = VecDestroy(&lnksP->ws2);CHKERRQ(ierr);
     ierr = VecDestroy(&lnksP->work_design);CHKERRQ(ierr);
     ierr = VecDestroy(&lnksP->work_design2);CHKERRQ(ierr);
     ierr = VecDestroy(&lnksP->work_design3);CHKERRQ(ierr);
+    ierr = VecDestroy(&lnksP->wd1);CHKERRQ(ierr);
+    ierr = VecDestroy(&lnksP->wd2);CHKERRQ(ierr);
     ierr = VecDestroy(&lnksP->rhs_design);CHKERRQ(ierr);
     ierr = VecDestroy(&lnksP->work_lamda);CHKERRQ(ierr);
+    ierr = VecDestroy(&lnksP->work_lamda2);CHKERRQ(ierr);
+    ierr = VecDestroy(&lnksP->wl1);CHKERRQ(ierr);
+    ierr = VecDestroy(&lnksP->wl2);CHKERRQ(ierr);
     ierr = VecDestroy(&lnksP->rhs_lamda);CHKERRQ(ierr);
     ierr = VecDestroy(&lnksP->ys);CHKERRQ(ierr);
     ierr = VecDestroy(&lnksP->yd);CHKERRQ(ierr);
     ierr = VecScatterDestroy(&lnksP->state_scatter_big);CHKERRQ(ierr);
     ierr = VecScatterDestroy(&lnksP->design_scatter_big);CHKERRQ(ierr);
     ierr = VecScatterDestroy(&lnksP->constraint_scatter_big);CHKERRQ(ierr);
-    /*    ierr = VecDestroy(&lnksP->X); CHKERRQ(ierr);
-     ierr = VecDestroy(&lnksP->G); CHKERRQ(ierr);
-     ierr = VecDestroy(&lnksP->S); CHKERRQ(ierr);*/
     ierr = KSPDestroy(&lnksP->as_ksp);CHKERRQ(ierr);
+    ierr = KSPDestroy(&lnksP->wz_ksp);CHKERRQ(ierr);
   }
   ierr = PetscFree(tao->data);
   tao->data = PETSC_NULL;
   ierr = PetscOptionsEList("-tao_lnks_pc_type","pc type", "", LNKS_PC, LNKS_PC_TYPES, LNKS_PC[lnksP->pc_type], &lnksP->pc_type,0);CHKERRQ(ierr);
 
   ierr = PetscOptionsTail(); CHKERRQ(ierr);
-  ierr = TaoLineSearchSetFromOptions(tao->linesearch); CHKERRQ(ierr);
+  ierr = KSPSetFromOptions(tao->ksp); CHKERRQ(ierr);
+  ierr = KSPSetFromOptions(lnksP->as_ksp); CHKERRQ(ierr);
 
   PetscFunctionReturn(0);
 }
     if (lnksP->kkt_ismatrixfree) {
       ierr = PetscViewerASCIIPrintf(viewer, "lnks matrix free\n");CHKERRQ(ierr);
     }
-    ierr = PetscViewerASCIIPrintf(viewer," lnks preconditioner %d\n",LNKS_PC[lnksP->pc_type]);
+    ierr = PetscViewerASCIIPrintf(viewer," lnks preconditioner %s\n",LNKS_PC[lnksP->pc_type]);
     ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
   }
   PetscFunctionReturn(0);
   ierr = VecDuplicate(tao->solution, &lnksP->work_state2); CHKERRQ(ierr);
   ierr = VecDuplicate(tao->solution, &lnksP->rhs_state); CHKERRQ(ierr);
   ierr = VecDuplicate(tao->solution, &lnksP->ys); CHKERRQ(ierr);
+  ierr = VecDuplicate(tao->solution, &lnksP->ws1); CHKERRQ(ierr);
+  ierr = VecDuplicate(tao->solution, &lnksP->ws2); CHKERRQ(ierr);
 
   ierr = VecDuplicate(tao->design_solution, &tao->design_gradient); CHKERRQ(ierr);
   ierr = VecDuplicate(tao->design_solution, &tao->design_stepdirection); CHKERRQ(ierr);
   ierr = VecDuplicate(tao->design_solution, &lnksP->work_design3); CHKERRQ(ierr);
   ierr = VecDuplicate(tao->design_solution, &lnksP->rhs_design); CHKERRQ(ierr);
   ierr = VecDuplicate(tao->design_solution, &lnksP->yd); CHKERRQ(ierr);
+  ierr = VecDuplicate(tao->design_solution, &lnksP->wd1); CHKERRQ(ierr);
+  ierr = VecDuplicate(tao->design_solution, &lnksP->wd2); CHKERRQ(ierr);
 
   ierr = VecDuplicate(tao->constraints_equality, &lnksP->lamda); CHKERRQ(ierr);
   ierr = VecDuplicate(tao->constraints_equality, &lnksP->lamda_stepdirection); CHKERRQ(ierr);
   ierr = VecDuplicate(tao->constraints_equality, &lnksP->work_lamda); CHKERRQ(ierr);
+  ierr = VecDuplicate(tao->constraints_equality, &lnksP->work_lamda2); CHKERRQ(ierr);
   ierr = VecDuplicate(tao->constraints_equality, &lnksP->rhs_lamda); CHKERRQ(ierr);
   ierr = VecDuplicate(tao->constraints_equality, &lnksP->yl); CHKERRQ(ierr);
+  ierr = VecDuplicate(tao->constraints_equality, &lnksP->wl1); CHKERRQ(ierr);
+  ierr = VecDuplicate(tao->constraints_equality, &lnksP->wl2); CHKERRQ(ierr);
 
   ierr = VecSet(lnksP->lamda,0.0); CHKERRQ(ierr);
 
   lnksP->bigsize = lnksP->nstateplusdesign+lnksP->nconstraints;
 
 
-  /* TODO duplicate work vectors here */
-
-  /* Create full space vectors for kkt solve */
-  /*
-  ierr = VecCreate(((PetscObject)tao)->comm,&lnksP->X); CHKERRQ(ierr);
-
-
-  ierr = VecGetOwnershipRange(tao->solution,&ulo,&uhi); CHKERRQ(ierr);
-  ierr = VecGetOwnershipRange(tao->design_solution,&vlo,&vhi); CHKERRQ(ierr);
-  ierr = VecGetOwnershipRange(tao->constraints_equality,&clo,&chi);CHKERRQ(ierr);
-  
-  ierr = VecSetSizes(lnksP->X,vhi-vlo + uhi-ulo + chi-clo, lnksP->nstateplusdesign); CHKERRQ(ierr);
-  ierr = VecSetType(lnksP->X,((PetscObject)(tao->solution))->type_name); CHKERRQ(ierr);
-  ierr = VecSetFromOptions(lnksP->X); CHKERRQ(ierr);
-  ierr = VecDuplicate(lnksP->X,&lnksP->G); CHKERRQ(ierr);
-  ierr = VecDuplicate(lnksP->X,&lnksP->S); CHKERRQ(ierr);
-
-  ierr = ISCreateStride(PETSC_COMM_SELF,uhi-ulo,ulo,1,&is_state); CHKERRQ(ierr);
-  ierr = ISCreateStride(PETSC_COMM_SELF,uhi-ulo,ulo+vlo,1,&is_state_in_full); CHKERRQ(ierr);
-  ierr = ISCreateStride(PETSC_COMM_SELF,vhi-vlo,vlo,1,&is_design); CHKERRQ(ierr);
-  ierr = ISCreateStride(PETSC_COMM_SELF,vhi-vlo,uhi+vlo,1,&is_design_in_full); CHKERRQ(ierr);
-
-  ierr = VecScatterCreate(lnksP->X,is_state_in_full,tao->solution,is_state,&lnksP->state_scatter_ls); CHKERRQ(ierr);
-  ierr = VecScatterCreate(lnksP->X,is_design_in_full,tao->design_solution,is_design,&lnksP->design_scatter_ls); CHKERRQ(ierr);
-   ierr = VecScatterCreate(lnksP->X,is_lamda_in_full,tao->constraints_equality,is_lamda,&lnksP->lamda_scatter); CHKERRQ(ierr);
-  ierr = ISDestroy(&is_state); CHKERRQ(ierr);
-  ierr = ISDestroy(&is_design); CHKERRQ(ierr);
-  ierr = ISDestroy(&is_lamda); CHKERRQ(ierr);
-  ierr = ISDestroy(&is_state_in_full); CHKERRQ(ierr);
-  ierr = ISDestroy(&is_design_in_full); CHKERRQ(ierr);
-   ierr = ISDestroy(&is_lamda_in_full); CHKERRQ(ierr);
-   */
-
 
+  lnksP->use_lmvm=PETSC_FALSE;
   switch(lnksP->pc_type) {
   case LNKS_PC_NONE:
     ierr = KSPGetPC(tao->ksp,&lnksP->pc);CHKERRQ(ierr);
     ierr = KSPGetPC(tao->ksp,&lnksP->pc);CHKERRQ(ierr);
     ierr = PCSetFromOptions(lnksP->pc);CHKERRQ(ierr);
     break;
+  case LNKS_PC_P2_LMVM:
+    lnksP->use_lmvm=PETSC_TRUE;
   case LNKS_PC_P2:
-  case LNKS_PC_P4:
     ierr = PCCreate(((PetscObject)tao)->comm,&lnksP->pc); CHKERRQ(ierr);
     ierr = PCSetType(lnksP->pc,PCSHELL);CHKERRQ(ierr);
-    ierr = PCShellSetApply(lnksP->pc,LNKSApplyPC);CHKERRQ(ierr);
+    ierr = PCShellSetApply(lnksP->pc,LNKSApplyPC_P2);CHKERRQ(ierr);
     ierr = PCShellSetName(lnksP->pc,"lnks_pc");CHKERRQ(ierr);
     ierr = PCShellSetContext(lnksP->pc,(void*)tao);CHKERRQ(ierr);
     ierr = KSPSetPC(tao->ksp,lnksP->pc);
+    break;
+  case LNKS_PC_P4_LMVM:
+    lnksP->use_lmvm=PETSC_FALSE;
+  case LNKS_PC_P4:
+    SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Not implemented yet");
   }
 
   /* Bz^-1 is approximated with an L-BFGS system */
   ierr = VecGetLocalSize(tao->design_solution,&nlocal); CHKERRQ(ierr);
-  ierr = MatCreateLMVM(((PetscObject)tao)->comm,nlocal,lnksP->ndesign,&lnksP->Bz); CHKERRQ(ierr);
-  ierr = MatLMVMAllocateVectors(lnksP->Bz,tao->design_solution); CHKERRQ(ierr);
+  if (lnksP->use_lmvm) {
+    ierr = MatCreateLMVM(((PetscObject)tao)->comm,nlocal,lnksP->ndesign,&lnksP->Bz); CHKERRQ(ierr);
+    ierr = MatLMVMAllocateVectors(lnksP->Bz,tao->design_solution); CHKERRQ(ierr);
+  } else {
+    /* TODO set correct local size (same as Hdd?) */
+    ierr = MatCreateShell(((PetscObject)tao)->comm,PETSC_DECIDE,PETSC_DECIDE,lnksP->ndesign,lnksP->ndesign,(void*)tao,&lnksP->Wz);CHKERRQ(ierr);
+    ierr = MatShellSetOperation(lnksP->Wz,MATOP_MULT,(void(*)(void))LNKSWMatMult);CHKERRQ(ierr);
+  }
+
 
   ierr = LNKSInitializeKKT(tao);CHKERRQ(ierr);
 
   TAO_LNKS *lnksP = (TAO_LNKS*)tao->data;
   PetscInt iter=0,its;
   TaoConvergedReason reason = TAO_CONTINUE_ITERATING;
-  //  TaoLineSearchConvergedReason ls_reason = TAOLINESEARCH_CONTINUE_ITERATING;
   PetscReal step=1.0,f;
   PetscReal cnorm, snorm, dnorm, gnorm;
   PetscErrorCode ierr;
     ierr = TaoComputeOptControlHessian(tao,tao->solution,tao->design_solution,tao->Hss,tao->Hss_pre,tao->Hsd,tao->Hsd_pre,tao->Hdd,tao->Hdd_pre);CHKERRQ(ierr);
     ierr = TaoComputeOptControlJacobian(tao,tao->solution,tao->design_solution,tao->jacobian_equality_state,tao->jacobian_equality_state_pre, tao->jacobian_equality_state_inv, tao->jacobian_equality_design,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
 
-    //ierr = VecCopy(lnksP->lamda,lnksP->lamda0); CHKERRQ(ierr);
-    //    ierr = LNKSComputeAugmentedLagrangianAndGradient(tao,tao->solution,tao->design_solution,&lnksP->aug,lnksP->GAugL_U,lnksP->GAugL_V); CHKERRQ(ierr);
-
 
     /* Evaluate constraint norm */
     ierr = VecNorm(tao->constraints_equality, NORM_2, &cnorm); CHKERRQ(ierr);
     ierr = VecNorm(tao->gradient,NORM_2,&snorm); CHKERRQ(ierr);
     ierr = VecNorm(tao->design_gradient,NORM_2,&dnorm); CHKERRQ(ierr);
     gnorm = PetscSqrtReal(snorm*snorm+dnorm*dnorm);
-    //    ierr = VecNorm(lnksP->GAugL_U, NORM_2, &mnorm); CHKERRQ(ierr);
-    //    ierr = VecNorm(lnksP->GAugL_V, NORM_2, &cnorm); CHKERRQ(ierr);
-    //    mnorm = PetscSqrtReal(cnorm*cnorm + mnorm*mnorm);
 
 
     /* Monitor convergence */
 PetscErrorCode TaoCreate_LNKS(Tao tao)
 {
   TAO_LNKS *lnksP;
-  PC pc;
   PetscErrorCode ierr;
   //const char *morethuente_type = TAOLINESEARCH_MT;
   PetscFunctionBegin;
   tao->gatol=1e-4;
   tao->grtol=1e-4;
   lnksP->KKT=0;
+  lnksP->pc_type=LNKS_PC_NONE;
   lnksP->as_symmetric=PETSC_FALSE;
   lnksP->kkt_ismatrixfree=PETSC_FALSE;
-  //ierr = TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch); CHKERRQ(ierr);
-  //ierr = TaoLineSearchSetType(tao->linesearch, morethuente_type); CHKERRQ(ierr);
-
-  //  ierr = TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch,LNKSComputeAugmentedLagrangianAndGradient_LS, tao); CHKERRQ(ierr);
   ierr = KSPCreate(((PetscObject)tao)->comm,&tao->ksp); CHKERRQ(ierr);
-  ierr = KSPGetPC(tao->ksp,&pc);CHKERRQ(ierr);
-  ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr);
-
   ierr = KSPCreate(((PetscObject)tao)->comm,&lnksP->as_ksp); CHKERRQ(ierr);
-  ierr = KSPSetFromOptions(tao->ksp); CHKERRQ(ierr);
-  ierr = KSPSetFromOptions(lnksP->as_ksp); CHKERRQ(ierr);
+  ierr = KSPCreate(((PetscObject)tao)->comm,&lnksP->wz_ksp); CHKERRQ(ierr);
+
 
   PetscFunctionReturn(0);
 
 EXTERN_C_END
 
 
-/* This function is used as a wrapper for
-         LNKSComputeAugmentedLagrangianAndGradient()
-   inside the line search. It extracts current U,V from X and calls the regular routine
- */
-/*
-#undef __FUNCT__
-#define __FUNCT__ "LNKSComputeAugmentedLagrangianAndGradient_LS"
-static PetscErrorCode LNKSComputeAugmentedLagrangianAndGradient_LS(TaoLineSearch ls, Vec X, PetscReal *f, Vec GAugL, void *ptr)
-{
-  Tao tao = (Tao)ptr;
-  TAO_LNKS *lnksP = (TAO_LNKS*)tao->data;
-  PetscErrorCode ierr;
-
-  PetscFunctionBegin;
-  ierr = LNKSScatterLS(lnksP,X,tao->solution,tao->design_solution); CHKERRQ(ierr);
-
-
-  ierr = LNKSComputeAugmentedLagrangianAndGradient(tao,tao->solution,tao->design_solution,f,lnksP->GAugL_U,lnksP->GAugL_V); CHKERRQ(ierr);
 
-  ierr = LNKSGatherLS(lnksP,lnksP->GAugL_U,lnksP->GAugL_V,GAugL); CHKERRQ(ierr);
-
-  PetscFunctionReturn(0);
-}
- */
-/*
-#undef __FUNCT__
-#define __FUNCT__ "LNKSScatterLS"
-PetscErrorCode LNKSScatterLS(TAO_LNKS *lnksP, Vec x, Vec s, Vec d)
-{
-  PetscErrorCode ierr;
-  PetscFunctionBegin;
-  ierr = VecScatterBegin(lnksP->state_scatter_ls,x,u,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
-  ierr = VecScatterEnd(lnksP->state_scatter_ls,x,u,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
-  ierr = VecScatterBegin(lnksP->design_scatter_ls,x,v,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
-  ierr = VecScatterEnd(lnksP->design_scatter_ls,x,v,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
 
-  PetscFunctionReturn(0);
-}
 #undef __FUNCT__
-#define __FUNCT__ "LNKSGatherLS"
-PetscErrorCode LNKSGatherLS(TAO_LNKS *lnksP, Vec u, Vec v, Vec x)
+#define __FUNCT__ "LNKSScatterBig"
+PetscErrorCode LNKSScatterBig(TAO_LNKS *lnksP, Vec big, Vec state, Vec design, Vec constraint)
 {
   PetscErrorCode ierr;
   PetscFunctionBegin;
-  ierr = VecScatterBegin(lnksP->state_scatter_ls,u,x,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
-  ierr = VecScatterEnd(lnksP->state_scatter_ls,u,x,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
-  ierr = VecScatterBegin(lnksP->design_scatter_ls,v,x,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
-  ierr = VecScatterEnd(lnksP->design_scatter_ls,v,x,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
+  ierr = VecScatterBegin(lnksP->state_scatter_big,big,state,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
+  ierr = VecScatterEnd(lnksP->state_scatter_big,big,state,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
+  ierr = VecScatterBegin(lnksP->design_scatter_big,big,design,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
+  ierr = VecScatterEnd(lnksP->design_scatter_big,big,design,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
+  ierr = VecScatterBegin(lnksP->constraint_scatter_big,big,constraint,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
+  ierr = VecScatterEnd(lnksP->constraint_scatter_big,big,constraint,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
 
   PetscFunctionReturn(0);
 }
- */
-
 
 #undef __FUNCT__
-#define __FUNCT__ "LNKSScatterBig"
-PetscErrorCode LNKSScatterBig(TAO_LNKS *lnksP, Vec big, Vec state, Vec design, Vec constraint)
+#define __FUNCT__ "LNKSGatherBig"
+PetscErrorCode LNKSGatherBig(TAO_LNKS *lnksP, Vec big, Vec state, Vec design, Vec constraint)
 {
   PetscErrorCode ierr;
   PetscFunctionBegin;
   PetscFunctionReturn(0);
 }
 
-#undef __FUNCT__
-#define __FUNCT__ "LNKSGatherBig"
-PetscErrorCode LNKSGatherBig(TAO_LNKS *lnksP, Vec big, Vec state, Vec design, Vec constraint)
-{
-  PetscErrorCode ierr;
-  PetscFunctionBegin;
-  ierr = VecScatterBegin(lnksP->state_scatter_big,big,state,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
-  ierr = VecScatterEnd(lnksP->state_scatter_big,big,state,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
-  ierr = VecScatterBegin(lnksP->design_scatter_big,big,design,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
-  ierr = VecScatterEnd(lnksP->design_scatter_big,big,design,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
-  ierr = VecScatterBegin(lnksP->constraint_scatter_big,big,constraint,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
-  ierr = VecScatterEnd(lnksP->constraint_scatter_big,big,constraint,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
-
-  PetscFunctionReturn(0);
-}
-
 
 #undef __FUNCT__
 #define __FUNCT__ "LNKSInitializeKKT"
                    [      c         ]   [g_d + Ad'*lamda ]
                                         [       c        ]   */
   if (tao->jacobian_equality_state) {
-    ierr = MatMultTranspose(tao->jacobian_equality_state,lnksP->lamda,lnksP->work_state);CHKERRQ(ierr);
+    ierr = MatMultTranspose(tao->jacobian_equality_state,lnksP->lamda,lnksP->rhs_state);CHKERRQ(ierr);
   } else {
-    ierr = VecSet(lnksP->work_state,0.0);CHKERRQ(ierr);
+    ierr = VecSet(lnksP->rhs_state,0.0);CHKERRQ(ierr);
   }
-  ierr = VecAXPY(lnksP->work_state,1.0,tao->gradient);CHKERRQ(ierr);
+  ierr = VecAXPY(lnksP->rhs_state,1.0,tao->gradient);CHKERRQ(ierr);
 
   if (tao->jacobian_equality_design) {
-    ierr = MatMultTranspose(tao->jacobian_equality_design,lnksP->lamda,lnksP->work_design);CHKERRQ(ierr);
+    ierr = MatMultTranspose(tao->jacobian_equality_design,lnksP->lamda,lnksP->rhs_design);CHKERRQ(ierr);
   } else {
-    ierr = VecSet(lnksP->work_design,0.0);CHKERRQ(ierr);
+    ierr = VecSet(lnksP->rhs_design,0.0);CHKERRQ(ierr);
   }
-  ierr = VecAXPY(lnksP->work_design,1.0,tao->design_gradient);CHKERRQ(ierr);
-  ierr = LNKSGatherBig(lnksP,lnksP->bigrhs,lnksP->work_state,lnksP->work_design,tao->constraints_equality);CHKERRQ(ierr);
+  ierr = VecAXPY(lnksP->rhs_design,1.0,tao->design_gradient);CHKERRQ(ierr);
+  ierr = VecCopy(tao->constraints_equality,lnksP->rhs_lamda);CHKERRQ(ierr);
+  PetscScalar snorm,dnorm,lnorm,bnorm,gsnorm,gdnorm;
+  if (lnksP->verbose) {
+    VecNorm(tao->gradient,NORM_2,&gsnorm);
+    VecNorm(tao->design_gradient,NORM_2,&gdnorm);
+    VecNorm(lnksP->rhs_state,NORM_2,&snorm);
+    VecNorm(lnksP->rhs_design,NORM_2,&dnorm);
+    VecNorm(lnksP->rhs_lamda,NORM_2,&lnorm);
+    VecNorm(lnksP->bigrhs,NORM_2,&bnorm);
+    PetscPrintf(PETSC_COMM_WORLD,"||gs||=%g, ||gd||=%g, ||rs||=%g, ||rd|=%g, ||rl||=%g, ||rhs||=%g\n",gsnorm,gdnorm,snorm,dnorm,lnorm,bnorm);
+  }
+  ierr = LNKSGatherBig(lnksP,lnksP->bigrhs,lnksP->rhs_state,lnksP->rhs_design,lnksP->rhs_lamda);CHKERRQ(ierr);
   ierr = VecScale(lnksP->bigrhs,-1.0);CHKERRQ(ierr);
-
+  if (lnksP->verbose) {
+    VecNorm(lnksP->rhs_state,NORM_2,&snorm);
+    VecNorm(lnksP->rhs_design,NORM_2,&dnorm);
+    VecNorm(lnksP->rhs_lamda,NORM_2,&lnorm);
+    VecNorm(lnksP->bigrhs,NORM_2,&bnorm);
+    PetscPrintf(PETSC_COMM_WORLD,"||rs||=%g, ||rd|=%g, ||rl||=%g, ||rhs||=%g\n",snorm,dnorm,lnorm,bnorm);
+  }
   PetscFunctionReturn(0);
 }
 
 /*
  P2 :
- P2^-1 =/ As^-   -As^-1 Ad Bz^-1   0     \ /   0           0         I \
+ P2^-1 =/ As^-1   -As^-1 Ad Wz^-1   0     \ /   0           0         I \
+        | 0           Wz^-1        0     | | -Ad^T As^-T   I         0 |
+        \ 0            0           As^-T / \   I           0         0 /
+
+ Wz = Ad^T As^-T Hss As^-T Ad - Ad^T As^-T Hsd - Hsd^T As^-1 Ad + Hdd
+
+ P2_LMVM :
+ P2^-1 =/ As^-1   -As^-1 Ad Bz^-1   0     \ /   0           0         I \
         | 0           Bz^-1        0     | | -Ad^T As^-T   I         0 |
         \ 0            0           As^-T / \   I           0         0 /
 
  Bz = lmvm approximation to Wz
+
+
  Need to solve As^-1 x1 (twice)
                As^-1 x2 (twice)
                Bz^-1 x3 (this uses lmvm)
  */
 #undef __FUNCT__
-#define __FUNCT__ "LNKSApplyPC"
-PetscErrorCode LNKSApplyPC(PC pc, Vec x, Vec y)
+#define __FUNCT__ "LNKSApplyPC_P2"
+PetscErrorCode LNKSApplyPC_P2(PC pc, Vec x, Vec y)
 {
   TAO_LNKS *lnksP;
   Tao      tao;
   ierr = PCShellGetContext(pc,(void**)&tao);CHKERRQ(ierr);
   lnksP = (TAO_LNKS*)(tao->data);
 
-  /* apply (wl)   /   0           0         I \    (Vs)
-           (wd) = | -Ad^T As^-T   I         0 |    (Vd)
-           (ws)   \   I           0         0 /    (Vl)
+  /* apply (wl)    /   0           0         I \    (Vs)
+           (wd3) = | -Ad^T As^-T   I         0 |    (Vd)
+           (ws)    \   I           0         0 /    (Vl)
   */
 
   /* apply permutation */
   ierr = LNKSScatterBig(lnksP,x,lnksP->work_lamda,lnksP->work_design,lnksP->work_state);CHKERRQ(ierr);
-  /* solve As^T wd = Vs */
+  /* Solve As^-T Vs = wl  (or As^T * wl = Vs ) */
+  ierr = KSPSetOperators(lnksP->as_ksp,tao->jacobian_equality_state,tao->jacobian_equality_state);CHKERRQ(ierr);
   if (lnksP->as_symmetric) {
-    ierr = KSPSolve(lnksP->as_ksp,lnksP->work_lamda,lnksP->work_design2);CHKERRQ(ierr);
+    ierr = KSPSolve(lnksP->as_ksp,lnksP->work_state,lnksP->work_lamda2);CHKERRQ(ierr);
   } else {
-    ierr = KSPSolveTranspose(lnksP->as_ksp,lnksP->work_lamda,lnksP->work_design2);CHKERRQ(ierr);
+    ierr = KSPSolveTranspose(lnksP->as_ksp,lnksP->work_state,lnksP->work_lamda2);CHKERRQ(ierr);
   }
   ierr = KSPGetIterationNumber(lnksP->as_ksp,&its); CHKERRQ(ierr);
   lnksP->as_ksp_its += its;
-  /* apply Ad^T wd */
-  ierr = MatMultTranspose(tao->jacobian_equality_design,lnksP->work_design2,lnksP->work_design3);CHKERRQ(ierr);
+  /* apply wd3 = -Ad^T wl2 */
+  ierr = MatMultTranspose(tao->jacobian_equality_design,lnksP->work_lamda2,lnksP->work_design3);CHKERRQ(ierr);
   ierr = VecAXPY(lnksP->work_design,-1.0,lnksP->work_design3);CHKERRQ(ierr);
 
 
 
   /* apply
    ys   / As^-1  -As^-1 Ad Bz^-1   0     \ (wl)
-   yd = | 0           Bz^-1        0     | (wd)
+   yd = | 0           Bz^-1        0     | (wd3)
    yl   \ 0            0           As^-T / (ws)
    */
   if (lnksP->as_symmetric) {
   }
   ierr = KSPGetIterationNumber(lnksP->as_ksp,&its); CHKERRQ(ierr);
   lnksP->as_ksp_its += its;
-  ierr = MatLMVMUpdate(lnksP->Bz,tao->design_solution,lnksP->work_design);
-  ierr = MatLMVMSolve(lnksP->Bz,lnksP->work_design,lnksP->yd); CHKERRQ(ierr);
+  if (lnksP->use_lmvm) {
+    ierr = MatLMVMUpdate(lnksP->Bz,tao->design_solution,lnksP->work_design);
+    ierr = MatLMVMSolve(lnksP->Bz,lnksP->work_design,lnksP->yd); CHKERRQ(ierr);
+  } else {
+    ierr = KSPSetOperators(lnksP->wz_ksp,lnksP->Wz,lnksP->Wz);CHKERRQ(ierr);
+    ierr = KSPSolve(lnksP->wz_ksp,lnksP->work_design,lnksP->yd);CHKERRQ(ierr);
+  }
   ierr = MatMult(tao->jacobian_equality_design,lnksP->yd,lnksP->work_state);CHKERRQ(ierr);
   ierr = KSPSolve(lnksP->as_ksp,lnksP->work_state,lnksP->work_state2);CHKERRQ(ierr);
   ierr = KSPGetIterationNumber(tao->ksp,&its); CHKERRQ(ierr);

src/tao/pde_constrained/impls/lnks/lnks.h

 
   PetscBool kkt_ismatrixfree;
   Mat KKT;
+  Mat Wz;
   Mat Bz; /* quasi-newton approximation for Wz */
 
   //  Vec X,S,G; /* full-space vectors for line search */
   Vec work_design3;
   Vec rhs_design;
   Vec work_lamda;
+  Vec work_lamda2;
   Vec rhs_lamda;
   Vec ys,yd,yl;
+  Vec wd1,wd2,ws1,ws2,wl1,wl2;
   Vec lamda;   /* Lagrange Multiplier */
   Vec lamda_stepdirection;
   Vec bigrhs;
   VecScatter state_scatter_big,design_scatter_big,constraint_scatter_big;
   PetscBool verbose;
   KSP as_ksp;
+  KSP wz_ksp;
   PC  pc;
   PetscInt pc_type;
   PetscInt as_ksp_its;
-  PetscBool as_symmetric;
+  PetscBool as_symmetric,use_lmvm;
 } TAO_LNKS;