System Preprocessors
flipsign.c
Go to the documentation of this file.
00001 /*! \file flipsign.c \ingroup linear
00002  */
00003 
00004 /*! \page flipsign Flip the sign of a matrix
00005 
00006   Most code for iterative methods and preconditioners assumes somewhere that
00007   the sign of a matrix is predominantely positive. Hence, we flip the sign
00008   of matrices that have no positive diagonal elements.
00009  */
00010 #include <stdlib.h>
00011 #include <stdio.h>
00012 #include "petsc.h"
00013 #include "syspro.h"
00014 #include "sysprotransform.h"
00015 #include "sysprolinear.h"
00016 #include "linear_impl.h"
00017 #include "anamod.h"
00018 
00019 #define PREPROCESSOR "flipsign"
00020 
00021 #undef __FUNCT__
00022 #define __FUNCT__ "flipsign"
00023 static PetscErrorCode flipsign
00024 (const char *type,int nopt,PetscTruth overwrite,
00025  NumericalProblem inproblem,NumericalProblem *outproblem,
00026  void *gctx,void **ctx,PetscTruth *success)
00027 {
00028   LinearSystem problem = (LinearSystem)inproblem, newproblem;
00029   PetscTruth flg; PetscErrorCode ierr; 
00030 
00031   PetscFunctionBegin;
00032 
00033   *success = PETSC_TRUE; *(int**)ctx = NULL;
00034   ierr = LinearSystemDuplicatePointers(problem,&newproblem); CHKERRQ(ierr);
00035 
00036   ierr = PetscStrcmp(type,"identity",&flg); CHKERRQ(ierr);
00037   if (flg) {
00038     goto done;
00039   }
00040 
00041 
00042   ierr = PetscStrcmp(type,"flip",&flg); CHKERRQ(ierr);
00043   if (flg) {
00044     Mat A,B,Ause,Buse; Vec X,V,I,Vuse; int ds;
00045     
00046     ierr = LinearSystemGetParts(newproblem,&A,&B,&V,&X,&I); CHKERRQ(ierr);
00047     ierr = SysProComputeQuantity
00048       (inproblem,"variance","diagonal-sign",(void*)&ds,PETSC_NULL,
00049        &flg); CHKERRQ(ierr);
00050     /*
00051      * Flip the sign if necessary
00052      */
00053     if (flg && ds<0) {
00054       ierr = MatDuplicate(A,MAT_COPY_VALUES,&Ause); CHKERRQ(ierr);
00055       ierr = MatScale(Ause,-1); CHKERRQ(ierr);
00056       if (B && B!=A) {
00057         ierr = MatDuplicate(B,MAT_COPY_VALUES,&Buse); CHKERRQ(ierr);
00058         ierr = MatScale(Buse,-1); CHKERRQ(ierr);
00059       } else Buse = Ause;
00060       ierr = VecDuplicate(V,&Vuse); CHKERRQ(ierr);
00061       ierr = VecCopy(V,Vuse); CHKERRQ(ierr);
00062       ierr = VecScale(Vuse,-1); CHKERRQ(ierr);
00063       ierr = LinearSystemSetParts
00064         (newproblem,Ause,Buse,Vuse,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
00065       
00066       {
00067         int* saveone;
00068         ierr = PetscMalloc(sizeof(int),&saveone); CHKERRQ(ierr);
00069         *saveone = 1;
00070         *(int**)ctx = saveone;
00071       }
00072     }
00073     goto done;
00074   }
00075 
00076   SETERRQ1(1,"Unknown flipsign preprocessor <%s>",type);
00077  done:
00078   *outproblem = (NumericalProblem)newproblem;
00079   PetscFunctionReturn(0);
00080 }
00081 
00082 #undef __FUNCT__
00083 #define __FUNCT__ "back_flipsign"
00084 static PetscErrorCode back_flipsign
00085 (const char *flipsign_type,PetscTruth overwrite,void *gctx,void *ctx,
00086  NumericalProblem nextproblem,NumericalProblem problem,
00087  NumericalSolution flipped,NumericalSolution straight)
00088 {
00089   /*LinearSystem sys = (LinearSystem)nextproblem;*/
00090   LinearSolution flippedsol = (LinearSolution)flipped,
00091     straightsol = (LinearSolution)straight;
00092   PetscErrorCode ierr;
00093 
00094   PetscFunctionBegin;
00095 
00096   if (ctx) {
00097     int *saveone = (int*)ctx;
00098     if (*saveone != 1)
00099       SETERRQ1(1,"Context should contain value 1; has %d",*(int*)ctx);
00100     ierr = PetscFree(saveone); CHKERRQ(ierr);
00101   }
00102 
00103   //ierr = DeleteLinearSystem(sys); CHKERRQ(ierr);
00104   ierr = LinearSolutionCopy(flippedsol,straightsol); CHKERRQ(ierr);
00105 
00106   PetscFunctionReturn(0);
00107 }
00108 
00109 #undef __FUNCT__
00110 #define __FUNCT__ "setup_flipsign_choices"
00111 /*! This routine is only called when the flipsign preprocessor
00112   is created by DeclarePreprocessor() inside DeclareFlipsignPreprocessor()
00113 */
00114 static PetscErrorCode setup_flipsign_choices()
00115 {
00116   SalsaTransform tf; SalsaTransformObject cur; PetscErrorCode ierr;
00117 
00118   PetscFunctionBegin;
00119   ierr = TransformGetByName("flipsign",&tf); CHKERRQ(ierr);
00120 
00121   ierr = NewTransformObject(PREPROCESSOR,"identity",&cur); CHKERRQ(ierr);
00122   ierr = TransformObjectSetExplanation
00123     (cur,"do not flip the sign"); CHKERRQ(ierr);
00124 
00125   ierr = NewTransformObject(PREPROCESSOR,"flip",&cur); CHKERRQ(ierr);
00126   ierr = TransformObjectSetExplanation
00127     (cur,"flip the sign of the system"); CHKERRQ(ierr);
00128 
00129   PetscFunctionReturn(0);
00130 }
00131 
00132 #undef __FUNCT__
00133 #define __FUNCT__ "specific_flipsign_choices"
00134 /*!
00135   This is the 'specific setup' phase of the flipsign preprocessor.
00136   See \ref modes for details.
00137 
00138   It disables either the identity or the flip routine, to
00139   leave only the one applicable to this particular system.
00140 
00141 */
00142 static PetscErrorCode specific_flipsign_choices
00143 (NumericalProblem theproblem,SalsaTransform flipsign)
00144 {
00145   SalsaTransformObject to;
00146   PetscTruth flg; int nsign; PetscErrorCode ierr;
00147   PetscFunctionBegin;
00148 
00149   /*
00150    * Mark one or the other; the other is then taken
00151    */
00152   ierr = SysProComputeQuantity
00153     (theproblem,"variance","diagonal-sign",(void*)&nsign,PETSC_NULL,
00154      &flg); CHKERRQ(ierr);
00155   if (flg && nsign<0) {
00156     ierr = TransformObjectGetByName(PREPROCESSOR,"identity",&to); CHKERRQ(ierr);
00157     ierr = TransformObjectMark(to); CHKERRQ(ierr);
00158   } else {
00159     ierr = TransformObjectGetByName(PREPROCESSOR,"flip",&to); CHKERRQ(ierr);
00160     ierr = TransformObjectMark(to); CHKERRQ(ierr);
00161   }
00162 
00163   PetscFunctionReturn(0);
00164 }
00165 
00166 #undef __FUNCT__
00167 #define __FUNCT__ "DeclareFlipsignPreprocessor"
00168 PetscErrorCode DeclareFlipsignPreprocessor(void)
00169 {
00170   PetscErrorCode ierr;
00171   PetscFunctionBegin;
00172   ierr = DeclarePreprocessor
00173     ("flipsign",
00174      &setup_flipsign_choices,&specific_flipsign_choices,0,0,
00175      0,0,
00176      &flipsign,&back_flipsign); CHKERRQ(ierr);
00177   ierr = PreprocessorSetPreservedCategories
00178     (PREPROCESSOR,"laplace,structure,simple"); CHKERRQ(ierr);
00179   PetscFunctionReturn(0);
00180 }