Commits

Peter Brune  committed b4a5846 Merge

Merge branch 'prbrune/snes-qnscaling'

  • Participants
  • Parent commits f5cedc2, 6a6e6c9

Comments (0)

Files changed (6)

File include/finclude/petscsnes.h

 ! SNESQNScaleType
 !
 
+      PetscEnum SNES_QN_SCALE_DEFAULT
       PetscEnum SNES_QN_SCALE_NONE
       PetscEnum SNES_QN_SCALE_SHANNO
       PetscEnum SNES_QN_SCALE_LINESEARCH
       PetscEnum SNES_QN_SCALE_JACOBIAN
 
-      parameter(SNES_QN_SCALE_NONE       = 0)
-      parameter(SNES_QN_SCALE_SHANNO     = 1)
-      parameter(SNES_QN_SCALE_LINESEARCH = 2)
-      parameter(SNES_QN_SCALE_JACOBIAN   = 3)
+      parameter(SNES_QN_SCALE_DEFAULT    = 0)
+      parameter(SNES_QN_SCALE_NONE       = 1)
+      parameter(SNES_QN_SCALE_SHANNO     = 2)
+      parameter(SNES_QN_SCALE_LINESEARCH = 3)
+      parameter(SNES_QN_SCALE_JACOBIAN   = 4)
 
 !
 ! SNESQNRestartType
 !
 
+      PetscEnum SNES_QN_RESTART_DEFAULT
       PetscEnum SNES_QN_RESTART_NONE
       PetscEnum SNES_QN_RESTART_POWELL
       PetscEnum SNES_QN_RESTART_PERIODIC
 
-      parameter(SNES_QN_RESTART_NONE     = 0)
-      parameter(SNES_QN_RESTART_POWELL   = 1)
-      parameter(SNES_QN_RESTART_PERIODIC = 2)
+      parameter(SNES_QN_RESTART_DEFAULT  = 0)
+      parameter(SNES_QN_RESTART_NONE     = 1)
+      parameter(SNES_QN_RESTART_POWELL   = 2)
+      parameter(SNES_QN_RESTART_PERIODIC = 3)
 
 !
 ! SNESNCGType

File include/finclude/petscsnesdef.h

 #define SNESLineSearch PetscFortranAddr
 #define SNESLineSearchOrder PetscEnum
 #define SNESNormSchedule PetscEnum
+#define SNESQNType PetscEnum
 #define SNESQNRestartType PetscEnum
 #define SNESQNCompositionType PetscEnum
 #define SNESQNScaleType PetscEnum

File include/petscsnes.h

 
 PETSC_EXTERN PetscErrorCode SNESNCGSetType(SNES, SNESNCGType);
 
-typedef enum {SNES_QN_SCALE_NONE       = 0,
-              SNES_QN_SCALE_SHANNO     = 1,
-              SNES_QN_SCALE_LINESEARCH = 2,
-              SNES_QN_SCALE_JACOBIAN   = 3} SNESQNScaleType;
+typedef enum {SNES_QN_SCALE_DEFAULT    = 0,
+              SNES_QN_SCALE_NONE       = 1,
+              SNES_QN_SCALE_SHANNO     = 2,
+              SNES_QN_SCALE_LINESEARCH = 3,
+              SNES_QN_SCALE_JACOBIAN   = 4} SNESQNScaleType;
 PETSC_EXTERN const char *const SNESQNScaleTypes[];
-typedef enum {SNES_QN_RESTART_NONE     = 0,
-              SNES_QN_RESTART_POWELL   = 1,
-              SNES_QN_RESTART_PERIODIC = 2} SNESQNRestartType;
+typedef enum {SNES_QN_RESTART_DEFAULT  = 0,
+              SNES_QN_RESTART_NONE     = 1,
+              SNES_QN_RESTART_POWELL   = 2,
+              SNES_QN_RESTART_PERIODIC = 3} SNESQNRestartType;
 PETSC_EXTERN const char *const SNESQNRestartTypes[];
+typedef enum {SNES_QN_LBFGS      = 0,
+              SNES_QN_BROYDEN    = 1,
+              SNES_QN_BADBROYDEN = 2
+             } SNESQNType;
+PETSC_EXTERN const char *const SNESQNTypes[];
 
+PETSC_EXTERN PetscErrorCode SNESQNSetType(SNES, SNESQNType);
 PETSC_EXTERN PetscErrorCode SNESQNSetScaleType(SNES, SNESQNScaleType);
 PETSC_EXTERN PetscErrorCode SNESQNSetRestartType(SNES, SNESQNRestartType);
 

File src/snes/examples/tutorials/makefile

 	   else  echo  ${PWD} "\nPossible problem with with ex5_5_qn, diffs above \n========================================="; fi; \
 	   ${RM} -f ex5_5_qn.tmp
 
+runex5_5_broyden:
+	-@${CSD_BASIC_COMMAND_LINE} -snes_type qn -snes_qn_type broyden -snes_qn_m ${N_RESTART} \
+        > ex5_5_broyden.tmp 2>&1; \
+	   if (${DIFF} output/ex5_5_broyden.out ex5_5_broyden.tmp) then true; \
+	   else  echo  ${PWD} "\nPossible problem with with ex5_5_broyden, diffs above \n========================================="; fi; \
+	   ${RM} -f ex5_5_broyden.tmp
+
 runex5_5_ls:
 	-@${CSD_BASIC_COMMAND_LINE} -snes_type newtonls \
         > ex5_5_ls.tmp 2>&1; \
 TESTEXAMPLES_C		       = ex1.PETSc runex1 runex1_2 ex1.rm ex2.PETSc runex2  runex2_3 ex2.rm ex3.PETSc runex3 \
                                  runex3_2 runex3_3 runex3_4 ex3.rm ex4.PETSc runex4 ex4.rm ex5.PETSc runex5 runex5_2 runex5_3 runex5_4 \
                                  runex5_5_ngmres runex5_5_ngmres_nrichardson runex5_5_ncg runex5_5_nrichardson \
-                                 runex5_5_ngmres_ngs runex5_5_qn runex5_5_ls \
+                                 runex5_5_ngmres_ngs runex5_5_qn runex5_5_broyden runex5_5_ls \
                                  runex5_5_fas runex5_5_ngmres_fas runex5_5_fas_additive \
                                  runex5_5_nasm runex5_5_newton_asm_dmda runex5_5_newton_gasm_dmda \
                                  runex5_6 ex5.rm ex7.PETSc runex7 ex7.rm\

File src/snes/examples/tutorials/output/ex5_5_broyden.out

+  0 SNES Function norm 1.11127 
+  1 SNES Function norm 2.11191 
+  2 SNES Function norm 1.14072 
+  3 SNES Function norm 0.428072 
+  4 SNES Function norm 0.349271 
+  5 SNES Function norm 0.515488 
+  6 SNES Function norm 0.312855 
+  7 SNES Function norm 0.362055 
+  8 SNES Function norm 1.03321 
+  9 SNES Function norm 0.45599 
+ 10 SNES Function norm 0.64581 
+ 11 SNES Function norm 0.15351 
+ 12 SNES Function norm 0.139291 
+ 13 SNES Function norm 0.138266 
+ 14 SNES Function norm 2.13085 
+ 15 SNES Function norm 0.30663 
+ 16 SNES Function norm 0.637933 
+ 17 SNES Function norm 0.432727 
+ 18 SNES Function norm 0.176073 
+ 19 SNES Function norm 0.102224 
+ 20 SNES Function norm 0.0906867 
+ 21 SNES Function norm 0.101961 
+ 22 SNES Function norm 0.234768 
+ 23 SNES Function norm 1.12379 
+ 24 SNES Function norm 0.629644 
+ 25 SNES Function norm 0.345374 
+ 26 SNES Function norm 0.253188 
+ 27 SNES Function norm 0.474962 
+ 28 SNES Function norm 0.0592169 
+ 29 SNES Function norm 0.0532744 
+ 30 SNES Function norm 0.0679556 
+ 31 SNES Function norm 0.0616865 
+ 32 SNES Function norm 0.0673367 
+ 33 SNES Function norm 0.193481 
+ 34 SNES Function norm 0.22362 
+ 35 SNES Function norm 0.0831414 
+ 36 SNES Function norm 0.0445144 
+ 37 SNES Function norm 0.042821 
+ 38 SNES Function norm 0.0448081 
+ 39 SNES Function norm 0.0412071 
+ 40 SNES Function norm 0.301366 
+ 41 SNES Function norm 0.0654126 
+ 42 SNES Function norm 0.0367439 
+ 43 SNES Function norm 0.0523366 
+ 44 SNES Function norm 0.036082 
+ 45 SNES Function norm 0.036637 
+ 46 SNES Function norm 0.0315853 
+ 47 SNES Function norm 0.0293608 
+ 48 SNES Function norm 0.0308405 
+ 49 SNES Function norm 0.0300748 
+ 50 SNES Function norm 0.0423938 

File src/snes/impls/qn/qn.c

 
 #define H(i,j)  qn->dXdFmat[i*qn->m + j]
 
-const char *const SNESQNScaleTypes[] =        {"NONE","SHANNO","LINESEARCH","JACOBIAN","SNESQNScaleType","SNES_QN_SCALING_",0};
-const char *const SNESQNRestartTypes[] =      {"NONE","POWELL","PERIODIC","SNESQNRestartType","SNES_QN_RESTART_",0};
-const char *const SNESQNTypes[] =             {"LBFGS","BROYDEN","BADBROYDEN","SNESQNType","SNES_QN_RESTART_",0};
-
-typedef enum {SNES_QN_LBFGS      = 0,
-              SNES_QN_BROYDEN    = 1,
-              SNES_QN_BADBROYDEN = 2
-             } SNESQNType;
+const char *const SNESQNScaleTypes[] =        {"DEFAULT","NONE","SHANNO","LINESEARCH","JACOBIAN","SNESQNScaleType","SNES_QN_SCALING_",0};
+const char *const SNESQNRestartTypes[] =      {"DEFAULT","NONE","POWELL","PERIODIC","SNESQNRestartType","SNES_QN_RESTART_",0};
+const char *const SNESQNTypes[] =             {"LBFGS","BROYDEN","BADBROYDEN","SNESQNType","SNES_QN_",0};
 
 typedef struct {
   Vec               *U;                   /* Stored past states (vary from method to method) */
   Vec               *V;                   /* Stored past states (vary from method to method) */
   PetscInt          m;                    /* The number of kept previous steps */
+  PetscReal         *lambda;              /* The line search history of the method */
+  PetscReal         *norm;                /* norms of the steps */
   PetscScalar       *alpha, *beta;
   PetscScalar       *dXtdF, *dFtdX, *YtdX;
   PetscBool         singlereduction;      /* Aggregated reduction implementation */
 
 #undef __FUNCT__
 #define __FUNCT__ "SNESQNApply_Broyden"
-PetscErrorCode SNESQNApply_Broyden(SNES snes,PetscInt it,Vec Y,Vec X,Vec Xold, Vec D, Vec Dold)
+PetscErrorCode SNESQNApply_Broyden(SNES snes,PetscInt it,Vec Y,Vec X,Vec Xold,Vec D)
 {
   PetscErrorCode     ierr;
   SNES_QN            *qn = (SNES_QN*)snes->data;
   Vec                W   = snes->work[3];
   Vec                *U  = qn->U;
-  Vec                *V  = qn->V;
   KSPConvergedReason kspreason;
-  PetscInt           k,i,lits;
   PetscInt           m = qn->m;
+  PetscInt           k,i,j,lits,l = m;
+  PetscReal          unorm,a,b;
+  PetscReal          *lambda=qn->lambda;
   PetscScalar        gdot;
-  PetscInt           l = m;
+  PetscReal          udot;
 
   PetscFunctionBegin;
   if (it < m) l = it;
+  if (it > 0) {
+    k = (it-1)%l;
+    ierr = SNESLineSearchGetLambda(snes->linesearch,&lambda[k]);CHKERRQ(ierr);
+    ierr = VecCopy(Xold,U[k]);CHKERRQ(ierr);
+    ierr = VecAXPY(U[k],-1.0,X);CHKERRQ(ierr);
+    if (qn->monitor) {
+      ierr = PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr);
+      ierr = PetscViewerASCIIPrintf(qn->monitor, "scaling vector %d of %d by lambda: %14.12e \n",k,l,lambda[k]);CHKERRQ(ierr);
+      ierr = PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr);
+    }
+  }
   if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
     ierr = KSPSolve(snes->ksp,D,W);CHKERRQ(ierr);
     ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr);
   }
 
   /* inward recursion starting at the first update and working forward */
-  if (it > 1) {
-    for (i = 0; i < l-1; i++) {
-      k    = (it+i-l)%l;
-      ierr = VecDot(U[k],Y,&gdot);CHKERRQ(ierr);
-      ierr = VecAXPY(Y,gdot,V[k]);CHKERRQ(ierr);
-      if (qn->monitor) {
-        ierr = PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr);
-        ierr = PetscViewerASCIIPrintf(qn->monitor, "it: %d k: %d gdot: %14.12e\n", it, k, PetscRealPart(gdot));CHKERRQ(ierr);
-        ierr = PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr);
-      }
+  for (i = 0; i < l-1; i++) {
+    j = (it+i-l)%l;
+    k = (it+i-l+1)%l;
+    ierr = VecNorm(U[j],NORM_2,&unorm);CHKERRQ(ierr);
+    ierr = VecDot(U[j],Y,&gdot);CHKERRQ(ierr);
+    unorm *= unorm;
+    udot = PetscRealPart(gdot);
+    a = (lambda[j]/lambda[k]);
+    b = -(1.-lambda[j]);
+    a *= udot/unorm;
+    b *= udot/unorm;
+    ierr = VecAXPBYPCZ(Y,a,b,1.,U[k],U[j]);CHKERRQ(ierr);
+
+    if (qn->monitor) {
+      ierr = PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr);
+      ierr = PetscViewerASCIIPrintf(qn->monitor, "using vector %d and %d, gdot: %14.12e\n",k,j,gdot);CHKERRQ(ierr);
+      ierr = PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr);
     }
   }
-  if (it < m) l = it;
-
-  /* set W to be the last step, post-linesearch */
-  ierr = VecCopy(Xold,W);CHKERRQ(ierr);
-  ierr = VecAXPY(W,-1.0,X);CHKERRQ(ierr);
-
   if (l > 0) {
-    k    = (it-1)%l;
-    ierr = VecCopy(W,U[k]);CHKERRQ(ierr);
-    ierr = VecAXPY(W,-1.0,Y);CHKERRQ(ierr);
-    ierr = VecDot(U[k],W,&gdot);CHKERRQ(ierr);
-    ierr = VecCopy(Y,V[k]);CHKERRQ(ierr);
-    ierr = VecScale(V[k],1.0/gdot);CHKERRQ(ierr);
+    k = (it-1)%l;
     ierr = VecDot(U[k],Y,&gdot);CHKERRQ(ierr);
-    ierr = VecAXPY(Y,gdot,V[k]);CHKERRQ(ierr);
+    ierr = VecNorm(U[k],NORM_2,&unorm);CHKERRQ(ierr);
+    unorm *= unorm;
+    udot = PetscRealPart(gdot);
+    a = unorm/(unorm-lambda[k]*udot);
+    b = -(1.-lambda[k])*udot/(unorm-lambda[k]*udot);
     if (qn->monitor) {
       ierr = PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr);
-      ierr = PetscViewerASCIIPrintf(qn->monitor, "update: %d k: %d gdot: %14.12e\n", it, k, PetscRealPart(gdot));CHKERRQ(ierr);
+      ierr = PetscViewerASCIIPrintf(qn->monitor, "using vector %d: a: %14.12e b: %14.12e \n",k,a,b);CHKERRQ(ierr);
       ierr = PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr);
     }
+    ierr = VecAXPBY(Y,b,a,U[k]);CHKERRQ(ierr);
+  }
+  l = m;
+  if (it+1<m)l=it+1;
+  k = it%l;
+  if (qn->monitor) {
+    ierr = PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr);
+    ierr = PetscViewerASCIIPrintf(qn->monitor, "setting vector %d of %d\n",k,l);CHKERRQ(ierr);
+    ierr = PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr);
   }
   PetscFunctionReturn(0);
 }
 
   /* ksp thing for jacobian scaling */
   KSPConvergedReason kspreason;
-  PetscInt           k, i, lits;
+  PetscInt           h,k,j,i,lits;
   PetscInt           m = qn->m;
-  PetscScalar        gdot;
+  PetscScalar        gdot,udot;
   PetscInt           l = m;
 
   PetscFunctionBegin;
   ierr = VecCopy(D,Y);CHKERRQ(ierr);
   if (l > 0) {
     k    = (it-1)%l;
+    ierr = SNESLineSearchGetLambda(snes->linesearch,&qn->lambda[k]);CHKERRQ(ierr);
     ierr = VecCopy(Dold,U[k]);CHKERRQ(ierr);
     ierr = VecAXPY(U[k],-1.0,D);CHKERRQ(ierr);
-    ierr = VecDot(U[k],U[k],&gdot);CHKERRQ(ierr);
-    ierr = VecCopy(D,T[k]);CHKERRQ(ierr);
-    ierr = VecScale(T[k],1.0/gdot);CHKERRQ(ierr);
-    ierr = VecDot(U[k],Y,&gdot);CHKERRQ(ierr);
-    ierr = VecAXPY(Y,gdot,T[k]);CHKERRQ(ierr);
-    if (qn->monitor) {
-      ierr = PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr);
-      ierr = PetscViewerASCIIPrintf(qn->monitor, "update: %d k: %d gdot: %14.12e\n", it, k, PetscRealPart(gdot));CHKERRQ(ierr);
-      ierr = PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr);
-    }
-  }
-
-  /* inward recursion starting at the first update and working forward */
-  if (it > 1) {
-    for (i = 0; i < l-1; i++) {
-      k    = (it+i-l)%l;
-      ierr = VecDot(U[k],Y,&gdot);CHKERRQ(ierr);
-      ierr = VecAXPY(Y,gdot,T[k]);CHKERRQ(ierr);
-      if (qn->monitor) {
-        ierr = PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr);
-        ierr = PetscViewerASCIIPrintf(qn->monitor, "it: %d k: %d gdot: %14.12e\n", it, k, PetscRealPart(gdot));CHKERRQ(ierr);
-        ierr = PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr);
-      }
-    }
+    ierr = VecCopy(Xold,T[k]);CHKERRQ(ierr);
+    ierr = VecAXPY(T[k],-1.0,X);CHKERRQ(ierr);
   }
 
   if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
   } else {
     ierr = VecScale(Y,qn->scaling);CHKERRQ(ierr);
   }
+
+  /* inward recursion starting at the first update and working forward */
+  if (l > 0) {
+    for (i = 0; i < l-1; i++) {
+      j    = (it+i-l)%l;
+      k    = (it+i-l+1)%l;
+      h    = (it-1)%l;
+      ierr = VecDotBegin(U[j],U[h],&gdot);CHKERRQ(ierr);
+      ierr = VecDotBegin(U[j],U[j],&udot);CHKERRQ(ierr);
+      ierr = VecDotEnd(U[j],U[h],&gdot);CHKERRQ(ierr);
+      ierr = VecDotEnd(U[j],U[j],&udot);CHKERRQ(ierr);
+      ierr = VecAXPY(Y,PetscRealPart(gdot)/PetscRealPart(udot),T[k]);CHKERRQ(ierr);
+      ierr = VecAXPY(Y,-(1.-qn->lambda[k])*PetscRealPart(gdot)/PetscRealPart(udot),T[j]);CHKERRQ(ierr);
+      if (qn->monitor) {
+        ierr = PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr);
+        ierr = PetscViewerASCIIPrintf(qn->monitor, "it: %d k: %d gdot: %14.12e\n", it, k, PetscRealPart(gdot));CHKERRQ(ierr);
+        ierr = PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);CHKERRQ(ierr);
+      }
+    }
+  }
   PetscFunctionReturn(0);
 }
 
     ierr = VecCopy(X, dX[k]);CHKERRQ(ierr);
     ierr = VecAXPY(dX[k], -1.0, Xold);CHKERRQ(ierr);
     if (qn->singlereduction) {
-      PetscScalar dFtdF;
       ierr = VecMDotBegin(dF[k],l,dX,dXtdF);CHKERRQ(ierr);
       ierr = VecMDotBegin(dX[k],l,dF,dFtdX);CHKERRQ(ierr);
       ierr = VecMDotBegin(Y,l,dX,YtdX);CHKERRQ(ierr);
-      if (qn->scale_type == SNES_QN_SCALE_SHANNO) {ierr = VecDotBegin(dF[k],dF[k],&dFtdF);CHKERRQ(ierr);}
       ierr = VecMDotEnd(dF[k],l,dX,dXtdF);CHKERRQ(ierr);
       ierr = VecMDotEnd(dX[k],l,dF,dFtdX);CHKERRQ(ierr);
       ierr = VecMDotEnd(Y,l,dX,YtdX);CHKERRQ(ierr);
-      if (qn->scale_type == SNES_QN_SCALE_SHANNO) {
-        ierr = VecDotEnd(dF[k],dF[k],&dFtdF);CHKERRQ(ierr);
-        qn->scaling = PetscRealPart(dXtdF[k]) / PetscRealPart(dFtdF);
-      }
       for (j = 0; j < l; j++) {
         H(k, j) = dFtdX[j];
         H(j, k) = dXtdF[j];
       for (j = 0; j < l; j++) dXtdF[j] = H(j, j);
     } else {
       ierr = VecDot(dX[k], dF[k], &dXtdF[k]);CHKERRQ(ierr);
-      if (qn->scale_type == SNES_QN_SCALE_SHANNO) {
-        PetscReal dFtdF;
-        ierr        = VecDotRealPart(dF[k],dF[k],&dFtdF);CHKERRQ(ierr);
-        qn->scaling = PetscRealPart(dXtdF[k])/dFtdF;
-      }
     }
     if (qn->scale_type == SNES_QN_SCALE_LINESEARCH) {
       ierr = SNESLineSearchGetLambda(snes->linesearch,&qn->scaling);CHKERRQ(ierr);
   PetscErrorCode      ierr;
   SNES_QN             *qn = (SNES_QN*) snes->data;
   Vec                 X,Xold;
-  Vec                 F;
+  Vec                 F,W;
   Vec                 Y,D,Dold;
   PetscInt            i, i_r;
   PetscReal           fnorm,xnorm,ynorm,gnorm;
   /* basically just a regular newton's method except for the application of the jacobian */
 
   PetscFunctionBegin;
-  F = snes->vec_func;                   /* residual vector */
-  Y = snes->vec_sol_update;             /* search direction generated by J^-1D*/
-
+  F    = snes->vec_func;                /* residual vector */
+  Y    = snes->vec_sol_update;          /* search direction generated by J^-1D*/
+  W    = snes->work[3];
   X    = snes->vec_sol;                 /* solution vector */
   Xold = snes->work[0];
 
   }
 
   for (i = 0, i_r = 0; i < snes->max_its; i++, i_r++) {
+    if (qn->scale_type == SNES_QN_SCALE_SHANNO && i_r > 0) {
+      PetscScalar ff,xf;
+      ierr = VecCopy(Dold,Y);CHKERRQ(ierr);
+      ierr = VecCopy(Xold,W);CHKERRQ(ierr);
+      ierr = VecAXPY(Y,-1.0,D);CHKERRQ(ierr);
+      ierr = VecAXPY(W,-1.0,X);CHKERRQ(ierr);
+      ierr = VecDotBegin(Y,Y,&ff);CHKERRQ(ierr);
+      ierr = VecDotBegin(W,Y,&xf);CHKERRQ(ierr);
+      ierr = VecDotEnd(Y,Y,&ff);CHKERRQ(ierr);
+      ierr = VecDotEnd(W,Y,&xf);CHKERRQ(ierr);
+      qn->scaling = PetscRealPart(xf)/PetscRealPart(ff);
+    }
     switch (qn->type) {
     case SNES_QN_BADBROYDEN:
       ierr = SNESQNApply_BadBroyden(snes,i_r,Y,X,Xold,D,Dold);CHKERRQ(ierr);
       break;
     case SNES_QN_BROYDEN:
-      ierr = SNESQNApply_Broyden(snes,i_r,Y,X,Xold,D,Dold);CHKERRQ(ierr);
+      ierr = SNESQNApply_Broyden(snes,i_r,Y,X,Xold,D);CHKERRQ(ierr);
       break;
     case SNES_QN_LBFGS:
       SNESQNApply_LBFGS(snes,i_r,Y,X,Xold,D,Dold);CHKERRQ(ierr);
     } else {
       ierr = VecCopy(F, D);CHKERRQ(ierr);
     }
-
     powell = PETSC_FALSE;
     if (qn->restart_type == SNES_QN_RESTART_POWELL) {
       /* check restart by Powell's Criterion: |F^T H_0 Fold| > 0.2 * |Fold^T H_0 Fold| */
-      ierr = VecDotBegin(Dold, Dold, &DolddotDold);CHKERRQ(ierr);
-      ierr = VecDotBegin(Dold, D, &DolddotD);CHKERRQ(ierr);
-      ierr = VecDotEnd(Dold, Dold, &DolddotDold);CHKERRQ(ierr);
-      ierr = VecDotEnd(Dold, D, &DolddotD);CHKERRQ(ierr);
+      if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
+        ierr = MatMult(snes->jacobian_pre,Dold,W);CHKERRQ(ierr);
+      } else {
+        ierr = VecCopy(Dold,W);CHKERRQ(ierr);
+      }
+      ierr = VecDotBegin(W, Dold, &DolddotDold);CHKERRQ(ierr);
+      ierr = VecDotBegin(W, D, &DolddotD);CHKERRQ(ierr);
+      ierr = VecDotEnd(W, Dold, &DolddotDold);CHKERRQ(ierr);
+      ierr = VecDotEnd(W, D, &DolddotD);CHKERRQ(ierr);
       if (PetscAbs(PetscRealPart(DolddotD)) > qn->powell_gamma*PetscAbs(PetscRealPart(DolddotDold))) powell = PETSC_TRUE;
     }
     periodic = PETSC_FALSE;
 
   PetscFunctionBegin;
   ierr = VecDuplicateVecs(snes->vec_sol, qn->m, &qn->U);CHKERRQ(ierr);
-  ierr = VecDuplicateVecs(snes->vec_sol, qn->m, &qn->V);CHKERRQ(ierr);
-  ierr = PetscMalloc3(qn->m, PetscScalar, &qn->alpha,
+  if (qn->type != SNES_QN_BROYDEN) ierr = VecDuplicateVecs(snes->vec_sol, qn->m, &qn->V);CHKERRQ(ierr);
+  ierr = PetscMalloc4(qn->m, PetscScalar, &qn->alpha,
                       qn->m, PetscScalar, &qn->beta,
-                      qn->m, PetscScalar, &qn->dXtdF);CHKERRQ(ierr);
+                      qn->m, PetscScalar, &qn->dXtdF,
+                      qn->m, PetscReal, &qn->lambda);CHKERRQ(ierr);
 
   if (qn->singlereduction) {
     ierr = PetscMalloc3(qn->m*qn->m, PetscScalar, &qn->dXdFmat,
                         qn->m, PetscScalar, &qn->YtdX);CHKERRQ(ierr);
   }
   ierr = SNESSetWorkVecs(snes,4);CHKERRQ(ierr);
+  /* set method defaults */
+  if (qn->scale_type == SNES_QN_SCALE_DEFAULT) {
+    if (qn->type == SNES_QN_BADBROYDEN) {
+      qn->scale_type = SNES_QN_SCALE_NONE;
+    } else {
+      qn->scale_type = SNES_QN_SCALE_SHANNO;
+    }
+  }
+  if (qn->restart_type == SNES_QN_RESTART_DEFAULT) {
+    if (qn->type == SNES_QN_LBFGS) {
+      qn->restart_type = SNES_QN_RESTART_POWELL;
+    } else {
+      qn->restart_type = SNES_QN_RESTART_PERIODIC;
+    }
+  }
 
-  /* set up the line search */
   if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
     ierr = SNESSetUpMatrices(snes);CHKERRQ(ierr);
   }
-
   if (snes->pcside == PC_LEFT && snes->functype == SNES_FUNCTION_DEFAULT) {snes->functype = SNES_FUNCTION_UNPRECONDITIONED;}
-
   PetscFunctionReturn(0);
 }
 
     if (qn->singlereduction) {
       ierr = PetscFree3(qn->dXdFmat, qn->dFtdX, qn->YtdX);CHKERRQ(ierr);
     }
-    ierr = PetscFree3(qn->alpha, qn->beta, qn->dXtdF);CHKERRQ(ierr);
+    ierr = PetscFree4(qn->alpha,qn->beta,qn->dXtdF,qn->lambda);CHKERRQ(ierr);
   }
   PetscFunctionReturn(0);
 }
   SNESLineSearch    linesearch;
   SNESQNRestartType rtype = qn->restart_type;
   SNESQNScaleType   stype = qn->scale_type;
+  SNESQNType        qtype = qn->type;
 
   PetscFunctionBegin;
   ierr = PetscOptionsHead("SNES QN options");CHKERRQ(ierr);
   if (flg) ierr = SNESQNSetRestartType(snes,rtype);CHKERRQ(ierr);
 
   ierr = PetscOptionsEnum("-snes_qn_type","Quasi-Newton update type","",SNESQNTypes,
-                          (PetscEnum)qn->type,(PetscEnum*)&qn->type,NULL);CHKERRQ(ierr);
+                          (PetscEnum)qtype,(PetscEnum*)&qtype,&flg);CHKERRQ(ierr);
+  if (flg) {ierr = SNESQNSetType(snes,qtype);CHKERRQ(ierr);}
   ierr = PetscOptionsTail();CHKERRQ(ierr);
   if (!snes->linesearch) {
     ierr = SNESGetLineSearch(snes, &linesearch);CHKERRQ(ierr);
-    ierr = SNESLineSearchSetType(linesearch, SNESLINESEARCHCP);CHKERRQ(ierr);
+    if (qn->type == SNES_QN_LBFGS) {
+      ierr = SNESLineSearchSetType(linesearch, SNESLINESEARCHCP);CHKERRQ(ierr);
+    } else if (qn->type == SNES_QN_BROYDEN) {
+      ierr = SNESLineSearchSetType(linesearch, SNESLINESEARCHBASIC);CHKERRQ(ierr);
+    } else {
+      ierr = SNESLineSearchSetType(linesearch, SNESLINESEARCHL2);CHKERRQ(ierr);
+    }
   }
   if (monflg) {
     qn->monitor = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes));CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 
+#undef __FUNCT__
+#define __FUNCT__ "SNESView_QN"
+static PetscErrorCode SNESView_QN(SNES snes, PetscViewer viewer)
+{
+  SNES_QN        *qn    = (SNES_QN*)snes->data;
+  PetscBool      iascii;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
+  if (iascii) {
+    ierr = PetscViewerASCIIPrintf(viewer,"  QN type is %s, restart type is %s, scale type is %s\n",SNESQNTypes[qn->type],SNESQNRestartTypes[qn->restart_type],SNESQNScaleTypes[qn->scale_type]);CHKERRQ(ierr);
+    ierr = PetscViewerASCIIPrintf(viewer,"  Stored subspace size: %d\n", qn->m);CHKERRQ(ierr);
+    if (qn->singlereduction) {
+      ierr = PetscViewerASCIIPrintf(viewer,"  Using the single reduction variant.\n");CHKERRQ(ierr);
+    }
+  }
+  PetscFunctionReturn(0);
+}
 
 #undef __FUNCT__
 #define __FUNCT__ "SNESQNSetRestartType"
   PetscFunctionReturn(0);
 }
 
+#undef __FUNCT__
+#define __FUNCT__ "SNESQNSetType"
+/*@
+    SNESQNSetType - Sets the quasi-Newton variant to be used in SNESQN.
+
+    Logically Collective on SNES
+
+    Input Parameters:
++   snes - the iterative context
+-   qtype - variant type
+
+    Options Database:
+.   -snes_qn_scale_type<lbfgs,broyden,badbroyden>
+
+    Level: beginner
+
+    SNESQNTypes:
++   SNES_QN_LBFGS - LBFGS variant
+.   SNES_QN_BROYDEN - Broyden variant
+-   SNES_QN_BADBROYDEN - Bad Broyden variant
+
+.keywords: SNES, SNESQN, type, set
+@*/
+
+PetscErrorCode SNESQNSetType(SNES snes, SNESQNType qtype)
+{
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(snes,SNES_CLASSID,1);
+  ierr = PetscTryMethod(snes,"SNESQNSetType_C",(SNES,SNESQNType),(snes,qtype));CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "SNESQNSetType_QN"
+PetscErrorCode SNESQNSetType_QN(SNES snes, SNESQNType qtype)
+{
+  SNES_QN *qn = (SNES_QN*)snes->data;
+
+  PetscFunctionBegin;
+  qn->type = qtype;
+  PetscFunctionReturn(0);
+}
+
 /* -------------------------------------------------------------------------- */
 /*MC
       SNESQN - Limited-Memory Quasi-Newton methods for the solution of nonlinear systems.
 
       References:
 
-      L-Broyden Methods: a generalization of the L-BFGS method to the limited memory Broyden family, M. B. Reed,
-      International Journal of Computer Mathematics, vol. 86, 2009.
+      Kelley, C.T., Iterative Methods for Linear and Nonlinear Equations, Chapter 8, SIAM, 1995.
+
+      R. Byrd, J. Nocedal, R. Schnabel, Representations of Quasi-Newton Matrices and their use in Limited Memory Methods,
+      Technical Report, Northwestern University, June 1992.
 
-      Experiments with Quasi-Newton Methods in Solving Stiff ODE Systems, Peter N. Brown, Alan C. Hindmarsh, Homer F. Walker
-      SIAM J. Sci. Stat. Comput. Vol 6(2), April 1985.
+      Peter N. Brown, Alan C. Hindmarsh, Homer F. Walker, Experiments with Quasi-Newton Methods in Solving Stiff ODE
+      Systems, SIAM J. Sci. Stat. Comput. Vol 6(2), April 1985.
 
 
       Level: beginner
   snes->ops->solve          = SNESSolve_QN;
   snes->ops->destroy        = SNESDestroy_QN;
   snes->ops->setfromoptions = SNESSetFromOptions_QN;
-  snes->ops->view           = 0;
+  snes->ops->view           = SNESView_QN;
   snes->ops->reset          = SNESReset_QN;
 
   snes->pcside = PC_LEFT;
   qn->monitor         = NULL;
   qn->singlereduction = PETSC_TRUE;
   qn->powell_gamma    = 0.9999;
-  qn->scale_type      = SNES_QN_SCALE_SHANNO;
-  qn->restart_type    = SNES_QN_RESTART_POWELL;
+  qn->scale_type      = SNES_QN_SCALE_DEFAULT;
+  qn->restart_type    = SNES_QN_RESTART_DEFAULT;
   qn->type            = SNES_QN_LBFGS;
 
   ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetScaleType_C",SNESQNSetScaleType_QN);CHKERRQ(ierr);
   ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetRestartType_C",SNESQNSetRestartType_QN);CHKERRQ(ierr);
+  ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetType_C",SNESQNSetType_QN);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }