System Preprocessors
|
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,PetscBool overwrite, 00025 NumericalProblem inproblem,NumericalProblem *outproblem, 00026 void *gctx,void **ctx,PetscBool *success) 00027 { 00028 LinearSystem problem = (LinearSystem)inproblem, newproblem; 00029 PetscBool 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(PETSC_COMM_SELF,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,PetscBool 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(PETSC_COMM_SELF,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 PetscBool 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 }