System Preprocessors
preprocess.c
Go to the documentation of this file.
00001 /*! \mainpage SysPro: a System preProcessors library
00002 
00003 
00004   \section Introduction
00005 
00006   This is a library of preprocessors for numerical problems, that is,
00007   mappings of one numerical problem into another, presumably more simple one,
00008   of the same type. For example, scaling a linear system.
00009   The SysPro library operates in two modes:
00010   - exhaustive mode: all possibly choices of a preprocessor are explored
00011   in sequence; this mode is set by command line options.
00012   - intelligent mode: based on problem properties, a suitable choice for
00013   each preprocessor is made.
00014 
00015   See \ref modes for more details.
00016 
00017   Each preprocessor has the following structure, which is executed
00018   in the PreprocessedSolution() routine:
00019   - a global setup is performed. This is a good place for computing
00020   problem features with the AnaMod library.
00021   - a specific setup is performed; this can for instance disable certain
00022   preprocessor choices based on the computed problem features.
00023   - a selection is made; this can be
00024     - a first choice, if the preprocessor is applied in exhaustive mode
00025     - an intelligent choice, if the exhaustive mode is off, and an 
00026     intelligent choice routine exists
00027     - some default choice otherwise
00028 
00029   In case of exhaustive mode, the following steps are executed inside a loop
00030   over all choices for this preprocessor, and possible all numerical option
00031   settings:
00032   - the start function transforms the problem into a preprocessed problem
00033   - if a next preprocessor is defined, it is applied;
00034   otherwise, the problem solving routine is applied 
00035   (see \ref declaration).
00036   - the end function backtransforms the solution of the preprocessed problem
00037   into that of the original problem.
00038 
00039   \section Topics
00040 
00041   \subpage use
00042 
00043   \subpage reporting and \subpage tracing
00044 
00045   \subpage linear
00046 
00047   \subpage external
00048 
00049   \subpage options
00050 
00051   \author Victor Eijkhout
00052   \version 1.3
00053   \date unreleased
00054 
00055   \section diary Change log
00056 
00057   1.3
00058   2008/08/20 : DeclarePreprocessor() has an extra argument for global
00059   deallocation at the end of a program run. Currently used in the ksp
00060   preprocessor.
00061 
00062   2008/05/10 : end function now has two NumericalProblem parameters;
00063   this is necessary for freeing the recursive problem.
00064 
00065   2007 : Accomodated the array length parameter in anamod/nmd
00066 
00067 */
00068 
00069 /*! \page use Use of the SysPro package
00070 
00071   \section sysprosetup Setup of the SysPro system
00072 
00073   The setup of SysPro has to be done once per program run.
00074 
00075   \subsection sysprogsetup Global declarations
00076 
00077   All use of SysPro has to be inside calls to SysProInitialize() and
00078   SysProFinalize(). These allocate and deallocate global data structures.
00079   You could place them right next to MPI_Initialize/Finalize or 
00080   PetscInitialize/Finalize calls. The SysProFinalize() call can also be used
00081   to deallocate data that was constructed during preprocessor setup.
00082 
00083   Another step in the global setup of SysPro is a call to 
00084   SysProDeclareFunctions(), This declares functions that are of use
00085   to all preprocessors that will be declared later.
00086 
00087   \subsection declaration Preprocessor declaration
00088 
00089   Preprocessors are declared with calls to DeclarePreprocessor(), which 
00090   installs the setup functions and the forward / backward transformations.
00091 
00092   Further specifications can be given for a specific preprocessor:
00093   - DeclarePreprocessorIntelligentChoice() for installing a routine that will
00094   intelligently pick a preprocessor choice.
00095   - see \ref anamodinterface for metadata category handling.
00096 
00097   \section setup Setup and invocation
00098 
00099   After the setup as described above,
00100   PreprocessorsOptionsHandling() can be called to provide the user with 
00101   runtime control (see section \ref modes) over the workings of SysPro.
00102 
00103   Preprocessed problem solving is activated by a call to
00104   PreprocessedProblemSolving(). This causes all declared preprocessors
00105   to be applied in sequence. Finally, the ultimately remaining problem is
00106   solved with the routine declared by PreprocessorsDeclareProblemSolver().
00107 
00108   \section modes Usage modes
00109 
00110   SysPro can be used in several modes:
00111   - one can specify the exact preprocessor values (or several values);
00112   - one can specify to test exhaustively the values of one preprocessor
00113   or all of them;
00114   - SysPro can intelligently pick the appropriate preprocessor.
00115 
00116   The intelligent preprocessor choice uses a model where each preprocessor
00117   has a measure of how applicable it is; the SysPro system then picks the 
00118   most appropriate preprocessor from a given class. See \ref suit.
00119 
00120   See \ref options for details on specific and exhaustive testing.
00121 
00122   \section modules Predefined preprocessors
00123 
00124   The SysPro package comes with a number of predefined preprocessors,
00125   mostly geared towards linear system processing.
00126 
00127   \page linear
00128 
00129   \page lineartransform
00130  */
00131 
00132 /*! \page external The interface to other packages
00133 
00134   SysPro by itself is not of a lot of use: it is a framework for 
00135   tying together transformations and operations. There is also an interface
00136   for computing metadata.
00137 
00138   \subpage anamodinterface
00139 
00140  */
00141 
00142 /*! \page linear The linear systems example package
00143   The mechanism of preprocessed solving through forward and backward
00144   transformations can be applied to all sorts of numerical problems.
00145   At present, we only supply code for linear system solving; 
00146   see \ref linearfile.
00147  */
00148 
00149 #include <stdlib.h>
00150 #include <stdio.h>
00151 #include <string.h>
00152 #include "petscmat.h"
00153 #include "syspro.h"
00154 #include "sysprotransform.h"
00155 #include "syspro_impl.h"
00156 
00157 #define NPREPROCESS 25
00158 SystemPreprocessor *preprocessors=NULL;
00159 const char **currentpreprocessors,**currentchoices=NULL;
00160 int *currentoptions=NULL;
00161 int npreprocess=0,preprocesslevel;
00162 static void **preprocessorcontexts,*solutioncontext;
00163 static PetscErrorCode(**unsetpreprocessor)();
00164 
00165 struct PreprocessorsGlobalInfo_ {
00166   PetscErrorCode(*problemmonitor)(NumericalProblem);
00167   PetscErrorCode(*classstaticsetup)(const char*);
00168   /** This routine is executed on the creation of a new
00169       preprocessor. It can be used to install standard options
00170       in the preprocessor transform object.
00171   */
00172   PetscErrorCode(*classdynamicsetup)(const char*,NumericalProblem);
00173   /** This routine is invoked at the start of each preprocessor
00174       class. It is not supposed to contain problem-dependent
00175       actions. It is useful for printing trace messages, and
00176       performing analysis on each incoming problem.
00177   */      
00178   PetscErrorCode(*classproblemcloner)
00179   (const char*,const char*,int,NumericalProblem,NumericalProblem);
00180   /** This routine is called everytime a new problem is created with
00181       a class/option pair. It can be used to copy preserved metadata elements
00182   */
00183   PetscErrorCode(*computecategory)(char*,NumericalProblem);
00184   /** This routine is called in sequence with the names of the
00185       required metadata categories.
00186    */
00187   PetscErrorCode(*metadatacomputer)(char*,char*,Mat,void*,PetscTruth*);
00188   PetscErrorCode(*clonecontext)(const char*,const char*,void*,void**);
00189   PetscErrorCode(*freecontext)(void*);
00190   PetscErrorCode(*problemsolver)(NumericalProblem,void*,NumericalSolution*);
00191   PetscErrorCode(*problemdelete)(NumericalProblem);
00192   PetscErrorCode(*errortracer)(NumericalProblem,NumericalSolution,const char*);
00193   PetscErrorCode(*solutioncreator)(NumericalProblem,NumericalSolution*);
00194   PetscErrorCode(*solutioncopy)(NumericalSolution,NumericalSolution);
00195   PetscErrorCode(*solutiondelete)(NumericalSolution);
00196   PetscErrorCode(*solutioncontextdelete)(NumericalSolution);
00197 } ; 
00198 typedef struct PreprocessorsGlobalInfo_* PreprocessorsGlobalInfo;
00199 static PreprocessorsGlobalInfo GlobalInfo = NULL;
00200 
00201 #undef __FUNCT__
00202 #define __FUNCT__ "CreateGlobalInfo"
00203 static PetscErrorCode CreateGlobalInfo()
00204 {
00205   PetscErrorCode ierr;
00206   PetscFunctionBegin;
00207   if (!GlobalInfo) {
00208     ierr = PetscMalloc
00209       (sizeof(struct PreprocessorsGlobalInfo_),&GlobalInfo); CHKERRQ(ierr);
00210     ierr = PetscMemzero
00211       (GlobalInfo,sizeof(struct PreprocessorsGlobalInfo_)); CHKERRQ(ierr);
00212   }
00213   CHKMEMQ;
00214   PetscFunctionReturn(0);
00215 }
00216 
00217 #undef __FUNCT__
00218 #define __FUNCT__ "SysProInitialize"
00219 /*! Allocate SysPro globals. See also SysProFinalize(). */
00220 PetscErrorCode SysProInitialize()
00221 {
00222   PetscErrorCode ierr;
00223   PetscFunctionBegin;
00224   ierr = CreateGlobalInfo(); CHKERRQ(ierr);
00225   ierr = PetscMalloc
00226     (NPREPROCESS*sizeof(SystemPreprocessor),&preprocessors); CHKERRQ(ierr);
00227   ierr = PetscMalloc
00228     (NPREPROCESS*sizeof(char*),&currentpreprocessors); CHKERRQ(ierr);
00229   ierr = PetscMalloc
00230     (NPREPROCESS*sizeof(void*),&preprocessorcontexts); CHKERRQ(ierr);
00231   ierr = PetscMalloc
00232     (NPREPROCESS*sizeof(char*),&currentchoices); CHKERRQ(ierr);
00233   ierr = PetscMalloc
00234     (NPREPROCESS*sizeof(int),&currentoptions); CHKERRQ(ierr);
00235   ierr = PetscMalloc
00236     (NPREPROCESS*sizeof(PetscErrorCode(*)()),&unsetpreprocessor); CHKERRQ(ierr);
00237   ierr = PetscMemzero
00238     (unsetpreprocessor,NPREPROCESS*sizeof(PetscErrorCode(*)())); CHKERRQ(ierr);
00239   CHKMEMQ;
00240   PetscFunctionReturn(0);
00241 }
00242 
00243 #undef __FUNCT__
00244 #define __FUNCT__ "SysProFinalize"
00245 PetscErrorCode SysProFinalize()
00246 {
00247   int ip; PetscErrorCode ierr;
00248   PetscFunctionBegin;
00249   for (ip=0; ip<npreprocess; ip++) {
00250     SystemPreprocessor thisone = preprocessors[ip];
00251     if (unsetpreprocessor[ip]) {
00252       ierr = (unsetpreprocessor[ip])(); CHKERRQ(ierr);
00253     }
00254     if (thisone->preserved) {
00255       ierr = PetscFree(thisone->preserved); CHKERRQ(ierr);
00256     }
00257     ierr = PetscFree(thisone->name); CHKERRQ(ierr);
00258     ierr = DeregisterTransform(thisone->transform); CHKERRQ(ierr);
00259     ierr = PetscFree(thisone); CHKERRQ(ierr);
00260   }
00261   ierr = PetscFree(preprocessors); CHKERRQ(ierr);
00262   ierr = PetscFree(currentpreprocessors); CHKERRQ(ierr);
00263   ierr = PetscFree(preprocessorcontexts); CHKERRQ(ierr);
00264   ierr = PetscFree(currentchoices); CHKERRQ(ierr);
00265   ierr = PetscFree(currentoptions); CHKERRQ(ierr);
00266   ierr = PetscFree(unsetpreprocessor); CHKERRQ(ierr);
00267   ierr = PetscFree(GlobalInfo); CHKERRQ(ierr);
00268   CHKMEMQ;
00269   PetscFunctionReturn(0);
00270 }
00271 
00272 #undef __FUNCT__
00273 #define __FUNCT__ "DeclarePreprocessor"
00274 /*!
00275   Declare a preprocessor class, by specifing its various members.
00276 
00277   The \c name argument should not contain the colon character.
00278 
00279   Here is an explanation of the various function arguments.
00280 
00281   \c this_preprocessor_setup() : this routine is called only once,
00282   inside this function. This is a good place for defining all the 
00283   preprocessors in this class
00284 
00285   \c specific_setup(NumericalProblem,SalsaTransform) : this is called
00286   at the start of a preprocessing stage; one could use this for
00287   computing matrix metadata.
00288 
00289   \c global_unset(void) : this is called in SysProFinalize().
00290 
00291   \c ctx_create(NumericalProblem,void**) : create an object that can
00292   be used for the duration of the application of this preprocessor
00293 
00294   \c ctxdelete(void*) : delete the context again
00295 
00296   \c start_function : this is the function that performs the forward
00297   transform of the problem. Prototype:
00298 \code
00299   PetscErrorCode start_function
00300     (char             *classmember,
00301      int               optionvalue,
00302      PetscTruth        overwrite,
00303      NumericalProblem  problem,
00304      NumericalProblem *transformedproblem,
00305      void             *globalcontext,
00306      void            **localcontext,
00307      PetscTruth       *success)
00308 \endcode
00309 
00310   \c end_function : this is the backtransform. Its main task is copying
00311   or backtransforming the preprocessed solution to the original solution.
00312 \code
00313   PetscErrorCode end_function
00314     (char             *classmember,
00315      PetscTruth        overwrite,
00316      void             *globalcontext,
00317      void             *localcontext,
00318      NumericalProblem  pproblem,
00319      NumericalProblem  oproblem,
00320      NumericalSolution psolution,
00321      NumericalSolution osolution)
00322 \endcode
00323 where \c pproblem and \c psolution are the preprocessed quantities, 
00324 the end function has to unprocess them and leave the result in
00325 \c oproblem, \c osolution. Actually, \c oproblem is only for reference.
00326 */
00327 PetscErrorCode DeclarePreprocessor
00328 (const char *name,
00329  PetscErrorCode(*this_preprocessor_setup)(),
00330  PetscErrorCode(*specific_setup)(NumericalProblem,SalsaTransform),
00331  PetscErrorCode(*specific_unset)(NumericalProblem),
00332  PetscErrorCode(*global_unset)(),
00333  PetscErrorCode(*ctxcreate)(NumericalProblem,void**),
00334  PetscErrorCode(*ctxdelete)(void*),
00335  PetscErrorCode(*start_function)(const char*,int,PetscTruth,NumericalProblem,NumericalProblem*,void*,void**,PetscTruth*),
00336  PetscErrorCode(*end_function)(const char*,PetscTruth,void*,void*,NumericalProblem,NumericalProblem,NumericalSolution,NumericalSolution)
00337  )
00338 {
00339   SystemPreprocessor snew;
00340   PetscErrorCode ierr;
00341 
00342   PetscFunctionBegin;
00343   if (!preprocessors) SETERRQ(1,"Forgotten SysProInitialize");
00344   if (npreprocess>=NPREPROCESS)
00345     SETERRQ(1,"Can not declare further system preprocessors");
00346   {
00347     int p;
00348     for (p=0; p<npreprocess; p++) 
00349       if (!strcmp(name,preprocessors[p]->name))
00350         SETERRQ1(1,"There already is a preprocessor called <%s>",name);
00351   }
00352   unsetpreprocessor[npreprocess] = global_unset;
00353   ierr = PetscMalloc(sizeof(struct SystemPreprocessor_),&snew); CHKERRQ(ierr);
00354   ierr = PetscMemzero(snew,sizeof(struct SystemPreprocessor_)); CHKERRQ(ierr);
00355   preprocessors[npreprocess++] = snew;
00356 
00357   snew->exhaustive = PETSC_FALSE;
00358   {
00359     char *name;
00360     ierr = PetscStrallocpy(name,&name); CHKERRQ(ierr);
00361     snew->name = name;
00362   }
00363   ierr = NewTransform(name,&(snew->transform)); CHKERRQ(ierr);
00364   if (GlobalInfo->classstaticsetup) {
00365     ierr = (GlobalInfo->classstaticsetup)(name); CHKERRQ(ierr);
00366   }
00367   if (!this_preprocessor_setup)
00368     SETERRQ(1,"No setup function given");
00369   ierr = (*this_preprocessor_setup)(); CHKERRQ(ierr);
00370   snew->setup = specific_setup;
00371   snew->unset = specific_unset;
00372   snew->ctxcreate = ctxcreate; snew->ctxdelete = ctxdelete;
00373   snew->start_function = start_function; snew->end_function = end_function;
00374   CHKMEMQ;
00375   PetscFunctionReturn(0);
00376 }
00377 
00378 #undef __FUNCT__
00379 #define __FUNCT__ "SysProDeclareProblemMonitor"
00380 PetscErrorCode SysProDeclareProblemMonitor
00381 (PetscErrorCode(*monitor)(NumericalProblem))
00382 {
00383   PetscFunctionBegin;
00384   if (monitor) GlobalInfo->problemmonitor = monitor;
00385   PetscFunctionReturn(0);
00386 }
00387 
00388 #undef __FUNCT__
00389 #define __FUNCT__ "SysProDeclareErrorTracer"
00390 PetscErrorCode SysProDeclareErrorTracer
00391 (PetscErrorCode(*tracer)(NumericalProblem,NumericalSolution,const char*))
00392 {
00393   PetscFunctionBegin;
00394   if (tracer) GlobalInfo->errortracer = tracer;
00395   PetscFunctionReturn(0);
00396 }
00397 
00398 #undef __FUNCT__
00399 #define __FUNCT__ "SysProGetErrorTracer"
00400 PetscErrorCode SysProGetErrorTracer
00401 (PetscErrorCode(**tracer)(NumericalProblem,NumericalSolution,const char*))
00402 {
00403   PetscFunctionBegin;
00404   if (tracer) *tracer = GlobalInfo->errortracer;
00405   PetscFunctionReturn(0);
00406 }
00407 
00408 #undef __FUNCT__
00409 #define __FUNCT__ "DeclarePreprocessorIntelligentChoice"
00410 /*! Install a function to pick the optimal choice for a preprocessor */
00411 PetscErrorCode DeclarePreprocessorIntelligentChoice
00412 (const char *name,PetscErrorCode(*picker)(NumericalProblem,const char**,const char**))
00413 {
00414   SystemPreprocessor pp;
00415   PetscErrorCode ierr;
00416   PetscFunctionBegin;
00417   ierr = SystemPreprocessorGetByName(name,&pp); CHKERRQ(ierr);
00418   pp->intelligence = picker;
00419   PetscFunctionReturn(0);
00420 }
00421 
00422 #undef __FUNCT__
00423 #define __FUNCT__ "NumericalProblemGetComm"
00424 PetscErrorCode NumericalProblemGetComm(NumericalProblem prob,MPI_Comm *comm)
00425 {
00426   PetscFunctionBegin;
00427   *comm = prob->comm;
00428   PetscFunctionReturn(0);
00429 }
00430 
00431 #undef __FUNCT__
00432 #define __FUNCT__ "SysProDeclareFunctions"
00433 /*! Install various functions
00434   - \c classstaticsetup : this function is called on each processor as
00435   it is being created; see DeclarePreprocessor().
00436   - \c classdynamicsetup : this function is called as any invocation
00437   of a preprocessor starts; see PreprocessedSolution();
00438   - \c classproblemcloner : a function to clone the context : optional
00439   see \ref tracing for more details.
00440   - \c problemsolver : the ultimate problem solver : required
00441   - \c problemdelete : delete a problem object
00442   - \c solutioncreator : creates a solution object; optional, but required
00443   a preprocessor has an endfunction.
00444   - \c solutioncopy : guess what this does; optional
00445   - \c solutiondelete : optional, but needed if solutioncopy is used
00446   - \c contextcloner : problems can carry a context; this clones the context
00447   if a problem is copied; otherwise the pointer is simply duplicated
00448   - \c contextfree : used to delete cloned contexts
00449   - \c solutioncontextdelete : hm.
00450 
00451 */
00452 PetscErrorCode SysProDeclareFunctions
00453 (
00454  PetscErrorCode(*classstaticsetup)(const char*),
00455  PetscErrorCode(*classdynamicsetup)(const char*,NumericalProblem),
00456  PetscErrorCode(*classproblemcloner)(const char*,const char*,int,NumericalProblem,NumericalProblem),
00457  PetscErrorCode(*problemsolver)(NumericalProblem,void*,NumericalSolution*),
00458  PetscErrorCode(*problemdelete)(NumericalProblem),
00459  PetscErrorCode(*solutioncreator)(NumericalProblem,NumericalSolution*),
00460  PetscErrorCode(*solutioncopy)(NumericalSolution,NumericalSolution),
00461  PetscErrorCode(*solutiondelete)(NumericalSolution),
00462  PetscErrorCode(*ctxcloner)(const char*,const char*,void*,void**),
00463  PetscErrorCode(*ctxfree)(void*),
00464  PetscErrorCode(*solutioncontextdelete)(NumericalSolution)
00465  )
00466 {
00467   PetscFunctionBegin;
00468 
00469   if (classstaticsetup) GlobalInfo->classstaticsetup = classstaticsetup;
00470   if (classdynamicsetup) GlobalInfo->classdynamicsetup = classdynamicsetup;
00471   if (classproblemcloner) GlobalInfo->classproblemcloner = classproblemcloner;
00472 
00473   if (!problemsolver) SETERRQ(1,"A solver routine is required");
00474   GlobalInfo->problemsolver = problemsolver;
00475   GlobalInfo->problemdelete = problemdelete;
00476 
00477   if (ctxcloner) GlobalInfo->clonecontext = ctxcloner;
00478   if (ctxfree) GlobalInfo->freecontext = ctxfree;
00479 
00480   if (solutioncreator) GlobalInfo->solutioncreator = solutioncreator;
00481   if (solutioncopy) GlobalInfo->solutioncopy = solutioncopy;
00482   if (solutiondelete) GlobalInfo->solutiondelete = solutiondelete;
00483   if (solutioncontextdelete)
00484     GlobalInfo->solutioncontextdelete = solutioncontextdelete;
00485   CHKMEMQ;
00486 
00487   PetscFunctionReturn(0);
00488 }
00489 
00490 /*! \page context Context handling
00491   
00492   In order to carry application-specific information and temporaries
00493   around, there are a few opaque handle contexts in syspro.
00494 
00495   A NumericalProblem structure is defined to have a context. This is
00496   cloned at the application of a preprocessor, and deleted when its
00497   application is finished.
00498 
00499   See SysProGetContextFunctions(), SysProProblemCloneContext(),
00500   SysProProblemDeleteContext().
00501  */
00502 #undef __FUNCT__
00503 #define __FUNCT__ "SysProGetContextFunctions"
00504 static PetscErrorCode SysProGetContextFunctions
00505 (PetscErrorCode(**ctxcloner)(const char*,const char*,void*,void**),
00506  PetscErrorCode(**ctxfree)(void*)
00507  )
00508 {
00509   PetscFunctionBegin;
00510   if (ctxcloner) *ctxcloner = GlobalInfo->clonecontext;
00511   if (ctxfree) *ctxfree = GlobalInfo->freecontext;
00512   CHKMEMQ;
00513   PetscFunctionReturn(0);
00514 }
00515 
00516 static PetscErrorCode SysProProblemCloneContext
00517 (const char *preprocessor,const char *type,NumericalProblem in,NumericalProblem out)
00518 {
00519   PetscErrorCode(*cloner)(const char*,const char*,void*,void**),ierr;
00520   PetscFunctionBegin;
00521   ierr = SysProGetContextFunctions(&cloner,PETSC_NULL); CHKERRQ(ierr);
00522   if (cloner) {
00523     ierr = (*cloner)(preprocessor,type,in->ctx,&(out->ctx)); CHKERRQ(ierr);
00524   }
00525   CHKMEMQ;
00526   PetscFunctionReturn(0);
00527 }
00528 
00529 static PetscErrorCode SysProProblemDeleteContext(NumericalProblem problem)
00530 {
00531   PetscErrorCode(*freer)(void*),ierr;
00532   PetscFunctionBegin;
00533   ierr = SysProGetContextFunctions(PETSC_NULL,&freer); CHKERRQ(ierr);
00534   /*
00535   if (freer && problem->ctx) {
00536     ierr = (*freer)(problem->ctx); CHKERRQ(ierr);
00537   }
00538   */
00539   CHKMEMQ;
00540   PetscFunctionReturn(0);
00541 }
00542 
00543 #undef __FUNCT__
00544 #define __FUNCT__ "RegisterPreprocessorSetting"
00545 static PetscErrorCode RegisterPreprocessorSetting
00546 (const char *preprocess,const char *type,int option)
00547 {
00548   PetscFunctionBegin;
00549   currentpreprocessors[preprocesslevel] = preprocess;
00550   currentchoices[preprocesslevel] = type;
00551   currentoptions[preprocesslevel] = option;
00552   CHKMEMQ;
00553   PetscFunctionReturn(0);
00554 }
00555 
00556 #undef __FUNCT__
00557 #define __FUNCT__ "PreprocessorGetSetting"
00558 PetscErrorCode PreprocessorGetSetting
00559 (const char *preprocess,const char **type,int *option)
00560 {
00561   int i;
00562   PetscFunctionBegin;
00563   for (i=0; i<=preprocesslevel; i++) {
00564     if (!strcmp(currentpreprocessors[i],preprocess)) {
00565       *type = currentchoices[i];
00566       *option = currentoptions[i];
00567       goto found;
00568     }
00569   }
00570   SETERRQ1(1,"Prprocessor <%s> is not active",preprocess);
00571  found:
00572   CHKMEMQ;
00573   PetscFunctionReturn(0);
00574 }
00575 
00576 
00577 #undef __FUNCT__
00578 #define __FUNCT__ "RetrievePreprocessorChoice"
00579 PetscErrorCode RetrievePreprocessorChoice(int idx,const char **type,int *option)
00580 {
00581   PetscFunctionBegin;
00582   *type = currentchoices[idx];
00583   *option = currentoptions[idx];
00584   CHKMEMQ;
00585   PetscFunctionReturn(0);
00586 }
00587 
00588 
00589 #undef __FUNCT__
00590 #define __FUNCT__ "PreprocessorGetIndex"
00591 PetscErrorCode PreprocessorGetIndex(const char *name,int *prenumber)
00592 {
00593   int ipre; PetscTruth flg; PetscErrorCode ierr;
00594   PetscFunctionBegin;
00595   *prenumber = -1;
00596   for (ipre=0; ipre<npreprocess; ipre++) {
00597     ierr = PetscStrcmp(name,preprocessors[ipre]->name,&flg); CHKERRQ(ierr);
00598     if (flg) {
00599       *prenumber = ipre;
00600       goto found;
00601     }
00602   }
00603   SETERRQ1(1,"Could not find preprocessor <%s>",name);
00604  found:
00605   CHKMEMQ;
00606   PetscFunctionReturn(0);
00607 }
00608 
00609 #undef __FUNCT__
00610 #define __FUNCT__ "SystemPreprocessorGetByName"
00611 PetscErrorCode SystemPreprocessorGetByName
00612 (const char *name,SystemPreprocessor *pp)
00613 {
00614   int ipre; PetscErrorCode ierr;
00615   PetscFunctionBegin;
00616   ierr = PreprocessorGetIndex(name,&ipre); CHKERRQ(ierr);
00617   *pp = preprocessors[ipre];
00618   CHKMEMQ;
00619   PetscFunctionReturn(0);
00620 }
00621 
00622 #undef __FUNCT__
00623 #define __FUNCT__ "TransformGetByName"
00624 PetscErrorCode TransformGetByName(const char *name,SalsaTransform *tf)
00625 {
00626   SystemPreprocessor pp; PetscErrorCode ierr;
00627   PetscFunctionBegin;
00628   ierr = SystemPreprocessorGetByName(name,&pp); CHKERRQ(ierr);
00629   *tf = pp->transform;
00630   CHKMEMQ;
00631   PetscFunctionReturn(0);
00632 }
00633 
00634 /*! \page context Preprocessor contexts
00635   In order to offer flexibility, there are several possibilities
00636   for user objects to be stored under opaque (\c void*) pointers.
00637 
00638   \section Global contexts
00639 
00640   Each preprocessor that has declared a \c contextcreate
00641   function, will create that context at the start of its traversal.
00642   This context is then globally registered with RegisterPreprocessorContext()
00643   under the name of this preprocessor, so that other preprocessor
00644   can have access to it with PreprocessorGetContext().
00645 
00646   There is one special context, which is stored under the \c solution
00647   handle. This one is not created by default. See the "linear" package
00648   for an example of how to use it.
00649 
00650   \section Local contexts
00651   
00652   The \c start_function can create a context, which is input
00653   to the \c end_function. This context serves to preser data that is
00654   necessary for the inverse transformation that is applied in the
00655   \c start_function
00656 */
00657 #undef __FUNCT__
00658 #define __FUNCT__ "RegisterPreprocessorContext"
00659 PetscErrorCode RegisterPreprocessorContext(const char *pre,void *ctx)
00660 {
00661   PetscErrorCode ierr;
00662   PetscFunctionBegin;
00663   if (!preprocessorcontexts) {
00664     SETERRQ(1,"Preprocessor contexts array should have been allocated");
00665   } else {
00666     PetscTruth flg; int ipre;
00667     ierr = PetscStrcmp(pre,"solution",&flg); CHKERRQ(ierr);
00668     if (flg) {
00669       solutioncontext = ctx;
00670     } else {
00671       ierr = PreprocessorGetIndex(pre,&ipre); CHKERRQ(ierr);
00672       preprocessorcontexts[ipre] = ctx;
00673     }
00674   }
00675   CHKMEMQ;
00676   PetscFunctionReturn(0);
00677 }
00678 
00679 #undef __FUNCT__
00680 #define __FUNCT__ "PreprocessorGetContext"
00681 PetscErrorCode PreprocessorGetContext(const char *pre,void **ctx)
00682 {
00683   PetscTruth flg; int ipre; PetscErrorCode ierr;
00684   PetscFunctionBegin;
00685   ierr = PetscStrcmp(pre,"solution",&flg); CHKERRQ(ierr);
00686   if (flg) {
00687     *ctx = solutioncontext;
00688   } else {
00689     ierr = PreprocessorGetIndex(pre,&ipre); CHKERRQ(ierr);
00690     if (!preprocessorcontexts[ipre])
00691       SETERRQ1(1,"No context set for preprocessor <%s>",pre);
00692     *ctx = preprocessorcontexts[ipre];
00693   }
00694   CHKMEMQ;
00695   PetscFunctionReturn(0);
00696 }
00697 
00698 /*! Setup actions that are particular to a specific class of preprocessors,
00699   such as scalings. This performs the following actions:
00700   - any user specified setup, for instance disqualifying preprocessors
00701   on purely logistical grounds.
00702   - subsequently, the suitability functions (\ref suit) are evaluated for each
00703   preprocessor and the current problem. In the current implementation
00704   this is only used to disqualify preprocessors. We'll get more
00705   sophisticated later.
00706 */
00707 #undef __FUNCT__
00708 #define __FUNCT__ "PreprocessorSpecificSetup"
00709 static PetscErrorCode PreprocessorSpecificSetup
00710 (const char *preprocess,NumericalProblem problem,PetscTruth user_choices)
00711 {
00712   SystemPreprocessor pp; SalsaTransform tf; PetscErrorCode ierr;
00713   PetscFunctionBegin;
00714   ierr = SystemPreprocessorGetByName(preprocess,&pp); CHKERRQ(ierr);
00715   ierr = TransformGetByName(preprocess,&tf); CHKERRQ(ierr);
00716 
00717   if (GlobalInfo->problemmonitor) {
00718     ierr = (GlobalInfo->problemmonitor)(problem); CHKERRQ(ierr);
00719   }
00720 
00721   /* if the user specified what choices to take, all others have
00722      already been disabled
00723   */
00724   if (!user_choices) {
00725     int i;
00726     /* declare all preprocessor choices as suitable, then 
00727        eliminate some if there is a setup function for it
00728     */
00729     ierr = PreprocessorApplyAprioriSelection(pp); CHKERRQ(ierr);
00730     /*ierr = TransformObjectsUnmarkAll(tf); CHKERRQ(ierr);*/
00731     if (pp->setup) {
00732       ierr = (pp->setup)(problem,tf); CHKERRQ(ierr);
00733     }
00734     /* Now evaluate any suitability functions, and mark preprocessors 
00735        accordingly.
00736     */
00737     for (i=0; i<tf->n_objects; i++) {
00738       SalsaTransformObject to = tf->transformobjects[i]; void *sctx;
00739       PetscErrorCode(*f)(NumericalProblem,void*,SuitabilityValue*);
00740       ierr = TransformObjectGetSuitabilityFunction(to,&sctx,&f); CHKERRQ(ierr);
00741       if (f) {
00742         SuitabilityValue v; const char *name;
00743         ierr = TestSuitabilityContext(tf,sctx); CHKERRQ(ierr);
00744         ierr = (*f)(problem,sctx,&v); CHKERRQ(ierr);
00745         ierr = TransformObjectGetName(to,&name); CHKERRQ(ierr);
00746         printf("method %s compute suitability as %e\n",name,v);
00747         if (v==0) {
00748           ierr = TransformObjectMark(to); CHKERRQ(ierr);
00749         }
00750       }
00751     }
00752   }
00753   ierr = ReportEnabledPreprocessors(preprocess); CHKERRQ(ierr);
00754   CHKMEMQ;
00755   PetscFunctionReturn(0);
00756 }
00757 
00758 #undef __FUNCT__
00759 #define __FUNCT__ "SysproPreprocessorStartFunction"
00760 static PetscErrorCode SysproPreprocessorStartFunction
00761 (SystemPreprocessor preprocessor,const char *type,int option,PetscTruth overwrite,
00762  NumericalProblem inproblem,NumericalProblem *outproblem,
00763  void *preprocessor_context,void **transform_context, PetscTruth *success)
00764 {
00765   PetscErrorCode ierr;
00766   PetscFunctionBegin;
00767   if (!(preprocessor->start_function)) 
00768     SETERRQ1(1,"Preprocessor <%s> has no start function",preprocessor->name);
00769   ierr = (preprocessor->start_function)
00770     (type,option,overwrite,inproblem,outproblem,
00771      preprocessor_context,transform_context,success); CHKERRQ(ierr);
00772   if (!*success) PetscFunctionReturn(0);
00773   /* what on earth ?!
00774   if (GlobalInfo->classproblemcloner) {
00775     ierr = (GlobalInfo->classproblemcloner)
00776       (preprocessor->name,type,option,inproblem,*outproblem); CHKERRQ(ierr);
00777   } else *outproblem = inproblem;
00778   */
00779   CHKMEMQ;
00780   PetscFunctionReturn(0);
00781 }
00782 
00783 #undef __FUNCT__
00784 #define __FUNCT__ "SysProPreprocessorEndFunction"
00785 static PetscErrorCode SysProPreprocessorEndFunction
00786 (SystemPreprocessor preprocessor,const char *pclassname,const char *type, PetscTruth do_not_keep_original,void *gctx,void *ctx,
00787  NumericalProblem next_problem,NumericalProblem problem, NumericalSolution next_solution,NumericalSolution *rsolution)
00788 {
00789   NumericalSolution solution; PetscErrorCode ierr;
00790   PetscFunctionBegin;
00791 
00792   /* if there are separate solutions on this level and the surrounding,
00793    * the outer one is created here and the nested one is deleted.
00794    */
00795   if (next_solution && (preprocessor->end_function || GlobalInfo->solutioncopy)) {
00796     if (!GlobalInfo->solutioncreator) SETERRQ1(1,"Preprocessor %s requires solution creator",pclassname);
00797     if (!GlobalInfo->solutiondelete)  SETERRQ1(1,"Preprocessor %s requires solution destructor",pclassname);
00798     ierr = (GlobalInfo->solutioncreator)(problem,&solution); CHKERRQ(ierr);
00799     if (preprocessor->end_function) {
00800       ierr = (preprocessor->end_function)
00801         (type,do_not_keep_original,gctx,ctx,next_problem,problem,next_solution,solution); CHKERRQ(ierr);
00802     } else if (GlobalInfo->solutioncopy && next_solution) {
00803       ierr = (GlobalInfo->solutioncopy)(next_solution,solution); CHKERRQ(ierr);
00804     }
00805     ierr = (GlobalInfo->solutiondelete)(next_solution); CHKERRQ(ierr);
00806   } else solution = next_solution;
00807 
00808   if (next_problem!=problem) {
00809     if (!GlobalInfo->problemdelete)  SETERRQ1(1,"Preprocessor %s requires problem destructor",pclassname);
00810     ierr = (GlobalInfo->problemdelete)(next_problem); CHKERRQ(ierr);
00811   }
00812   ierr = SysProProblemDeleteContext(next_problem); CHKERRQ(ierr);
00813   if (!do_not_keep_original) { /* deallocate next_problem */ }
00814   if (GlobalInfo->errortracer) {
00815     char caption[200];
00816     sprintf(caption,"After preprocessor <%s>=<%s>",pclassname,type);
00817     ierr = (GlobalInfo->errortracer)(problem,solution,caption); CHKERRQ(ierr);
00818   }
00819   *rsolution  = solution;
00820   PetscFunctionReturn(0);
00821 }
00822 
00823 #undef __FUNCT__
00824 #define __FUNCT__ "ChooseFirstTransform"
00825 static PetscErrorCode ChooseFirstTransform
00826 (NumericalProblem problem,const char *preprocess,
00827  PetscTruth user_choices,PetscTruth exhaustive,
00828  SalsaTransform transform,SystemPreprocessor preprocessor,
00829  SalsaTransformObject *transformitem,PetscTruth *moretransform,
00830  const char **type)
00831 {
00832   PetscErrorCode ierr;
00833   PetscFunctionBegin;
00834   if (user_choices || exhaustive || !preprocessor->intelligence ) {
00835     ierr = TransformGetNextUnmarkedItem
00836       (transform,NULL,transformitem,moretransform); CHKERRQ(ierr);
00837     if (!*moretransform)
00838       SETERRQ1(1,"Can not find any <%s>\n",preprocessor->name);
00839     ierr = TransformObjectGetName(*transformitem,type); CHKERRQ(ierr);
00840     ierr = SysProTraceMessage
00841       ("Preprocessor choice: %s=%s\n",preprocess,*type); CHKERRQ(ierr);
00842   } else {
00843     const char *reason;
00844     /* 
00845      * Intelligent setting:
00846      * Compute required categories
00847      */
00848     *moretransform = PETSC_TRUE;
00849     if (GlobalInfo->computecategory) {
00850       char *c;
00851       c = strtok((char*)preprocessor->required,",");
00852       while (c) {
00853         ierr = (GlobalInfo->computecategory)(c,problem); CHKERRQ(ierr);
00854         c = strtok(NULL,",");
00855       }
00856     }
00857     ierr = (preprocessor->intelligence)(problem,type,&reason); CHKERRQ(ierr);
00858     ierr = SysProTraceMessage
00859         ("Intelligent choice <%s=%s> because <%s>\n",preprocess,*type,reason); CHKERRQ(ierr);
00860   }
00861   CHKMEMQ;
00862   PetscFunctionReturn(0);
00863 }
00864 
00865 #undef __FUNCT__
00866 #define __FUNCT__ "PreprocessedSolution"
00867 /*!
00868   This routine handles the application of one preprocessor. Depending on
00869   the runtime setup (see section \ref modes), one choice is applied, or
00870   a sequence of choices is applied consecutively. The forward and backward
00871   transformation of the preprocessor are done here, and if necessary, backup
00872   copies of the system are kept around.
00873  */
00874 PetscErrorCode PreprocessedSolution
00875 (const char *pclassname,NumericalProblem problem,void *prevctx,NumericalSolution *rsolution)
00876 {
00877   NumericalSolution solution=NULL; SystemPreprocessor preprocessor;
00878   SalsaTransform transform; SalsaTransformObject transformitem;
00879   PetscTruth exhaustive,do_not_keep_original,success,moretransform,user_choices;
00880   const char *type;
00881   void *preprocessor_context;
00882   PetscErrorCode ierr;
00883 
00884   PetscFunctionBegin;
00885 
00886   ierr = ReportSysProCallStackState(pclassname); CHKERRQ(ierr);
00887 
00888   ierr = SystemPreprocessorGetByName(pclassname,&preprocessor); CHKERRQ(ierr);
00889   transform = preprocessor->transform; exhaustive = preprocessor->exhaustive;
00890   ierr = TransformGetUserChoices(transform,&user_choices); CHKERRQ(ierr);
00891   CHKMEMQ;
00892 
00893   /*
00894    * Do a global setup that is identical for all preprocessors,
00895    * for instance, compute problem characteristics
00896    */
00897   if (GlobalInfo->classdynamicsetup) {ierr = (GlobalInfo->classdynamicsetup)(pclassname,problem); CHKERRQ(ierr);}
00898 
00899   /*
00900    * Do setup based on the current matrix and preprocessor,
00901    * for instance, disable certain preprocessor choices.
00902    */
00903   ierr = PreprocessorSpecificSetup(pclassname,problem,user_choices); CHKERRQ(ierr);
00904 
00905   /*
00906    * Create a context object that will exist for the duration of
00907    * this preprocessor
00908    */
00909   if (preprocessor->ctxcreate) {
00910     ierr = (preprocessor->ctxcreate)(problem,&preprocessor_context); CHKERRQ(ierr);
00911     ierr = RegisterPreprocessorContext(pclassname,preprocessor_context); CHKERRQ(ierr)
00912   }
00913 
00914   /*
00915    * Pick the first or only choice
00916    */
00917   ierr = ChooseFirstTransform(problem,pclassname,user_choices,exhaustive,transform,preprocessor,
00918                               &transformitem,&moretransform,&type); CHKERRQ(ierr);
00919 
00920   /*
00921    * Do we need to back up the original system?
00922    */
00923   do_not_keep_original = TRUTH(!exhaustive && !GlobalInfo->errortracer);
00924   CHKMEMQ;
00925 
00926   /*
00927    * Loop over preprocessor choices
00928    */
00929   while (moretransform) {
00930     const char *next_solver; int option;
00931     PetscTruth optionsloop,optionflg,solution_is_allocated = PETSC_FALSE;
00932 
00933     optionsloop = exhaustive;
00934 
00935     /*
00936      * Loop over option values; if there are no options, force one run through.
00937      */
00938     ierr = TransformItemGetFirstOption(pclassname,type,&option,&optionflg); CHKERRQ(ierr);
00939     if (!optionflg) {optionflg = PETSC_TRUE; option = 0;}
00940 
00941     while (optionflg) {
00942       void *transform_context; PetscTruth didaniteration=PETSC_TRUE;
00943       NumericalProblem next_problem; NumericalSolution next_solution;
00944 
00945       ierr = SysProTraceMessage("Preprocessor <%s=%s,%d>\n",pclassname,type,option); CHKERRQ(ierr);
00946       ierr = RegisterPreprocessorSetting(pclassname,type,option); CHKERRQ(ierr);
00947 
00948       /*
00949        * Apply the preprocessor
00950        */
00951       solution_is_allocated = PETSC_FALSE;
00952       ierr = SysproPreprocessorStartFunction(preprocessor,type,option,do_not_keep_original,problem,&next_problem,
00953                                              preprocessor_context,&transform_context,&success); CHKERRQ(ierr);
00954       if (!success) {
00955         ierr = SysProTraceMessage("Preprocessor breakdown <%s=%s,%d>\n",pclassname,type,option); CHKERRQ(ierr);
00956         didaniteration = PETSC_FALSE; goto next;
00957       }
00958       ierr = SysProProblemCloneContext(pclassname,type,problem,next_problem); CHKERRQ(ierr);
00959       CHKMEMQ;
00960 
00961       /* 
00962        * Here is where we recursively call the next preprocessor or the final problem solver
00963        */
00964       ierr = SuccessorPreprocessor(pclassname,&next_solver); CHKERRQ(ierr);
00965       preprocesslevel++;
00966       if (next_solver) {
00967         ierr = PreprocessedSolution(next_solver,next_problem,transform_context,&next_solution); CHKERRQ(ierr);
00968       } else if (GlobalInfo->problemsolver) {
00969         ierr = (GlobalInfo->problemsolver)(next_problem,transform_context,&next_solution); CHKERRQ(ierr);
00970       } else SETERRQ(1,"No problem solving routine set");
00971       preprocesslevel--;
00972       CHKMEMQ;
00973 
00974       /*
00975        * Un-preprocess, and delete inner solution
00976        */
00977       ierr = SysProPreprocessorEndFunction
00978         (preprocessor,pclassname,type,do_not_keep_original,preprocessor_context,transform_context,
00979          next_problem,problem,next_solution,&solution); CHKERRQ(ierr);
00980       solution_is_allocated = (PetscTruth)(solution!=next_solution);
00981     next:
00982       CHKMEMQ;
00983       optionflg = optionsloop;
00984       if (optionflg) {ierr = TransformItemGetNextOption(pclassname,type,&option,&optionflg); CHKERRQ(ierr);}
00985       if (didaniteration && optionflg && solution) {
00986         /* if we do another iteration, delete the solution we just created;
00987            the null solution case hapens if the inner solve crashed.
00988         */
00989         if (GlobalInfo->solutiondelete) {
00990           ierr = (GlobalInfo->solutiondelete)(solution); CHKERRQ(ierr);
00991         } else
00992           SETERRQ1(1,"Option looping (for %s) requires solution destructor",pclassname);
00993         if (GlobalInfo->solutioncontextdelete) {
00994           /* ierr = (GlobalInfo->solutioncontextdelete)(next_solution); CHKERRQ(ierr); */
00995         }
00996       }
00997     } /* end of options loops */
00998     CHKMEMQ;
00999 
01000     if (exhaustive) {
01001       ierr = TransformGetNextUnmarkedItem(transform,type,&transformitem,&moretransform); CHKERRQ(ierr);
01002       if (moretransform) {ierr = TransformObjectGetName(transformitem,&type); CHKERRQ(ierr);}
01003     } else moretransform = PETSC_FALSE;
01004     if (solution && moretransform) {
01005       /* if we do another iteration, delete temporary stuff; solution can not exist if the preprocessor broke down */
01006       if (GlobalInfo->solutioncontextdelete) {ierr = (GlobalInfo->solutioncontextdelete)(solution); CHKERRQ(ierr);}
01007       if (GlobalInfo->solutiondelete) {ierr = (GlobalInfo->solutiondelete)(solution); CHKERRQ(ierr);
01008       } else
01009         SETERRQ1(1,"Preprocessor looping (for %s) requires solution destructor",preprocessor);
01010     }
01011   } /* end of preprocessor loop */
01012 
01013   if (preprocessor->unset) {ierr = (preprocessor->unset)(problem); CHKERRQ(ierr);}
01014   if (preprocessor->ctxdelete) {ierr = (preprocessor->ctxdelete)(preprocessor_context); CHKERRQ(ierr);}
01015   CHKMEMQ;
01016 
01017   *rsolution = solution;
01018 
01019   PetscFunctionReturn(0);
01020 }
01021 
01022 #undef __FUNCT__
01023 #define __FUNCT__ "PreprocessedProblemSolving"
01024 /*! Invoking this routine starts the preprocessing and ultimate solution
01025   of the numerical problem.
01026 */
01027 PetscErrorCode PreprocessedProblemSolving
01028 (NumericalProblem problem,NumericalSolution *solution)
01029 {
01030   const char *preprocess; PetscErrorCode ierr;
01031   PetscFunctionBegin;
01032   preprocesslevel = 0;
01033   ierr = GetFirstPreprocessor(&preprocess); CHKERRQ(ierr);
01034   if (!preprocess) {
01035     if (GlobalInfo->problemsolver) {
01036       ierr = (GlobalInfo->problemsolver)(problem,NULL,solution); CHKERRQ(ierr);
01037       if (GlobalInfo->errortracer) {
01038         ierr = (GlobalInfo->errortracer)
01039           (problem,*solution,"After main solve"); CHKERRQ(ierr);
01040       }
01041     } else SETERRQ(1,"No problem solving routine set");
01042   } else {
01043     ierr = PreprocessedSolution
01044       (preprocess,problem,NULL,solution); CHKERRQ(ierr);
01045   }
01046   CHKMEMQ;
01047   PetscFunctionReturn(0);
01048 }