System Preprocessors
|
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(PETSC_COMM_SELF,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(PETSC_COMM_SELF,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,PetscBool 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,PetscBool *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(PETSC_COMM_SELF,1,"Linear system has no rhs"); 00342 ierr = VecDuplicate(problem->Rhs,&(lnew->Rhs)); CHKERRQ(ierr); 00343 if (!problem->Sol) 00344 SETERRQ(PETSC_COMM_SELF,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(PETSC_COMM_SELF,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 PetscBool 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 PetscBool 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