00001
00002
00003
00004
00005
00006
00007
00008
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 (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
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 (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
00112
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
00136
00137
00138
00139
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
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 }