1: #include <petsc-private/kspimpl.h>  /*I "petscksp.h" I*/

  3: typedef struct {
  4:   PetscViewer viewer;
  5:   PetscInt    neigs;
  6:   PetscReal   *eigi;
  7:   PetscReal   *eigr;
  8: #ifdef PETSC_HAVE_AMS
  9:   AMS_Memory  amem;
 10: #endif
 11: } KSPMonitor_AMS;

 15: /*@C
 16:    KSPMonitorAMSCreate - create an AMS monitor context

 18:    Collective

 20:    Input Arguments:
 21: +  ksp - KSP to monitor
 22: -  amscommname - name of AMS communicator to use

 24:    Output Arguments:
 25: .  ctx - context for monitor

 27:    Level: developer

 29: .seealso: KSPMonitorAMS(), KSPMonitorAMSDestroy()
 30: @*/
 31: PetscErrorCode KSPMonitorAMSCreate(KSP ksp,const char *amscommname,void **ctx)
 32: {
 34:   KSPMonitor_AMS *mon;

 37:   PetscNewLog(ksp,KSPMonitor_AMS,&mon);
 38: #ifdef PETSC_HAVE_AMS
 39:   PetscViewerAMSOpen(((PetscObject)ksp)->comm,amscommname,&mon->viewer);
 40:   mon->amem = -1;
 41: #endif
 42:   *ctx = (void*)mon;
 43:   return(0);
 44: }

 48: /*@C
 49:    KSPMonitorAMSDestroy - destroy a monitor context created with KSPMonitorAMSCreate()

 51:    Collective

 53:    Input Arguments:
 54: .  ctx - monitor context

 56:    Level: developer

 58: .seealso: KSPMonitorAMSCreate()
 59: @*/
 60: PetscErrorCode KSPMonitorAMSDestroy(void **ctx)
 61: {
 62:   KSPMonitor_AMS *mon = (KSPMonitor_AMS*)*ctx;

 66: #ifdef PETSC_HAVE_AMS
 67:   if (mon->amem != -1) {
 68:     AMS_Memory_destroy(mon->amem);
 69:     mon->amem = -1;
 70:   }
 71: #endif
 72:   PetscViewerDestroy(&mon->viewer);
 73:   PetscFree(mon->eigr);
 74:   mon->eigi = PETSC_NULL;
 75:   PetscFree(*ctx);
 76:   return(0);
 77: }

 81: /*@C
 82:    KSPMonitorAMS - monitor solution using AMS

 84:    Logically Collective on KSP

 86:    Input Parameters:
 87: +  ksp   - iterative context
 88: .  n     - iteration number
 89: .  rnorm - 2-norm (preconditioned) residual value (may be estimated).  
 90: -  ctx -  PetscViewer of type AMS

 92:    Level: advanced

 94: .keywords: KSP, CG, monitor, AMS, singular values

 96: .seealso: KSPMonitorSingularValue(), KSPComputeExtremeSingularValues(), PetscViewerAMSOpen()
 97: @*/
 98: PetscErrorCode KSPMonitorAMS(KSP ksp,PetscInt n,PetscReal rnorm,void *ctx)
 99: {
100: #if defined(PETSC_HAVE_AMS)
102:   KSPMonitor_AMS *mon = (KSPMonitor_AMS*)ctx;
103:   PetscViewer viewer = mon->viewer;
104:   PetscReal emax,emin;;
105:   AMS_Comm acomm;

110:   KSPComputeExtremeSingularValues(ksp,&emax,&emin);

112:   /* UnPublish  */
113:   if (mon->amem != -1) {AMS_Memory_destroy(mon->amem);}
114:   mon->amem = -1;

116:   PetscFree(mon->eigr);
117:   PetscMalloc(2*n*sizeof(PetscReal),&mon->eigr);
118:   mon->eigi = mon->eigr + n;
119:   if (n) {KSPComputeEigenvalues(ksp,n,mon->eigr,mon->eigi,&mon->neigs);}

121:   PetscViewerAMSGetAMSComm(viewer,&acomm);
122:   AMS_Memory_create(acomm,"ksp_monitor_ams",&mon->amem);
123:   AMS_Memory_take_access(mon->amem);

125:   AMS_Memory_add_field(mon->amem,"rnorm",&ksp->rnorm,1,AMS_DOUBLE,AMS_READ,AMS_COMMON,AMS_REDUCT_UNDEF);
126:   AMS_Memory_add_field(mon->amem,"neigs",&mon->neigs,1,AMS_INT,AMS_READ,AMS_COMMON,AMS_REDUCT_UNDEF);
127:   if (mon->neigs > 0) {
128:     AMS_Memory_add_field(mon->amem,"eigr",&mon->eigr,mon->neigs,AMS_DOUBLE,AMS_READ,AMS_COMMON,AMS_REDUCT_UNDEF);
129:     AMS_Memory_add_field(mon->amem,"eigi",&mon->eigr,mon->neigs,AMS_DOUBLE,AMS_READ,AMS_COMMON,AMS_REDUCT_UNDEF);
130:   }
131:   AMS_Memory_publish(mon->amem);
132:   AMS_Memory_grant_access(mon->amem);

134:   PetscInfo2(ksp,"KSP extreme singular values min=%G max=%G\n",emin,emax);
135:   return(0);
136: #else
138:   SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"Missing package AMS");
139:   return(0);
140: #endif
141: }