Actual source code: petscprom.c

  1: /*  $Id: petsc_prom.c,v 1.7 2005/04/07 17:43:33 adams Exp $ */
  2: /*  Author:             Mark F. Adams    */
  3: /*  Copyright (c) 2004 by Mark F. Adams  */
  4: /*  Filename:           petsc_prom.c     */
  5: 
  6: /*  -------------------------------------------------------------------- 
  7:     
  8:      This file implements a Prometheus preconditioner for matrices that use
  9:      the Mat interface (various matrix formats).  This wraps the Prometheus
 10:      class - this is a C intercace to a C++ code.

 12:      Prometheus assumes that 'PetscScalar' is 'double'.  Prometheus does 
 13:      have a complex-valued solver, but this is runtime parameter, not a 
 14:      compile time parameter.

 16:      The following basic routines are required for each preconditioner.
 17:           PCCreate_XXX()          - Creates a preconditioner context
 18:           PCSetFromOptions_XXX()  - Sets runtime options
 19:           PCApply_XXX()           - Applies the preconditioner
 20:           PCDestroy_XXX()         - Destroys the preconditioner context
 21:      where the suffix "_XXX" denotes a particular implementation, in
 22:      this case we use _Prometheus (e.g., PCCreate_Prometheus, PCApply_Prometheus).
 23:      These routines are actually called via the common user interface
 24:      routines PCCreate(), PCSetFromOptions(), PCApply(), and PCDestroy(), 
 25:      so the application code interface remains identical for all 
 26:      preconditioners.  
 27:  
 28:      Another key routine is:
 29:           PCSetUp_XXX()           - Prepares for the use of a preconditioner
 30:      by setting data structures and options.   The interface routine PCSetUp()
 31:      is not usually called directly by the user, but instead is called by
 32:      PCApply() if necessary.
 33:  
 34:      Additional basic routines are:
 35:           PCView_XXX()            - Prints details of runtime options that
 36:           have actually been used.
 37:           These are called by application codes via the interface routines
 38:           PCView().
 39:  
 40:      The various types of solvers (preconditioners, Krylov subspace methods,
 41:      nonlinear solvers, timesteppers) are all organized similarly, so the
 42:      above description applies to these categories also.  One exception is
 43:      that the analogues of PCApply() for these components are KSPSolve(), 
 44:      SNESSolve(), and TSSolve().

 46:      Additional optional functionality unique to preconditioners is left and
 47:      right symmetric preconditioner application via PCApplySymmetricLeft() 
 48:      and PCApplySymmetricRight().  The Prometheus implementation is 
 49:      PCApplySymmetricLeftOrRight_Prometheus().

 51:     -------------------------------------------------------------------- */

 53:  #include private/pcimpl.h
 54: #include "petscpromproto.h"

 56: /* -------------------------------------------------------------------------- */
 57: /*
 58:    PCCreate_Prometheus - Creates a Prometheus preconditioner context, Prometheus, 
 59:    and sets this as the private data within the generic preconditioning 
 60:    context, PC, that was created within PCCreate().

 62:    Input Parameter:
 63: .  pc - the preconditioner context

 65:    Application Interface Routine: PCCreate()
 66: */

 68: /*MC
 69:      PCPROMETHEUS - Prometheus (i.e. diagonal scaling preconditioning)

 71:    Options Database Key:
 72: .    -pc_prometheus_rowmax - use the maximum absolute value in each row as the scaling factor,
 73:                         rather than the diagonal

 75:    Level: beginner

 77:   Concepts: Prometheus, diagonal scaling, preconditioners

 79:   Notes: By using KSPSetPreconditionerSide(ksp,PC_SYMMETRIC) or -ksp_symmetric_pc you 
 80:          can scale each side of the matrix by the squareroot of the diagonal entries.

 82:          Zero entries along the diagonal are replaced with the value 1.0

 84: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC
 85: M*/


 91: PetscErrorCode PCCreate_Prometheus(PC pc)
 92: {
 94: 

 97:   PCCreate_Prometheus_private( pc );
 98: 
 99:   /*
100:     Set the pointers for the functions that are provided above.
101:     Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
102:     are called, they will automatically call these functions.  Note we
103:     choose not to provide a couple of these functions since they are
104:     not needed.
105:   */
106:   pc->ops->apply               = PCApply_Prometheus;
107:   pc->ops->applytranspose      = PCApply_Prometheus;
108:   pc->ops->setup               = PCSetUp_Prometheus;
109:   pc->ops->destroy             = PCDestroy_Prometheus;
110:   pc->ops->setfromoptions      = PCSetFromOptions_Prometheus;
111:   pc->ops->view                = PCView_Prometheus;
112:   pc->ops->applyrichardson     = 0;
113:   pc->ops->applysymmetricleft  = 0;/* PCApplySymmetricLeftOrRight_Prometheus; */
114:   pc->ops->applysymmetricright = 0;/* PCApplySymmetricLeftOrRight_Prometheus; */
115:   PetscObjectComposeFunctionDynamic( (PetscObject)pc,
116:                                             "PCSetCoordinates_C",
117:                                             "PCSetCoordinates_Prometheus",
118:                                             PCSetCoordinates_Prometheus);
119:   PetscObjectComposeFunctionDynamic( (PetscObject)pc,
120:                                             "PCSASetVectors_C",
121:                                             "PCSASetVectors_Prometheus",
122:                                             PCSASetVectors_Prometheus);
123:   return(0);
124: }

129: /*@
130:    PCSetCoordinates - sets the coordinates of all the nodes on the local process

132:    Collective on PC

134:    Input Parameters:
135: +  pc - the solver context
136: .  dim - the dimension of the coordinates 1, 2, or 3
137: -  coords - the coordinates

139:    Level: intermediate

141:    Notes: coords is an array of the 3D coordinates for the nodes on
142:    the local processor.  So if there are 108 equation on a processor
143:    for a displacement finite element discretization of elasticity (so
144:    that there are 36 = 108/3 nodes) then the array must have 108
145:    double precision values (ie, 3 * 36).  These x y z coordinates
146:    should be ordered for nodes 0 to N-1 like so: [ 0.x, 0.y, 0.z, 1.x,
147:    ... , N-1.z ].

149: .seealso: PCPROMETHEUS
150: @*/
151: PetscErrorCode  PCSetCoordinates(PC pc,PetscInt dim,PetscReal *coords)
152: {

156:   PetscTryMethod(pc,PCSetCoordinates_C,(PC,PetscInt,PetscReal*),(pc,dim,coords));
157:   return(0);
158: }

162: /*@
163:    PCSASetVectors - sets the vectors of all the nodes on the local process

165:    Collective on PC

167:    Input Parameters:
168: +  pc - the solver context
169: .  nects - the number of vectors
170: -  vects - the vectors

172:    Level: intermediate

174:    Notes: 'vects' is a dense tall skinny matrix with 'nvects' columns and 
175:    the number of local equations rows.  'vects' is stored in row major order.

177: .seealso: PCPROMETHEUS
178: @*/
179: PetscErrorCode  PCSASetVectors(PC pc,PetscInt nvects,PetscReal *vects)
180: {

184:   PetscTryMethod(pc,PCSASetVectors_C,(PC,PetscInt,PetscReal*),(pc,nvects,vects));
185:   return(0);
186: }