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