System Preprocessors
linear.c
Go to the documentation of this file.
00001 /*! \file linear.c \anchor linearfile 
00002   
00003   \page linear The Linear package for SysPro
00004 
00005   A linear system is a special case of a numerical problem. This file
00006   contains the routines for creating, deleting, and duplicating linear 
00007   systems.
00008 
00009   \ref lineardef
00010 
00011   \ref lineartransform
00012 
00013   \section lineartransform Transformations for linear systems
00014 
00015   \subpage flipsign
00016 
00017   \subpage singleton
00018 
00019   \subpage distribution
00020 
00021   \subpage scaling
00022 
00023   \subpage approximation
00024 
00025   \subpage ksp
00026 
00027   \subpage pc
00028 
00029   The linear package for SysPro can use the NMD and AnaMod libraries,
00030   but does not require them. Any dependencies on NMD and Anamod should
00031   all be restricted to \ref computec.
00032  */
00033 
00034 #include <stdlib.h>
00035 #include "syspro.h"
00036 #include "syspro_impl.h"
00037 #include "sysprolinear.h"
00038 #include "linear_impl.h"
00039 #include "petscmat.h"
00040 #include "nmd.h"
00041 #include "anamod.h"
00042 
00043 /*! \page lineardef Definition of a linear system
00044 
00045   A linear system has the following components:
00046   - \c A : coefficient matrix
00047   - \c B : a matrix from which to build a preconditioner. Often this
00048   will just be \c A.
00049   - \c rhs : the right hand side
00050   - \c sol : a known solution, if any; there is a boolean to indicate
00051   - \c init : a starting guess for iterative methods
00052   whether a solution is known.
00053   - \c ctx : a void pointer for storing an arbitrary data item; this
00054   can be used by the user.
00055 
00056 
00057   See CreateLinearSystem(), DeleteLinearSystem(), LinearSystemSetParts(),
00058   LinearSystemGetParts(), LinearSystemInheritParts().
00059 
00060   See also \ref linearsolution.
00061  */
00062 
00063 #undef __FUNCT__
00064 #define __FUNCT__ "LinearPackageSetUp"
00065 PetscErrorCode LinearPackageSetUp()
00066 {
00067   PetscFunctionBegin;
00068   PetscFunctionReturn(0);
00069 }
00070 
00071 #undef __FUNCT__
00072 #define __FUNCT__ "CreateLinearSystem"
00073 /*! Allocate the structure for a linear system */
00074 PetscErrorCode CreateLinearSystem(MPI_Comm comm,LinearSystem *system)
00075 {
00076   LinearSystem lnew; PetscErrorCode ierr;
00077   PetscFunctionBegin;
00078   ierr = PetscMalloc(sizeof(struct LinearSystem_),&lnew); CHKERRQ(ierr);
00079   ierr = PetscMemzero(lnew,sizeof(struct LinearSystem_)); CHKERRQ(ierr);
00080   ((NumericalProblem)lnew)->comm = comm;
00081   lnew->cookie = LINSYSCOOKIE; lnew->partsoriginal = 0;
00082   *system = lnew;
00083   CHKMEMQ;
00084   PetscFunctionReturn(0);
00085 }
00086 
00087 #undef __FUNCT__
00088 #define __FUNCT__ "LinearSystemDelete"
00089 PetscErrorCode LinearSystemDelete(LinearSystem system)
00090 {
00091   PetscErrorCode ierr;
00092   PetscFunctionBegin;
00093   SYSPROCHECKVALIDLINSYS(system);
00094   /*printf("deleting parts: %d\n",system->partsoriginal);*/
00095   if (system->A!=system->B && system->B) {
00096     if (system->partsoriginal & 8) {
00097       /*      ierr = PurgeAttachedArrays(system->B); CHKERRQ(ierr);*/
00098       ierr = MatDestroy(system->B); CHKERRQ(ierr); system->B = NULL;}}
00099   if (system->partsoriginal & 16) {
00100     /*    ierr = PurgeAttachedArrays(system->A); CHKERRQ(ierr);*/
00101     ierr = MatDestroy(system->A); CHKERRQ(ierr); system->A = NULL;}
00102   if (system->partsoriginal & 4) {
00103     ierr = VecDestroy(system->Rhs); CHKERRQ(ierr); system->Rhs = NULL;}
00104   if (system->partsoriginal & 2) {
00105     ierr = VecDestroy(system->Sol); CHKERRQ(ierr); system->Sol = NULL;}
00106   if (system->Init && system->partsoriginal & 1) {
00107     ierr = VecDestroy(system->Init); CHKERRQ(ierr); system->Init = NULL;}
00108   if (system->Tmp) {
00109     ierr = VecDestroy(system->Tmp); CHKERRQ(ierr); system->Tmp = NULL;}
00110   CHKMEMQ;
00111   ierr = PetscFree(system); CHKERRQ(ierr);
00112   CHKMEMQ;
00113   PetscFunctionReturn(0);
00114 }
00115 
00116 #undef __FUNCT__
00117 #define __FUNCT__ "LinearDeleteNumericalProblem"
00118 /*! This is like LinearSystemDelete(), except that the argument
00119   has been cast so that this routine can be used as the \c problemdelete
00120   argument of SysProDeclareFunctions().
00121 */
00122 PetscErrorCode LinearDeleteNumericalProblem(NumericalProblem sys)
00123 {
00124   PetscErrorCode ierr;
00125   PetscFunctionBegin;
00126   ierr = LinearSystemDelete((LinearSystem)sys); CHKERRQ(ierr);
00127   CHKMEMQ;
00128   PetscFunctionReturn(0);
00129 }
00130 
00131 #undef __FUNCT__
00132 #define __FUNCT__ "LinearSystemSetParts"
00133 /*! Declare the matrices and vectors for a linear system.
00134 
00135   Arguments:
00136   - \c system
00137   - \c A : the matrix
00138   - \c B : operator to construct the preconditioner from; if NULL, 
00139   (or identical to A), A will be used
00140   - rhs : right hand side
00141   - sol : storage for the computed solution
00142   - init : (optional) nontrivial starting vector for iterative solution
00143  */
00144 PetscErrorCode LinearSystemSetParts
00145 (LinearSystem system,Mat A,Mat B,Vec Rhs,Vec Sol,Vec Init)
00146 {
00147   PetscFunctionBegin;
00148   SYSPROCHECKVALIDLINSYS(system);
00149   if (A && A!=system->A) {
00150     system->A = A; system->partsoriginal |= 16;}
00151   if (B && B!=system->B) {
00152     system->B = B;  system->partsoriginal |= 8;}
00153   if (Rhs && Rhs!=system->Rhs) {
00154 #if PETSC_USE_DEBUG==1
00155     PetscReal nrm;
00156     PetscErrorCode ierr = VecNorm(Rhs,NORM_2,&nrm); CHKERRQ(ierr);
00157     if (nrm==0.) SETERRQ(1,"Setting zero Rhs");
00158 #endif
00159     system->Rhs = Rhs; system->partsoriginal |= 4;}
00160   if (Sol && Sol!=system->Sol) {
00161     system->Sol = Sol; system->partsoriginal |= 2;}
00162   if (Init && Init!=system->Init) {
00163     system->Init = Init; system->partsoriginal |= 1;}
00164   CHKMEMQ;
00165   PetscFunctionReturn(0);
00166 }
00167 
00168 #undef __FUNCT__
00169 #define __FUNCT__ "LinearSystemInheritParts"
00170 /*! Declare the matrices and vectors for a linear system. 
00171   Unlike in LinearSystemSetParts(), here the parts are marked as
00172   not original, so they will not be deleted in DeleteLinearSystem().
00173 */
00174 PetscErrorCode LinearSystemInheritParts
00175 (LinearSystem system,Mat A,Mat B,Vec Rhs,Vec Sol,Vec Init)
00176 {
00177   PetscFunctionBegin;
00178   SYSPROCHECKVALIDLINSYS(system);
00179   if (A && A!=system->A) {
00180     system->A = A; system->partsoriginal &= ALLPARTSNEW-16;}
00181   if (B && B!=system->B) {
00182     system->B = B;  system->partsoriginal &= ALLPARTSNEW-8;}
00183   if (Rhs && Rhs!=system->Rhs) {
00184 #if PETSC_USE_DEBUG==1
00185     PetscReal nrm;
00186     PetscErrorCode ierr = VecNorm(Rhs,NORM_2,&nrm); CHKERRQ(ierr);
00187     if (nrm==0.) SETERRQ(1,"Inheriting zero Rhs");
00188 #endif
00189     system->Rhs = Rhs; system->partsoriginal &= ALLPARTSNEW-4;}
00190   if (Sol && Sol!=system->Sol) {
00191     system->Sol = Sol; system->partsoriginal &= ALLPARTSNEW-2;}
00192   if (Init && Init!=system->Init) {
00193     system->Init = Init; system->partsoriginal &= ALLPARTSNEW-1;}
00194   CHKMEMQ;
00195   PetscFunctionReturn(0);
00196 }
00197 
00198 #undef __FUNCT__
00199 #define __FUNCT__ "LinearSystemGetParts"
00200 /*! Get the matrices and vectors of the system */
00201 PetscErrorCode LinearSystemGetParts
00202 (LinearSystem system,Mat *A,Mat *B,Vec *Rhs,Vec *Sol,Vec *Init)
00203 {
00204   PetscFunctionBegin;
00205   SYSPROCHECKVALIDLINSYS(system);
00206   if (A) *A = system->A;
00207   if (B) *B = system->B;
00208   if (Rhs) *Rhs = system->Rhs;
00209   if (Sol) *Sol = system->Sol;
00210   if (Init) *Init = system->Init;
00211   CHKMEMQ;
00212   PetscFunctionReturn(0);
00213 }
00214 
00215 #undef __FUNCT__
00216 #define __FUNCT__ "LinearSystemSetContext"
00217 PetscErrorCode LinearSystemSetContext(LinearSystem system,void *ctx)
00218 {
00219   PetscFunctionBegin;
00220   SYSPROCHECKVALIDLINSYS(system);
00221   system->ctx = ctx;
00222   PetscFunctionReturn(0);
00223 }
00224 
00225 #undef __FUNCT__
00226 #define __FUNCT__ "LinearSystemGetContext"
00227 PetscErrorCode LinearSystemGetContext(LinearSystem system,void **ctx)
00228 {
00229   PetscFunctionBegin;
00230   SYSPROCHECKVALIDLINSYS(system);
00231   *ctx = system->ctx;
00232   PetscFunctionReturn(0);
00233 }
00234 
00235 #undef __FUNCT__
00236 #define __FUNCT__ "LinearSystemSetKnownSolution"
00237 PetscErrorCode LinearSystemSetKnownSolution(LinearSystem sys,PetscTruth sol)
00238 {
00239   PetscFunctionBegin;
00240   SYSPROCHECKVALIDLINSYS(sys);
00241   sys->known_solution = sol;
00242   PetscFunctionReturn(0);
00243 }
00244 
00245 #undef __FUNCT__
00246 #define __FUNCT__ "LinearSystemGetKnownSolution"
00247 PetscErrorCode LinearSystemGetKnownSolution(LinearSystem sys,PetscTruth *sol)
00248 {
00249   PetscFunctionBegin;
00250   SYSPROCHECKVALIDLINSYS(sys);
00251   *sol = sys->known_solution;
00252   PetscFunctionReturn(0);
00253 }
00254 
00255 #undef __FUNCT__
00256 #define __FUNCT__ "LinearSystemSetMetadata"
00257 PetscErrorCode LinearSystemSetMetadata(LinearSystem system,NMD_metadata nmd)
00258 {
00259   PetscFunctionBegin;
00260   SYSPROCHECKVALIDLINSYS(system);
00261   system->metadata = nmd;
00262   PetscFunctionReturn(0);
00263 }
00264 
00265 #undef __FUNCT__
00266 #define __FUNCT__ "LinearSystemGetMetadata"
00267 PetscErrorCode LinearSystemGetMetadata(LinearSystem system,NMD_metadata *nmd)
00268 {
00269   PetscFunctionBegin;
00270   SYSPROCHECKVALIDLINSYS(system);
00271   *nmd = system->metadata;
00272   PetscFunctionReturn(0);
00273 }
00274 
00275 #undef __FUNCT__
00276 #define __FUNCT__ "LinearSystemGetTmpVector"
00277 PetscErrorCode LinearSystemGetTmpVector(LinearSystem sys,Vec *tmp)
00278 {
00279   PetscErrorCode ierr;
00280   PetscFunctionBegin;
00281   SYSPROCHECKVALIDLINSYS(sys);
00282   if (!sys->Tmp) {ierr = VecDuplicate(sys->Rhs,&(sys->Tmp)); CHKERRQ(ierr);}
00283   *tmp = sys->Tmp;
00284   CHKMEMQ;
00285   PetscFunctionReturn(0);
00286 }
00287 
00288 #undef __FUNCT__
00289 #define __FUNCT__ "LinearSystemDuplicatePointers"
00290 /*! Allocate a new linear system and give it the components of the old
00291   by pointer duplication.
00292 */
00293 PetscErrorCode LinearSystemDuplicatePointers
00294 (LinearSystem problem,LinearSystem *newproblem)
00295 {
00296   LinearSystem lnew; PetscErrorCode ierr;
00297   SYSPROCHECKVALIDLINSYS(problem);
00298   PetscFunctionBegin;
00299   ierr = CreateLinearSystem
00300     (((NumericalProblem)problem)->comm,&lnew); CHKERRQ(ierr);
00301   lnew->partsoriginal = 0;
00302   lnew->A = problem->A; lnew->B = problem->B;
00303   lnew->Rhs = problem->Rhs;
00304   lnew->Sol = problem->Sol;
00305   lnew->Init = problem->Init;
00306   lnew->ctx = problem->ctx;
00307   lnew->known_solution = problem->known_solution;
00308   lnew->metadata = problem->metadata;
00309   *newproblem = lnew;
00310   CHKMEMQ;
00311   PetscFunctionReturn(0);
00312 }
00313 
00314 #undef __FUNCT__
00315 #define __FUNCT__ "LinearSystemDuplicate"
00316 /*! Allocate a new linear system, and create copies in it
00317   of the data structure, but not the values,
00318   of the components of the old system.
00319 
00320   See also LinearSystemCopy().
00321 */
00322 PetscErrorCode LinearSystemDuplicate
00323 (LinearSystem problem,LinearSystem *newproblem)
00324 {
00325   LinearSystem lnew; PetscErrorCode ierr;
00326   PetscFunctionBegin;
00327   SYSPROCHECKVALIDLINSYS(problem);
00328   ierr = CreateLinearSystem
00329     (((NumericalProblem)problem)->comm,&lnew); CHKERRQ(ierr);
00330   lnew->partsoriginal = ALLPARTSNEW;
00331   ierr = MatDuplicate
00332     (problem->A,MAT_DO_NOT_COPY_VALUES,&(lnew->A)); CHKERRQ(ierr);
00333   if (problem->A==problem->B)
00334     lnew->B = lnew->A;
00335   else if (problem->B) {
00336     ierr = MatDuplicate
00337       (problem->B,MAT_DO_NOT_COPY_VALUES,&(lnew->B)); CHKERRQ(ierr);
00338   }
00339   CHKMEMQ;
00340   if (!problem->Rhs)
00341     SETERRQ(1,"Linear system has no rhs");
00342   ierr = VecDuplicate(problem->Rhs,&(lnew->Rhs)); CHKERRQ(ierr);
00343   if (!problem->Sol)
00344     SETERRQ(1,"Linear system has no sol");
00345   ierr = VecDuplicate(problem->Sol,&(lnew->Sol)); CHKERRQ(ierr);
00346   if (problem->Init) {
00347     ierr = VecDuplicate(problem->Init,&(lnew->Init)); CHKERRQ(ierr);
00348   }
00349   lnew->known_solution = problem->known_solution;
00350   lnew->ctx = problem->ctx;
00351 #if 0
00352   {
00353     ierr = NMDCloneObject(problem->nmd,&lnew->nmd); CHKERRQ(ierr);
00354   }
00355 #endif
00356   *newproblem = lnew;
00357   CHKMEMQ;
00358   PetscFunctionReturn(0);
00359 }
00360 
00361 #undef __FUNCT__
00362 #define __FUNCT__ "LinearSystemCopy"
00363 /*! Copy the values of the components of an old linear system into a new.
00364   The new system has to have been created with LinearSystemDuplicate()
00365   because this routine assumes that the data structures are already
00366   in place.
00367 */
00368 PetscErrorCode LinearSystemCopy(LinearSystem old,LinearSystem lnew)
00369 {
00370   PetscErrorCode ierr;
00371   PetscFunctionBegin;
00372   SYSPROCHECKVALIDLINSYSa(old,1);
00373   SYSPROCHECKVALIDLINSYSa(lnew,2);
00374   if (lnew->partsoriginal!=ALLPARTSNEW)
00375     SETERRQ(1,"Attempting in-place copy of linear system");
00376   ierr = MatCopy(old->A,lnew->A,SAME_NONZERO_PATTERN); CHKERRQ(ierr);
00377   if (old->A==old->B)
00378     lnew->B = lnew->A;
00379   else if (old->B) {
00380     ierr = MatCopy(old->B,lnew->B,SAME_NONZERO_PATTERN); CHKERRQ(ierr);
00381   }
00382   CHKMEMQ;
00383   if (old->metadata) {
00384     ierr = NMDCloneObjectStructure
00385       (old->metadata,&lnew->metadata); CHKERRQ(ierr);
00386     ierr = NMDCloneObject(old->metadata,lnew->metadata); CHKERRQ(ierr);
00387   }
00388   ierr = VecCopy(old->Rhs,lnew->Rhs); CHKERRQ(ierr);
00389   ierr = VecCopy(old->Sol,lnew->Sol); CHKERRQ(ierr);
00390   if (old->Init) {
00391     ierr = VecCopy(old->Init,lnew->Init); CHKERRQ(ierr);}
00392   lnew->known_solution = old->known_solution;
00393   CHKMEMQ;
00394   PetscFunctionReturn(0);
00395 }
00396 
00397 /*! \page linearsolution
00398 
00399   There is an object to store the solution of a linear system.
00400 
00401   The solution of a linear system is stored in a data structure that contains
00402   - \c out : the computed output vector
00403   - \c statistics : an NMD object. See LinearSolutionCreateStatistics().
00404 
00405   See CreateLinearSolution(), LinearSolutionDelete(), 
00406   LinearSolutionCopy(), LinearCopyNumericalSolution(), 
00407  */
00408 
00409 #undef __FUNCT__
00410 #define __FUNCT__ "CreateLinearSolution"
00411 PetscErrorCode CreateLinearSolution(LinearSolution *sol)
00412 {
00413   LinearSolution lnew; PetscErrorCode ierr;
00414   PetscFunctionBegin;
00415   ierr = PetscMalloc(sizeof(struct LinearSolution_),&lnew); CHKERRQ(ierr);
00416   ierr = PetscMemzero(lnew,sizeof(struct LinearSolution_)); CHKERRQ(ierr);
00417   ierr = NMDCreateObject(&(lnew->statistics)); CHKERRQ(ierr);
00418   lnew->cookie = LINSOLCOOKIE; *sol = lnew;
00419   CHKMEMQ;
00420   PetscFunctionReturn(0);
00421 }
00422 
00423 #undef __FUNCT__
00424 #define __FUNCT__ "LinearCreateNumericalSolution"
00425 /*!
00426   Shell routine around CreateLinearSolution() to save you some type casting.
00427 
00428   If the first argument is not NULL, its matrix is extracted and used
00429   to create the vector of the solution object.
00430 */
00431 PetscErrorCode LinearCreateNumericalSolution
00432 (NumericalProblem prob,NumericalSolution *sol)
00433 {
00434   LinearSolution linsol; PetscErrorCode ierr; 
00435   PetscFunctionBegin;
00436   ierr = CreateLinearSolution(&linsol); CHKERRQ(ierr);
00437   if (prob) {
00438     MPI_Comm comm; Mat A; Vec V; int localsize;
00439     ierr = LinearSystemGetParts
00440       ((LinearSystem)prob,&A,NULL,NULL,NULL,NULL); CHKERRQ(ierr);
00441     ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr);
00442     ierr = MatGetLocalSize(A,&localsize,NULL); CHKERRQ(ierr);
00443     ierr = VecCreateMPI(comm,localsize,PETSC_DECIDE,&V); CHKERRQ(ierr);
00444     ierr = LinearSolutionSetVector(linsol,V); CHKERRQ(ierr);
00445   }
00446   *sol = (NumericalSolution)linsol;
00447   PetscFunctionReturn(0);
00448 }
00449 
00450 #undef __FUNCT__
00451 #define __FUNCT__ "LinearSolutionDelete"
00452 /*! Delete a linear solution.
00453 
00454   This does not affect the context stored in the solution. That needs
00455   a special purpose routine.
00456 
00457   See also LinearDeleteNumericalSolution().
00458 */
00459 PetscErrorCode LinearSolutionDelete(LinearSolution sol)
00460 {
00461   PetscErrorCode ierr;
00462   PetscFunctionBegin;
00463   SYSPROCHECKVALIDLINSOL(sol);
00464   ierr = VecDestroy(sol->Out); CHKERRQ(ierr);
00465   ierr = NMDDestroyObject(sol->statistics); CHKERRQ(ierr);
00466   ierr = PetscFree(sol); CHKERRQ(ierr);
00467   CHKMEMQ;
00468   PetscFunctionReturn(0);
00469 }
00470 
00471 #undef __FUNCT__
00472 #define __FUNCT__ "LinearDeleteNumericalSolution"
00473 /*! This is like LinearSolutionDelete(), except that the argument
00474   has been cast so that this routine can be used as the \c solutiondelete
00475   argument of SysProDeclareFunctions().
00476 */
00477 PetscErrorCode LinearDeleteNumericalSolution(NumericalSolution sol)
00478 {
00479   PetscErrorCode ierr;
00480   PetscFunctionBegin;
00481   ierr = LinearSolutionDelete((LinearSolution)sol); CHKERRQ(ierr);
00482   CHKMEMQ;
00483   PetscFunctionReturn(0);
00484 }
00485 
00486 #undef __FUNCT__
00487 #define __FUNCT__ "LinearSolutionCopy"
00488 /*! Copy one linear solution object into another.
00489   This clearly only works if their vectors are similarly layed out.
00490 
00491   The context pointer is blindly copied. We may have to think about this 
00492   a bit more.
00493 
00494   See also LinearCopyNumericalSolution().
00495 */
00496 PetscErrorCode LinearSolutionCopy(LinearSolution old,LinearSolution lnew)
00497 {
00498   PetscErrorCode ierr;
00499   PetscFunctionBegin;
00500   SYSPROCHECKVALIDLINSOLa(old,1);
00501   SYSPROCHECKVALIDLINSOLa(lnew,2);
00502   ierr = NMDCloneObject(old->statistics,lnew->statistics); CHKERRQ(ierr);
00503   if (old->Out!=lnew->Out) {ierr = VecCopy(old->Out,lnew->Out); CHKERRQ(ierr);}
00504   lnew->ctx = old->ctx;
00505   CHKMEMQ;
00506   PetscFunctionReturn(0);
00507 }
00508 
00509 #undef __FUNCT__
00510 #define __FUNCT__ "LinearCopyNumericalSolution"
00511 /*! This routine is essentially LinearSolutionCopy(), except that it does
00512   casts of the arguments so that it can be used as the \c solutioncopy
00513   member of SysProDeclareFunctions()
00514 */
00515 PetscErrorCode LinearCopyNumericalSolution(NumericalSolution old,NumericalSolution nnew)
00516 {
00517   PetscErrorCode ierr;
00518   PetscFunctionBegin;
00519   ierr = LinearSolutionCopy
00520     ((LinearSolution)old,(LinearSolution)nnew); CHKERRQ(ierr);
00521   CHKMEMQ;
00522   PetscFunctionReturn(0);
00523 }
00524 
00525 #undef __FUNCT__
00526 #define __FUNCT__ "CreateDefaultLinearSolution"
00527 PetscErrorCode CreateDefaultLinearSolution
00528 (NumericalProblem problem,NumericalSolution *rsol)
00529 {
00530   LinearSystem sys = (LinearSystem)problem;
00531   LinearSolution solution; Vec Rhs,Out;
00532   PetscErrorCode ierr;
00533   PetscFunctionBegin;
00534   SYSPROCHECKVALIDLINSYS(sys);
00535   ierr = CreateLinearSolution(&solution); CHKERRQ(ierr);
00536   ierr = LinearSystemGetParts
00537     (sys,PETSC_NULL,PETSC_NULL,&Rhs,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
00538   ierr = VecDuplicate(Rhs,&Out); CHKERRQ(ierr);
00539   ierr = LinearSolutionSetVector(solution,Out); CHKERRQ(ierr);
00540   *rsol = (NumericalSolution)solution;
00541   CHKMEMQ;
00542   PetscFunctionReturn(0);
00543 }
00544 
00545 #undef __FUNCT__
00546 #define __FUNCT__ "LinearSolutionSetVector"
00547 PetscErrorCode LinearSolutionSetVector(LinearSolution sol,Vec out)
00548 {
00549   PetscFunctionBegin;
00550   SYSPROCHECKVALIDLINSOL(sol);
00551   sol->Out = out;
00552   PetscFunctionReturn(0);
00553 }
00554 
00555 #undef __FUNCT__
00556 #define __FUNCT__ "LinearSolutionGetVector"
00557 PetscErrorCode LinearSolutionGetVector(LinearSolution sol,Vec *out)
00558 {
00559   PetscFunctionBegin;
00560   SYSPROCHECKVALIDLINSOL(sol);
00561   *out = sol->Out;
00562   PetscFunctionReturn(0);
00563 }
00564 
00565 /*! \page solutionstats Solution Statistics
00566   
00567   The \c LinearSolution object carries an \c NMD_metadata object that 
00568   contains performance measurements. This object is initially empty,
00569   so we need to build its content.
00570  */
00571 #undef __FUNCT__
00572 #define __FUNCT__ "LinearSolutionCreateStatistics"
00573 PetscErrorCode LinearSolutionCreateStatistics(LinearSolution sol)
00574 {
00575   NMD_metadata stats; NMD_metadata_category data; PetscErrorCode ierr;
00576   PetscFunctionBegin;
00577   SYSPROCHECKVALIDLINSOL(sol);
00578   stats = sol->statistics;
00579   ierr = NMDObjectAllocateNewCategory(stats,"statistics",&data); CHKERRQ(ierr);
00580   CHKMEMQ;
00581   ierr = NMDCategoryAllocateNewComponent
00582     (data,"converged",NMDInt,NULL); CHKERRQ(ierr);
00583   ierr = NMDCategoryAllocateNewComponent
00584     (data,"abort",NMDInt,NULL); CHKERRQ(ierr);
00585   ierr = NMDCategoryAllocateNewComponent
00586     (data,"iterations",NMDInt,NULL); CHKERRQ(ierr);
00587   ierr = NMDCategoryAllocateNewComponent
00588     (data,"setup-time",NMDReal,NULL); CHKERRQ(ierr);
00589   ierr = NMDCategoryAllocateNewComponent
00590     (data,"solve-time",NMDReal,NULL); CHKERRQ(ierr);
00591   ierr = NMDCategoryAllocateNewComponent
00592     (data,"convergence-type",NMDInt,NULL); CHKERRQ(ierr);
00593   ierr = NMDCategoryAllocateNewComponent
00594     (data,"error",NMDReal,NULL); CHKERRQ(ierr);
00595   CHKMEMQ;
00596   PetscFunctionReturn(0);
00597 }
00598 
00599 #undef __FUNCT__
00600 #define __FUNCT__ "LinearSolutionGetStatistics"
00601 PetscErrorCode LinearSolutionGetStatistics(LinearSolution sol,NMD_metadata *s)
00602 {
00603   PetscFunctionBegin;
00604   SYSPROCHECKVALIDLINSOL(sol);
00605   *s = sol->statistics;
00606   PetscFunctionReturn(0);
00607 }
00608 
00609 #undef __FUNCT__
00610 #define __FUNCT__ "LinearSolutionCopyStats"
00611 PetscErrorCode LinearSolutionCopyStats(LinearSolution in,LinearSolution out)
00612 {
00613   PetscErrorCode ierr;
00614   PetscFunctionBegin;
00615   SYSPROCHECKVALIDLINSOLa(in,1);
00616   SYSPROCHECKVALIDLINSOLa(out,2);
00617   ierr = NMDCloneObject(in->statistics,out->statistics); CHKERRQ(ierr);
00618   CHKMEMQ;
00619   PetscFunctionReturn(0);
00620 }
00621 
00622 /*! \page linearsolutioncontext The linear solution context
00623   We use the context pointer in a LinearSolution object to store diagnostics.
00624   This pointer is blindly copied in LinearSolutionCopy() (unlike the solution 
00625   vector, which is replicated) so we have to be careful with deallocating.
00626  */
00627 
00628 #undef __FUNCT__
00629 #define __FUNCT__ "LinearSolutionSetContext"
00630 PetscErrorCode LinearSolutionSetContext(LinearSolution sol,void *ctx)
00631 {
00632   PetscFunctionBegin;
00633   SYSPROCHECKVALIDLINSOL(sol);
00634   sol->ctx = ctx;
00635   PetscFunctionReturn(0);
00636 }
00637 
00638 #undef __FUNCT__
00639 #define __FUNCT__ "LinearSolutionGetContext"
00640 PetscErrorCode LinearSolutionGetContext(LinearSolution sol,void **ctx)
00641 {
00642   PetscFunctionBegin;
00643   SYSPROCHECKVALIDLINSOL(sol);
00644   *ctx = sol->ctx;
00645   PetscFunctionReturn(0);
00646 }
00647 
00648 #undef __FUNCT__
00649 #define __FUNCT__ "LinearDeleteNumericalSolutionContext"
00650 PetscErrorCode LinearDeleteNumericalSolutionContext(NumericalSolution sol)
00651 {
00652   LinearSolution solution = (LinearSolution)sol;
00653   PetscFunctionBegin;
00654   SYSPROCHECKVALIDLINSOL(solution);
00655   if (solution->ctx) {
00656     /*ierr = delete_diagnostics(solution->ctx); CHKERRQ(ierr);*/
00657   }
00658   PetscFunctionReturn(0);
00659 }
00660 
00661 #undef __FUNCT__
00662 #define __FUNCT__ "LinearSystemTrueDistance"
00663 PetscErrorCode LinearSystemTrueDistance
00664 (LinearSystem system,LinearSolution linsol,PetscReal *rnrm)
00665 {
00666   Vec tmp,out,sol;
00667   PetscReal nrm; PetscErrorCode ierr;
00668   PetscFunctionBegin;
00669   ierr = LinearSystemGetParts
00670     (system,PETSC_NULL,PETSC_NULL,PETSC_NULL,&sol,PETSC_NULL); CHKERRQ(ierr);
00671   ierr = LinearSystemGetTmpVector(system,&tmp); CHKERRQ(ierr);
00672   ierr = LinearSolutionGetVector(linsol,&out); CHKERRQ(ierr);
00673   ierr = VecCopy(sol,tmp); CHKERRQ(ierr);
00674   ierr = VecAXPY(tmp,-1,out); CHKERRQ(ierr);
00675   ierr = VecNorm(tmp,NORM_2,&nrm); CHKERRQ(ierr);
00676   *rnrm = nrm;
00677   CHKMEMQ;
00678   PetscFunctionReturn(0);
00679 }
00680 
00681 #undef __FUNCT__
00682 #define __FUNCT__ "LinearSystemTrueDistancePrint"
00683 PetscErrorCode LinearSystemTrueDistancePrint
00684 (NumericalProblem problem,NumericalSolution solution,char *caption)
00685 {
00686   LinearSystem system = (LinearSystem)problem;
00687   LinearSolution linsol = (LinearSolution)solution;
00688   Mat A; Vec rhs,out,sol,tmp; PetscReal nrm;
00689   PetscTruth known; PetscErrorCode ierr;
00690 
00691   PetscFunctionBegin;
00692   SYSPROCHECKVALIDLINSYS(system);
00693   SYSPROCHECKVALIDLINSOL(linsol);
00694   ierr = LinearSystemGetParts
00695     (system,&A,PETSC_NULL,&rhs,&sol,PETSC_NULL); CHKERRQ(ierr);
00696   ierr = LinearSystemGetTmpVector(system,&tmp); CHKERRQ(ierr);
00697   ierr = LinearSolutionGetVector(linsol,&out); CHKERRQ(ierr);
00698   ierr = MatMult(A,out,tmp); CHKERRQ(ierr);
00699   ierr = VecAXPY(tmp,-1,rhs); CHKERRQ(ierr);
00700   ierr = VecNorm(tmp,NORM_2,&nrm); CHKERRQ(ierr);
00701   CHKMEMQ;
00702   printf("== %s:\n",caption);
00703   printf("Residual 2 norm: %e\n",nrm);
00704   ierr = LinearSystemGetKnownSolution(system,&known); CHKERRQ(ierr);
00705   if (known) {
00706     PetscReal nrm; 
00707     ierr = LinearSystemTrueDistance(system,linsol,&nrm); CHKERRQ(ierr);
00708     printf("Distance to true solution: %e\n",nrm);
00709   }
00710   CHKMEMQ;
00711   PetscFunctionReturn(0);
00712 }
00713 
00714 #undef __FUNCT__
00715 #define __FUNCT__ "PreprocessedLinearSystemSolution"
00716 PetscErrorCode PreprocessedLinearSystemSolution
00717 (LinearSystem sys,LinearSolution *sol)
00718 {
00719   PetscLogDouble *tsetup; PetscErrorCode ierr;
00720   PetscFunctionBegin;
00721   SYSPROCHECKVALIDLINSYS(sys);
00722   ierr = PetscMalloc(sizeof(PetscLogDouble),&tsetup); CHKERRQ(ierr);
00723   ierr = RegisterPreprocessorContext("solution",(void*)tsetup); CHKERRQ(ierr);
00724   CHKMEMQ;
00725   ierr = PreprocessedProblemSolving
00726     ((NumericalProblem)sys,(NumericalSolution*)sol); CHKERRQ(ierr);
00727   CHKMEMQ;
00728   PetscFunctionReturn(0);
00729 }
00730 
00731 #if 0
00732 /*
00733  * ancillary function
00734  */
00735 #undef __FUNCT__
00736 #define __FUNCT__ "PurgeAttachedArrays"
00737 PetscErrorCode PurgeAttachedArrays(Mat A)
00738 {
00739   int icat,icmp; PetscErrorCode ierr;
00740   PetscFunctionBegin;
00741   for (icat=0; icat<ncategories; icat++) {
00742     for (icmp=0; icmp<ncomponents[icat]; icmp++) {
00743       int id=ids[icat][icmp],*isv; PetscScalar *rsv; 
00744       PetscTruth f=PETSC_FALSE;
00745       switch (types[icat][icmp]) {
00746       case ANALYSISINTARRAY :
00747         ierr = PetscObjectComposedDataGetIntstar
00748           ((PetscObject)A,id,isv,f); CHKERRQ(ierr);
00749         if (f) {
00750           ierr = PetscFree(isv); CHKERRQ(ierr);
00751           ierr = PetscObjectComposedDataSetIntstar
00752             ((PetscObject)A,id,0); CHKERRQ(ierr);
00753         }
00754         break;
00755       case ANALYSISDBLARRAY :
00756         ierr = PetscObjectComposedDataGetRealstar
00757           ((PetscObject)A,id,rsv,f); CHKERRQ(ierr);
00758         if (f) {
00759           ierr = PetscFree(rsv); CHKERRQ(ierr);
00760           ierr = PetscObjectComposedDataSetRealstar
00761             ((PetscObject)A,id,0); CHKERRQ(ierr);
00762         }
00763         break;
00764       default : ;
00765       }
00766     }
00767   }
00768   PetscFunctionReturn(0);
00769 }
00770 #endif