[cig-commits] r5762 - in long/3D/Gale/trunk/src/StgFEM: .
SLE/LinearAlgebra/src
walter at geodynamics.org
walter at geodynamics.org
Thu Jan 11 21:03:53 PST 2007
Author: walter
Date: 2007-01-11 21:03:52 -0800 (Thu, 11 Jan 2007)
New Revision: 5762
Added:
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMGSolver.c
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMGSolver.h
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMGSolver.meta
Modified:
long/3D/Gale/trunk/src/StgFEM/
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/Init.c
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/LinearAlgebra.h
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MGOpGenerator.c
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MGOpGenerator.h
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MultigridSolver.c
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrix.c
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrixSolver.c
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrixSolver.h
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/SROpGenerator.c
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/SROpGenerator.h
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/types.h
Log:
r946 at earth (orig r705): LukeHodkinson | 2007-01-09 23:16:27 -0800
Adding a multigrid solver that uses PETSc's
multigrid interface.
Property changes on: long/3D/Gale/trunk/src/StgFEM
___________________________________________________________________
Name: svk:merge
- 38867592-cf10-0410-9e16-a142ea72ac34:/cig:880
db209038-57f2-0310-97fa-b160e0ae9d04:/branches/decomp3d:704
db209038-57f2-0310-97fa-b160e0ae9d04:/trunk:669
+ 38867592-cf10-0410-9e16-a142ea72ac34:/cig:880
db209038-57f2-0310-97fa-b160e0ae9d04:/branches/decomp3d:705
db209038-57f2-0310-97fa-b160e0ae9d04:/trunk:669
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/Init.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/Init.c 2007-01-12 05:03:46 UTC (rev 5761)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/Init.c 2007-01-12 05:03:52 UTC (rev 5762)
@@ -70,10 +70,13 @@
PETScMatrix_Type, "0", (Stg_Component_DefaultConstructorFunction*)PETScMatrix_New );
Stg_ComponentRegister_Add( Stg_ComponentRegister_Get_ComponentRegister(),
PETScMatrixSolver_Type, "0", (Stg_Component_DefaultConstructorFunction*)PETScMatrixSolver_New );
+ Stg_ComponentRegister_Add( Stg_ComponentRegister_Get_ComponentRegister(),
+ PETScMGSolver_Type, "0", (Stg_Component_DefaultConstructorFunction*)PETScMGSolver_New );
RegisterParent( PETScVector_Type, Vector_Type );
RegisterParent( PETScMatrix_Type, Matrix_Type );
RegisterParent( PETScMatrixSolver_Type, MatrixSolver_Type );
+ RegisterParent( PETScMGSolver_Type, PETScMatrixSolver_Type );
{
PetscErrorCode ec;
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/LinearAlgebra.h
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/LinearAlgebra.h 2007-01-12 05:03:46 UTC (rev 5761)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/LinearAlgebra.h 2007-01-12 05:03:52 UTC (rev 5762)
@@ -66,10 +66,12 @@
#include <petscvec.h>
#include <petscmat.h>
#include <petscksp.h>
+ #include <petscmg.h>
#include "PETScErrorChecking.h"
#include "PETScVector.h"
#include "PETScMatrix.h"
#include "PETScMatrixSolver.h"
+ #include "PETScMGSolver.h"
#endif
#include "Init.h"
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MGOpGenerator.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MGOpGenerator.c 2007-01-12 05:03:46 UTC (rev 5761)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MGOpGenerator.c 2007-01-12 05:03:52 UTC (rev 5762)
@@ -55,7 +55,7 @@
self = (MGOpGenerator*)_Stg_Component_New( STG_COMPONENT_PASSARGS );
/* Virtual info */
- self->setMultigridSolverFunc = setMultigridSolverFunc;
+ self->setNumLevelsFunc = setNumLevelsFunc;
self->hasExpiredFunc = hasExpiredFunc;
self->generateFunc = generateFunc;
@@ -68,7 +68,7 @@
void _MGOpGenerator_Init( MGOpGenerator* self ) {
assert( self && Stg_CheckType( self, MGOpGenerator ) );
- self->solver = NULL;
+ self->nLevels = 0;
}
@@ -118,14 +118,12 @@
void _MGOpGenerator_Destroy( void* mgOpGenerator, void* data ) {
}
-void _MGOpGenerator_SetMultigridSolver( void* mgOpGenerator, void* _multigridSolver ) {
+void _MGOpGenerator_SetNumLevels( void* mgOpGenerator, unsigned nLevels ) {
MGOpGenerator* self = (MGOpGenerator*)mgOpGenerator;
- MultigridSolver* multigridSolver = (MultigridSolver*)_multigridSolver;
assert( self && Stg_CheckType( self, MGOpGenerator ) );
- assert( !multigridSolver || Stg_CheckType( multigridSolver, MultigridSolver ) );
- self->solver = multigridSolver;
+ self->nLevels = nLevels;
}
@@ -133,12 +131,22 @@
** Public Functions
*/
-MultigridSolver* MGOpGenerator_GetMultigridSolver( void* mgOpGenerator ) {
+void MGOpGenerator_SetMatrixSolver( void* mgOpGenerator, void* _solver ) {
+ MGOpGenerator* self = (MGOpGenerator*)mgOpGenerator;
+ MatrixSolver* solver = (MatrixSolver*)_solver;
+
+ assert( self && Stg_CheckType( self, MGOpGenerator ) );
+ assert( solver && Stg_CheckType( solver, MatrixSolver ) );
+
+ self->solver = solver;
+}
+
+unsigned MGOpGenerator_GetNumLevels( void* mgOpGenerator ) {
MGOpGenerator* self = (MGOpGenerator*)mgOpGenerator;
assert( self && Stg_CheckType( self, MGOpGenerator ) );
- return self->solver;
+ return self->nLevels;
}
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MGOpGenerator.h
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MGOpGenerator.h 2007-01-12 05:03:46 UTC (rev 5761)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MGOpGenerator.h 2007-01-12 05:03:52 UTC (rev 5762)
@@ -45,22 +45,23 @@
extern const Type MGOpGenerator_Type;
/** Virtual function types */
- typedef void (MGOpGenerator_SetMultigridSolverFunc)( void* mgOpGenerator, void* multigridSolver );
+ typedef void (MGOpGenerator_SetNumLevelsFunc)( void* mgOpGenerator, unsigned nLevels );
typedef Bool (MGOpGenerator_HasExpiredFunc)( void* mgOpGenerator );
- typedef void (MGOpGenerator_GenerateFunc)( void* mgOpGenerator );
+ typedef void (MGOpGenerator_GenerateFunc)( void* mgOpGenerator, Matrix*** pOps, Matrix*** rOps );
/** MGOpGenerator class contents */
- #define __MGOpGenerator \
- /* General info */ \
- __Stg_Component \
- \
- /* Virtual info */ \
- MGOpGenerator_SetMultigridSolverFunc* setMultigridSolverFunc; \
- MGOpGenerator_HasExpiredFunc* hasExpiredFunc; \
- MGOpGenerator_GenerateFunc* generateFunc; \
- \
- /* MGOpGenerator info */ \
- MultigridSolver* solver;
+ #define __MGOpGenerator \
+ /* General info */ \
+ __Stg_Component \
+ \
+ /* Virtual info */ \
+ MGOpGenerator_SetNumLevelsFunc* setNumLevelsFunc; \
+ MGOpGenerator_HasExpiredFunc* hasExpiredFunc; \
+ MGOpGenerator_GenerateFunc* generateFunc; \
+ \
+ /* MGOpGenerator info */ \
+ MatrixSolver* solver; \
+ unsigned nLevels;
struct MGOpGenerator { __MGOpGenerator };
@@ -70,14 +71,14 @@
#define MGOPGENERATOR_DEFARGS \
STG_COMPONENT_DEFARGS, \
- MGOpGenerator_SetMultigridSolverFunc* setMultigridSolverFunc, \
+ MGOpGenerator_SetNumLevelsFunc* setNumLevelsFunc, \
MGOpGenerator_HasExpiredFunc* hasExpiredFunc, \
MGOpGenerator_GenerateFunc* generateFunc
#define MGOPGENERATOR_PASSARGS \
STG_COMPONENT_PASSARGS, \
- setMultigridSolverFunc, \
+ setNumLevelsFunc, \
hasExpiredFunc, \
generateFunc
@@ -97,22 +98,23 @@
void _MGOpGenerator_Execute( void* mgOpGenerator, void* data );
void _MGOpGenerator_Destroy( void* mgOpGenerator, void* data );
- void _MGOpGenerator_SetMultigridSolver( void* mgOpGenerator, void* multigridSolver );
+ void _MGOpGenerator_SetNumLevels( void* mgOpGenerator, unsigned nLevels );
- #define MGOpGenerator_SetMultigridSolver( self, multigridSolver ) \
- VirtualCall( self, setMultigridSolverFunc, self, multigridSolver )
+ #define MGOpGenerator_SetNumLevels( self, nLevels ) \
+ VirtualCall( self, setNumLevelsFunc, self, nLevels )
- #define MGOpGenerator_HasExpired( self ) \
+ #define MGOpGenerator_HasExpired( self ) \
VirtualCall( self, hasExpiredFunc, self )
- #define MGOpGenerator_Generate( self ) \
- VirtualCall( self, generateFunc, self )
+ #define MGOpGenerator_Generate( self, pOps, rOps ) \
+ VirtualCall( self, generateFunc, self, pOps, rOps )
/*--------------------------------------------------------------------------------------------------------------------------
** Public functions
*/
- MultigridSolver* MGOpGenerator_GetMultigridSolver( void* mgOpGenerator );
+ void MGOpGenerator_SetMatrixSolver( void* mgOpGenerator, void* solver );
+ unsigned MGOpGenerator_GetNumLevels( void* mgOpGenerator );
/*--------------------------------------------------------------------------------------------------------------------------
** Private Member functions
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MultigridSolver.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MultigridSolver.c 2007-01-12 05:03:46 UTC (rev 5761)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MultigridSolver.c 2007-01-12 05:03:52 UTC (rev 5762)
@@ -156,7 +156,8 @@
MultigridSolver_SetLevels( self, nLevels );
self->opGen = Stg_ComponentFactory_ConstructByKey( cf, self->name, "opGenerator", MGOpGenerator,
True, data );
- MGOpGenerator_SetMultigridSolver( self->opGen, self );
+ MGOpGenerator_SetMatrixSolver( self->opGen, self );
+ MGOpGenerator_SetNumLevels( self->opGen, nLevels );
}
void _MultigridSolver_Build( void* matrixSolver, void* data ) {
@@ -240,6 +241,7 @@
Vector* rhs = (Vector*)_rhs;
Vector* solution = (Vector*)_solution;
Bool rebuildOps;
+ Matrix **pOps, **rOps;
unsigned l_i;
assert( self && Stg_CheckType( self, MultigridSolver ) );
@@ -255,8 +257,13 @@
if( self->ready && !rebuildOps )
return;
- if( !self->ready || rebuildOps )
- MGOpGenerator_Generate( self->opGen );
+ MGOpGenerator_Generate( self->opGen, &pOps, &rOps );
+ for( l_i = 1; l_i < self->nLevels; l_i++ ) {
+ MultigridSolver_SetProlongation( self, l_i, pOps[l_i] );
+ MultigridSolver_SetRestriction( self, l_i, pOps[l_i] );
+ }
+ FreeArray( pOps );
+ FreeArray( rOps );
for( l_i = self->nLevels - 1; l_i < self->nLevels; l_i-- ) {
MultigridSolver_Level* level = self->levels + l_i;
@@ -370,7 +377,7 @@
for( l_i = 0; l_i < nLevels; l_i++ ) {
MultigridSolver_Level* level = self->levels + l_i;
- level->nDownIts = 1;
+ level->nDownIts = 5;
level->nCycles = 1;
if( l_i == 0 ) {
@@ -389,12 +396,13 @@
PETScMatrixSolver_SetPCType( level->downSolver,
PETScMatrixSolver_PCType_RedundantLU );
}
+ PETScMatrixSolver_EnableShifting( level->downSolver, True );
#else
self->downSolver = NULL;
#endif
}
else {
- level->nUpIts = 1;
+ level->nUpIts = 5;
/* Create smoother. */
#ifdef HAVE_PETSC
Added: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMGSolver.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMGSolver.c 2007-01-12 05:03:46 UTC (rev 5761)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMGSolver.c 2007-01-12 05:03:52 UTC (rev 5762)
@@ -0,0 +1,407 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** Copyright (C), 2003, Victorian Partnership for Advanced Computing (VPAC) Ltd, 110 Victoria Street, Melbourne, 3053, Australia.
+**
+** Authors:
+** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+** Patrick D. Sunter, Software Engineer, VPAC. (pds at vpac.org)
+** Luke J. Hodkinson, Computational Engineer, VPAC. (lhodkins at vpac.org)
+** Siew-Ching Tan, Software Engineer, VPAC. (siew at vpac.org)
+** Alan H. Lo, Computational Engineer, VPAC. (alan at vpac.org)
+** Raquibul Hassan, Computational Engineer, VPAC. (raq at vpac.org)
+**
+** This library is free software; you can redistribute it and/or
+** modify it under the terms of the GNU Lesser General Public
+** License as published by the Free Software Foundation; either
+** version 2.1 of the License, or (at your option) any later version.
+**
+** This library is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+** Lesser General Public License for more details.
+**
+** You should have received a copy of the GNU Lesser General Public
+** License along with this library; if not, write to the Free Software
+** Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+**
+** $Id: PETScMGSolver.c 3584 2006-05-16 11:11:07Z PatrickSunter $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <assert.h>
+
+#include <mpi.h>
+#include "StGermain/StGermain.h"
+#include "Discretisation/Discretisation.h"
+#include "LinearAlgebra.h"
+
+
+/* Textual name of this class */
+const Type PETScMGSolver_Type = "PETScMGSolver";
+
+
+/*----------------------------------------------------------------------------------------------------------------------------------
+** Constructors
+*/
+
+PETScMGSolver* PETScMGSolver_New( Name name ) {
+ return _PETScMGSolver_New( sizeof(PETScMGSolver),
+ PETScMGSolver_Type,
+ _PETScMGSolver_Delete,
+ _PETScMGSolver_Print,
+ NULL,
+ (void* (*)(Name))_PETScMGSolver_New,
+ _PETScMGSolver_Construct,
+ _PETScMGSolver_Build,
+ _PETScMGSolver_Initialise,
+ _PETScMGSolver_Execute,
+ _PETScMGSolver_Destroy,
+ name,
+ NON_GLOBAL,
+ PETScMGSolver_SetComm,
+ PETScMatrixSolver_SetMatrix,
+ PETScMatrixSolver_SetMaxIterations,
+ PETScMatrixSolver_SetRelativeTolerance,
+ PETScMatrixSolver_SetAbsoluteTolerance,
+ PETScMatrixSolver_SetUseInitialSolution,
+ PETScMatrixSolver_Solve,
+ PETScMGSolver_Setup,
+ PETScMatrixSolver_GetSolveStatus,
+ PETScMatrixSolver_GetIterations,
+ PETScMatrixSolver_GetMaxIterations,
+ PETScMatrixSolver_GetResidualNorm );
+}
+
+PETScMGSolver* _PETScMGSolver_New( MULTIGRIDSOLVER_DEFARGS ) {
+ PETScMGSolver* self;
+
+ /* Allocate memory */
+ assert( sizeOfSelf >= sizeof(PETScMGSolver) );
+ self = (PETScMGSolver*)_PETScMatrixSolver_New( PETSCMATRIXSOLVER_PASSARGS );
+
+ /* Virtual info */
+
+ /* PETScMGSolver info */
+ _PETScMGSolver_Init( self );
+
+ return self;
+}
+
+void _PETScMGSolver_Init( PETScMGSolver* self ) {
+ assert( self && Stg_CheckType( self, PETScMGSolver ) );
+
+ self->nLevels = 0;
+ self->levels = NULL;
+ self->opGen = NULL;
+ self->ready = False;
+}
+
+
+/*----------------------------------------------------------------------------------------------------------------------------------
+** Virtual functions
+*/
+
+void _PETScMGSolver_Delete( void* matrixSolver ) {
+ PETScMGSolver* self = (PETScMGSolver*)matrixSolver;
+
+ assert( self && Stg_CheckType( self, PETScMGSolver ) );
+
+ PETScMGSolver_Destruct( self );
+
+ /* Delete the parent. */
+ _PETScMatrixSolver_Delete( self );
+}
+
+void _PETScMGSolver_Print( void* matrixSolver, Stream* stream ) {
+ PETScMGSolver* self = (PETScMGSolver*)matrixSolver;
+
+ /* Set the Journal for printing informations */
+ Stream* matrixSolverStream;
+ matrixSolverStream = Journal_Register( InfoStream_Type, "PETScMGSolverStream" );
+
+ assert( self && Stg_CheckType( self, PETScMGSolver ) );
+
+ /* Print parent */
+ Journal_Printf( stream, "PETScMGSolver (ptr): (%p)\n", self );
+ _PETScMatrixSolver_Print( self, stream );
+}
+
+void _PETScMGSolver_Construct( void* matrixSolver, Stg_ComponentFactory* cf, void* data ) {
+ PETScMGSolver* self = (PETScMGSolver*)matrixSolver;
+ unsigned nLevels;
+
+ assert( self && Stg_CheckType( self, PETScMGSolver ) );
+ assert( cf );
+
+ _PETScMatrixSolver_Construct( self, cf, data );
+
+ nLevels = Stg_ComponentFactory_GetUnsignedInt( cf, self->name, "levels", 1 );
+ PETScMGSolver_SetLevels( self, nLevels );
+ self->opGen = Stg_ComponentFactory_ConstructByKey( cf, self->name, "opGenerator", MGOpGenerator,
+ True, data );
+ MGOpGenerator_SetMatrixSolver( self->opGen, self );
+ MGOpGenerator_SetNumLevels( self->opGen, nLevels );
+}
+
+void _PETScMGSolver_Build( void* matrixSolver, void* data ) {
+ PETScMGSolver* self = (PETScMGSolver*)matrixSolver;
+
+ assert( self && Stg_CheckType( self, PETScMGSolver ) );
+
+ if( self->opGen )
+ Build( self->opGen, data, False );
+}
+
+void _PETScMGSolver_Initialise( void* matrixSolver, void* data ) {
+ PETScMGSolver* self = (PETScMGSolver*)matrixSolver;
+
+ assert( self && Stg_CheckType( self, PETScMGSolver ) );
+
+ if( self->opGen )
+ Initialise( self->opGen, data, False );
+}
+
+void _PETScMGSolver_Execute( void* matrixSolver, void* data ) {
+}
+
+void _PETScMGSolver_Destroy( void* matrixSolver, void* data ) {
+}
+
+void PETScMGSolver_SetComm( void* matrixSolver, MPI_Comm comm ) {
+ PETScMatrixSolver* self = (PETScMatrixSolver*)matrixSolver;
+
+ assert( self && Stg_CheckType( self, PETScMatrixSolver ) );
+
+ PETScMatrixSolver_SetComm( self, comm );
+
+ PETScMatrixSolver_SetKSPType( self, PETScMatrixSolver_KSPType_Richardson );
+ PETScMatrixSolver_SetPCType( self, PETScMatrixSolver_PCType_Multigrid );
+}
+
+void PETScMGSolver_Setup( void* matrixSolver, void* rhs, void* _solution ) {
+ PETScMGSolver* self = (PETScMGSolver*)matrixSolver;
+ PETScVector* solution = (PETScVector*)_solution;
+ Bool rebuildOps;
+ PETScMatrix **pOps, **rOps;
+ KSP downKSP;
+ PC downPC;
+ PC pc;
+ PetscErrorCode ec;
+ unsigned l_i;
+
+ assert( self && Stg_CheckType( self, PETScMGSolver ) );
+ assert( solution && Stg_CheckType( solution, PETScVector ) );
+
+ _MatrixSolver_Setup( self, rhs, solution );
+
+ if( self->opGen )
+ rebuildOps = MGOpGenerator_HasExpired( self->opGen );
+ else
+ rebuildOps = False;
+ if( self->ready && !rebuildOps )
+ return;
+
+ ec = KSPGetPC( self->ksp, &pc );
+ CheckPETScError( ec );
+
+ ec = PCMGSetLevels( pc, self->nLevels, PETSC_NULL );
+ CheckPETScError( ec );
+ ec = PCMGSetType( pc, PC_MG_MULTIPLICATIVE );
+ CheckPETScError( ec );
+ ec = PCMGSetGalerkin( pc );
+ CheckPETScError( ec );
+
+ for( l_i = 0; l_i < self->nLevels; l_i++ ) {
+ PCMGGetSmoother( pc, l_i, &downKSP );
+ KSPSetType( downKSP, KSPRICHARDSON );
+ KSPGetPC( downKSP, &downPC );
+ PCSetType( downPC, PCSOR );
+ }
+
+ MGOpGenerator_Generate( self->opGen, (Matrix***)&pOps, (Matrix***)&rOps );
+ for( l_i = 1; l_i < self->nLevels; l_i++ ) {
+ ec = PCMGSetInterpolate( pc, l_i, pOps[l_i]->petscMat );
+ CheckPETScError( ec );
+ ec = PCMGSetRestriction( pc, l_i, rOps[l_i]->petscMat );
+ CheckPETScError( ec );
+
+ if( l_i == self->nLevels - 1 ) {
+ ec = VecDuplicate( solution->petscVec, &self->levels[l_i].workRes );
+ CheckPETScError( ec );
+ }
+ else {
+ unsigned resSize;
+
+ Matrix_GetLocalSize( pOps[l_i], &resSize, NULL );
+
+ ec = VecCreate( self->comm, &self->levels[l_i].workRes );
+ CheckPETScError( ec );
+ ec = VecSetSizes( self->levels[l_i].workRes, resSize, PETSC_DECIDE );
+ CheckPETScError( ec );
+ ec = VecSetFromOptions( self->levels[l_i].workRes );
+ CheckPETScError( ec );
+ }
+
+ ec = PCMGSetR( pc, l_i, self->levels[l_i].workRes );
+ CheckPETScError( ec );
+ }
+ FreeArray( pOps );
+ FreeArray( rOps );
+
+ self->ready = True;
+}
+
+
+/*--------------------------------------------------------------------------------------------------------------------------
+** Public Functions
+*/
+
+void PETScMGSolver_SetLevels( void* matrixSolver, unsigned nLevels ) {
+ PETScMGSolver* self = (PETScMGSolver*)matrixSolver;
+ unsigned l_i;
+
+ assert( self && Stg_CheckType( self, PETScMGSolver ) );
+
+ PETScMGSolver_DestructLevels( self );
+
+ self->nLevels = nLevels;
+ self->levels = AllocArray( PETScMGSolver_Level, nLevels );
+ for( l_i = 0; l_i < nLevels; l_i++ ) {
+ PETScMGSolver_Level* level = self->levels + l_i;
+
+ level->nDownIts = 1;
+ level->nUpIts = (l_i == 0) ? 0 : 1;
+ level->nCycles = 1;
+ level->R = NULL;
+ level->P = NULL;
+ level->workRes = PETSC_NULL;
+ }
+}
+
+void PETScMGSolver_SetRestriction( void* matrixSolver, unsigned level, void* _R ) {
+ PETScMGSolver* self = (PETScMGSolver*)matrixSolver;
+ Matrix* R = (Matrix*)_R;
+
+ assert( self && Stg_CheckType( self, PETScMGSolver ) );
+ assert( level < self->nLevels && level > 0 );
+ assert( !R || Stg_CheckType( R, Matrix ) );
+
+ if( self->levels[level].R != R )
+ self->ready = False;
+ if( self->levels[level].R )
+ Stg_Class_RemoveRef( self->levels[level].R );
+ self->levels[level].R = R;
+ if( R )
+ Stg_Class_AddRef( R );
+}
+
+void PETScMGSolver_SetProlongation( void* matrixSolver, unsigned level, void* _P ) {
+ PETScMGSolver* self = (PETScMGSolver*)matrixSolver;
+ Matrix* P = (Matrix*)_P;
+
+ assert( self && Stg_CheckType( self, PETScMGSolver ) );
+ assert( level < self->nLevels && level > 0 );
+ assert( !P || Stg_CheckType( P, Matrix ) );
+
+ if( self->levels[level].P != P )
+ self->ready = False;
+ if( self->levels[level].P )
+ Stg_Class_RemoveRef( self->levels[level].P );
+ self->levels[level].P = P;
+ if( P )
+ Stg_Class_AddRef( P );
+}
+
+void PETScMGSolver_SetLevelDownIterations( void* matrixSolver, unsigned level, unsigned nIts ) {
+ PETScMGSolver* self = (PETScMGSolver*)matrixSolver;
+
+ assert( self && Stg_CheckType( self, PETScMGSolver ) );
+ assert( level < self->nLevels );
+
+ self->levels[level].nDownIts = nIts;
+ self->ready = False;
+}
+
+void PETScMGSolver_SetLevelUpIterations( void* matrixSolver, unsigned level, unsigned nIts ) {
+ PETScMGSolver* self = (PETScMGSolver*)matrixSolver;
+
+ assert( self && Stg_CheckType( self, PETScMGSolver ) );
+ assert( level < self->nLevels );
+ assert( level > 0 );
+
+ self->levels[level].nUpIts = nIts;
+ self->ready = False;
+}
+
+void PETScMGSolver_SetLevelCycles( void* matrixSolver, unsigned level, unsigned nCycles ) {
+ PETScMGSolver* self = (PETScMGSolver*)matrixSolver;
+
+ assert( self && Stg_CheckType( self, PETScMGSolver ) );
+ assert( level < self->nLevels );
+ assert( level > 0 );
+
+ self->levels[level].nCycles = nCycles;
+ self->ready = False;
+}
+
+void PETScMGSolver_SetAllDownIterations( void* matrixSolver, unsigned level, unsigned nIts ) {
+ PETScMGSolver* self = (PETScMGSolver*)matrixSolver;
+ unsigned l_i;
+
+ assert( self && Stg_CheckType( self, PETScMGSolver ) );
+
+ for( l_i = 1; l_i < self->nLevels; l_i++ )
+ PETScMGSolver_SetLevelDownIterations( self, l_i, nIts );
+}
+
+void PETScMGSolver_SetAllUpIterations( void* matrixSolver, unsigned level, unsigned nIts ) {
+ PETScMGSolver* self = (PETScMGSolver*)matrixSolver;
+ unsigned l_i;
+
+ assert( self && Stg_CheckType( self, PETScMGSolver ) );
+
+ for( l_i = 1; l_i < self->nLevels; l_i++ )
+ PETScMGSolver_SetLevelUpIterations( self, l_i, nIts );
+}
+
+unsigned PETScMGSolver_GetNumLevels( void* matrixSolver ) {
+ PETScMGSolver* self = (PETScMGSolver*)matrixSolver;
+
+ assert( self && Stg_CheckType( self, PETScMGSolver ) );
+
+ return self->nLevels;
+}
+
+
+/*----------------------------------------------------------------------------------------------------------------------------------
+** Private Functions
+*/
+
+void PETScMGSolver_Destruct( PETScMGSolver* self ) {
+ assert( self && Stg_CheckType( self, PETScMGSolver ) );
+
+ PETScMGSolver_DestructLevels( self );
+ KillObject( self->opGen );
+}
+
+void PETScMGSolver_DestructLevels( PETScMGSolver* self ) {
+ unsigned l_i;
+
+ assert( self && Stg_CheckType( self, PETScMGSolver ) );
+
+ for( l_i = 0; l_i < self->nLevels; l_i++ ) {
+ PETScMGSolver_Level* level = self->levels + l_i;
+
+ if( level->R )
+ Stg_Class_RemoveRef( level->R );
+ if( level->P )
+ Stg_Class_RemoveRef( level->P );
+ }
+
+ KillArray( self->levels );
+ self->nLevels = 0;
+ self->ready = False;
+}
Added: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMGSolver.h
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMGSolver.h 2007-01-12 05:03:46 UTC (rev 5761)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMGSolver.h 2007-01-12 05:03:52 UTC (rev 5762)
@@ -0,0 +1,132 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** Copyright (C), 2003, Victorian Partnership for Advanced Computing (VPAC) Ltd, 110 Victoria Street, Melbourne, 3053, Australia.
+**
+** Authors:
+** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+** Patrick D. Sunter, Software Engineer, VPAC. (pds at vpac.org)
+** Luke J. Hodkinson, Computational Engineer, VPAC. (lhodkins at vpac.org)
+** Siew-Ching Tan, Software Engineer, VPAC. (siew at vpac.org)
+** Alan H. Lo, Computational Engineer, VPAC. (alan at vpac.org)
+** Raquibul Hassan, Computational Engineer, VPAC. (raq at vpac.org)
+**
+** This library is free software; you can redistribute it and/or
+** modify it under the terms of the GNU Lesser General Public
+** License as published by the Free Software Foundation; either
+** version 2.1 of the License, or (at your option) any later version.
+**
+** This library is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+** Lesser General Public License for more details.
+**
+** You should have received a copy of the GNU Lesser General Public
+** License along with this library; if not, write to the Free Software
+** Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+**
+*/
+/** \file
+** Role:
+**
+** Assumptions:
+**
+** Invariants:
+**
+** Comments:
+**
+** $Id: PETScMGSolver.h 3584 2006-05-16 11:11:07Z PatrickSunter $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#ifndef __StgFEM_SLE_LinearAlgebra_PETScMGSolver_h__
+#define __StgFEM_SLE_LinearAlgebra_PETScMGSolver_h__
+
+ /** Textual name of this class */
+ extern const Type PETScMGSolver_Type;
+
+ /** Virtual function types */
+
+ /** PETScMGSolver class contents */
+ typedef struct {
+ unsigned nDownIts;
+ unsigned nUpIts;
+ unsigned nCycles;
+
+ Matrix* R;
+ Matrix* P;
+
+ Vec workRes;
+ } PETScMGSolver_Level;
+
+ #define __PETScMGSolver \
+ /* General info */ \
+ __PETScMatrixSolver \
+ \
+ /* Virtual info */ \
+ \
+ /* PETScMGSolver info */ \
+ unsigned nLevels; \
+ PETScMGSolver_Level* levels; \
+ MGOpGenerator* opGen; \
+ Bool ready;
+
+ struct PETScMGSolver { __PETScMGSolver };
+
+ /*--------------------------------------------------------------------------------------------------------------------------
+ ** Constructors
+ */
+
+ #define PETSCMGSOLVER_DEFARGS \
+ PETSCMATRIXSOLVER_DEFARGS
+
+ #define PETSCMGSOLVER_PASSARGS \
+ PETSCMATRIXSOLVER_PASSARGS
+
+ PETScMGSolver* PETScMGSolver_New( Name name );
+ PETScMGSolver* _PETScMGSolver_New( PETSCMATRIXSOLVER_DEFARGS );
+ void _PETScMGSolver_Init( PETScMGSolver* self );
+
+ /*--------------------------------------------------------------------------------------------------------------------------
+ ** Virtual functions
+ */
+
+ void _PETScMGSolver_Delete( void* matrixSolver );
+ void _PETScMGSolver_Print( void* matrixSolver, Stream* stream );
+ void _PETScMGSolver_Construct( void* matrixSolver, Stg_ComponentFactory* cf, void* data );
+ void _PETScMGSolver_Build( void* matrixSolver, void* data );
+ void _PETScMGSolver_Initialise( void* matrixSolver, void* data );
+ void _PETScMGSolver_Execute( void* matrixSolver, void* data );
+ void _PETScMGSolver_Destroy( void* matrixSolver, void* data );
+
+ void PETScMGSolver_SetComm( void* matrixSolver, MPI_Comm comm );
+ void PETScMGSolver_Setup( void* matrixSolver, void* rhs, void* solution );
+
+ /*--------------------------------------------------------------------------------------------------------------------------
+ ** Public functions
+ */
+
+ void PETScMGSolver_SetLevels( void* matrixSolver, unsigned nLevels );
+ void PETScMGSolver_SetRestriction( void* matrixSolver, unsigned level, void* R );
+ void PETScMGSolver_SetProlongation( void* matrixSolver, unsigned level, void* P );
+ void PETScMGSolver_SetLevelDownSolver( void* matrixSolver, unsigned level, MatrixSolver* solver );
+ void PETScMGSolver_SetLevelDownIterations( void* matrixSolver, unsigned level, unsigned nIts );
+ void PETScMGSolver_SetLevelUpSolver( void* matrixSolver, unsigned level, MatrixSolver* solver );
+ void PETScMGSolver_SetLevelUpIterations( void* matrixSolver, unsigned level, unsigned nIts );
+ void PETScMGSolver_SetLevelCycles( void* matrixSolver, unsigned level, unsigned nCycles );
+ void PETScMGSolver_SetAllDownSolver( void* matrixSolver, MatrixSolver* solver );
+ void PETScMGSolver_SetAllDownIterations( void* matrixSolver, unsigned level, unsigned nIts );
+ void PETScMGSolver_SetAllUpSolver( void* matrixSolver, MatrixSolver* solver );
+ void PETScMGSolver_SetAllUpIterations( void* matrixSolver, unsigned level, unsigned nIts );
+ void PETScMGSolver_SetAllSolver( void* matrixSolver, MatrixSolver* solver );
+ void PETScMGSolver_SetCoarseSolver( void* matrixSolver, MatrixSolver* solver );
+
+ unsigned PETScMGSolver_GetNumLevels( void* matrixSolver );
+
+ /*--------------------------------------------------------------------------------------------------------------------------
+ ** Private Member functions
+ */
+
+ void PETScMGSolver_Destruct( PETScMGSolver* self );
+ void PETScMGSolver_DestructLevels( PETScMGSolver* self );
+
+#endif /* __StgFEM_SLE_LinearAlgebra_PETScMGSolver_h__ */
Added: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMGSolver.meta
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMGSolver.meta 2007-01-12 05:03:46 UTC (rev 5761)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMGSolver.meta 2007-01-12 05:03:52 UTC (rev 5762)
@@ -0,0 +1,15 @@
+<?xml version="1.0"?>
+<!DOCTYPE StGermainData SYSTEM "stgermain.dtd">
+<StGermainData xmlns="http://www.vpac.org/StGermain/XML_IO_Handler/Jun2003">
+
+<param name="Name">PETScMGSolver</param>
+<param name="Organisation">VPAC</param>
+<param name="Project">StgFEM</param>
+<param name="Location">StgFEM/SLE/LinearAlgebra/src/</param>
+<param name="Project Web">https://csd.vpac.org/twiki/bin/view/StgFEM/WebHome</param>
+<param name="Copyright">StGermain Framework. Copyright (C) 2003-2005 VPAC.</param>
+<param name="License">The Gnu Lesser General Public License http://www.gnu.org/licenses/lgpl.html</param>
+<param name="Parent"></param>
+<param name="Description">...</param>
+
+</StGermainData>
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrix.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrix.c 2007-01-12 05:03:46 UTC (rev 5761)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrix.c 2007-01-12 05:03:52 UTC (rev 5762)
@@ -554,14 +554,15 @@
PetscInt nr, nc;
assert( self && Stg_CheckType( self, PETScMatrix ) );
- assert( nRows );
- assert( nColumns );
+ assert( nRows || nColumns );
ec = MatGetLocalSize( self->petscMat, &nr, &nc );
CheckPETScError( ec );
- *nRows = (unsigned)nr;
- *nColumns = (unsigned)nc;
+ if( nRows )
+ *nRows = (unsigned)nr;
+ if( nColumns )
+ *nColumns = (unsigned)nc;
}
void PETScMatrix_GetRow( void* matrix, unsigned rowIndex,
@@ -694,7 +695,17 @@
CheckPETScError( ec );
}
+void PETScMatrix_View( void* matrix ) {
+ PETScMatrix* self = (PETScMatrix*)matrix;
+ PetscErrorCode ec;
+ assert( self && Stg_CheckType( self, PETScMatrix ) );
+
+ ec = MatView( self->petscMat, PETSC_VIEWER_STDOUT_SELF );
+ CheckPETScError( ec );
+}
+
+
/*----------------------------------------------------------------------------------------------------------------------------------
** Private Functions
*/
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrixSolver.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrixSolver.c 2007-01-12 05:03:46 UTC (rev 5761)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrixSolver.c 2007-01-12 05:03:52 UTC (rev 5762)
@@ -68,7 +68,7 @@
PETScMatrixSolver_SetAbsoluteTolerance,
PETScMatrixSolver_SetUseInitialSolution,
PETScMatrixSolver_Solve,
- _MatrixSolver_Setup,
+ PETScMatrixSolver_Setup,
PETScMatrixSolver_GetSolveStatus,
PETScMatrixSolver_GetIterations,
PETScMatrixSolver_GetMaxIterations,
@@ -132,6 +132,8 @@
assert( self && Stg_CheckType( self, PETScMatrixSolver ) );
assert( cf );
+
+ _MatrixSolver_Construct( self, cf, data );
}
void _PETScMatrixSolver_Build( void* matrixSolver, void* data ) {
@@ -158,74 +160,98 @@
KSPDestroy( self->ksp );
ec = KSPCreate( self->comm, &self->ksp );
CheckPETScError( ec );
- ec = KSPSetFromOptions( self->ksp );
- CheckPETScError( ec );
}
void PETScMatrixSolver_SetMatrix( void* matrixSolver, void* _matrix ) {
PETScMatrixSolver* self = (PETScMatrixSolver*)matrixSolver;
PETScMatrix* matrix = (PETScMatrix*)_matrix;
+ PetscErrorCode ec;
assert( self && Stg_CheckType( self, PETScMatrixSolver ) );
assert( matrix && Stg_CheckType( matrix, PETScMatrix ) );
_MatrixSolver_SetMatrix( self, matrix );
- KSPSetOperators( self->ksp, matrix->petscMat, matrix->petscMat, DIFFERENT_NONZERO_PATTERN );
+ ec = KSPSetOperators( self->ksp, matrix->petscMat, matrix->petscMat, DIFFERENT_NONZERO_PATTERN );
+ CheckPETScError( ec );
}
void PETScMatrixSolver_SetMaxIterations( void* matrixSolver, unsigned nIterations ) {
PETScMatrixSolver* self = (PETScMatrixSolver*)matrixSolver;
+ PetscErrorCode ec;
assert( self && Stg_CheckType( self, PETScMatrixSolver ) );
- KSPSetTolerances( self->ksp, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, (PetscInt)nIterations );
+ ec = KSPSetTolerances( self->ksp, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, (PetscInt)nIterations );
+ CheckPETScError( ec );
}
void PETScMatrixSolver_SetRelativeTolerance( void* matrixSolver, double tolerance ) {
PETScMatrixSolver* self = (PETScMatrixSolver*)matrixSolver;
+ PetscErrorCode ec;
assert( self && Stg_CheckType( self, PETScMatrixSolver ) );
- KSPSetTolerances( self->ksp, tolerance, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT );
+ ec = KSPSetTolerances( self->ksp, tolerance, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT );
+ CheckPETScError( ec );
}
void PETScMatrixSolver_SetAbsoluteTolerance( void* matrixSolver, double tolerance ) {
PETScMatrixSolver* self = (PETScMatrixSolver*)matrixSolver;
+ PetscErrorCode ec;
assert( self && Stg_CheckType( self, PETScMatrixSolver ) );
- KSPSetTolerances( self->ksp, PETSC_DEFAULT, tolerance, PETSC_DEFAULT, PETSC_DEFAULT );
+ ec = KSPSetTolerances( self->ksp, PETSC_DEFAULT, tolerance, PETSC_DEFAULT, PETSC_DEFAULT );
+ CheckPETScError( ec );
}
void PETScMatrixSolver_SetUseInitialSolution( void* matrixSolver, Bool state ) {
PETScMatrixSolver* self = (PETScMatrixSolver*)matrixSolver;
+ PetscErrorCode ec;
assert( self && Stg_CheckType( self, PETScMatrixSolver ) );
- KSPSetInitialGuessNonzero( self->ksp, (PetscTruth)state );
+ ec = KSPSetInitialGuessNonzero( self->ksp, (PetscTruth)state );
+ CheckPETScError( ec );
}
void PETScMatrixSolver_Solve( void* matrixSolver, void* _rhs, void* _solution ) {
PETScMatrixSolver* self = (PETScMatrixSolver*)matrixSolver;
PETScVector* rhs = (PETScVector*)_rhs;
PETScVector* solution = (PETScVector*)_solution;
+ PetscErrorCode ec;
assert( self && Stg_CheckType( self, PETScMatrixSolver ) );
assert( solution && Stg_CheckType( solution, PETScVector ) );
assert( rhs && Stg_CheckType( rhs, PETScVector ) );
MatrixSolver_Setup( self, rhs, solution );
- KSPSolve( self->ksp, rhs->petscVec, solution->petscVec );
+ ec = KSPSolve( self->ksp, rhs->petscVec, solution->petscVec );
+ CheckPETScError( ec );
}
+void PETScMatrixSolver_Setup( void* matrixSolver, void* rhs, void* solution ) {
+ PETScMatrixSolver* self = (PETScMatrixSolver*)matrixSolver;
+ PetscErrorCode ec;
+
+ assert( self && Stg_CheckType( self, PETScMatrixSolver ) );
+
+ _MatrixSolver_Setup( self, rhs, solution );
+
+ ec = KSPSetFromOptions( self->ksp );
+ CheckPETScError( ec );
+}
+
MatrixSolver_Status PETScMatrixSolver_GetSolveStatus( void* matrixSolver ) {
PETScMatrixSolver* self = (PETScMatrixSolver*)matrixSolver;
KSPConvergedReason reason;
+ PetscErrorCode ec;
assert( self && Stg_CheckType( self, PETScMatrixSolver ) );
- KSPGetConvergedReason( self->ksp, &reason );
+ ec = KSPGetConvergedReason( self->ksp, &reason );
+ CheckPETScError( ec );
return reason;
}
@@ -233,10 +259,12 @@
unsigned PETScMatrixSolver_GetIterations( void* matrixSolver ) {
PETScMatrixSolver* self = (PETScMatrixSolver*)matrixSolver;
PetscInt nIts;
+ PetscErrorCode ec;
assert( self && Stg_CheckType( self, PETScMatrixSolver ) );
- KSPGetIterationNumber( self->ksp, &nIts );
+ ec = KSPGetIterationNumber( self->ksp, &nIts );
+ CheckPETScError( ec );
return (unsigned)nIts;
}
@@ -244,10 +272,12 @@
unsigned PETScMatrixSolver_GetMaxIterations( void* matrixSolver ) {
PETScMatrixSolver* self = (PETScMatrixSolver*)matrixSolver;
PetscInt nIts;
+ PetscErrorCode ec;
assert( self && Stg_CheckType( self, PETScMatrixSolver ) );
- KSPGetTolerances( self->ksp, PETSC_NULL, PETSC_NULL, PETSC_NULL, &nIts );
+ ec = KSPGetTolerances( self->ksp, PETSC_NULL, PETSC_NULL, PETSC_NULL, &nIts );
+ CheckPETScError( ec );
return (unsigned)nIts;
}
@@ -255,10 +285,12 @@
double PETScMatrixSolver_GetResidualNorm( void* matrixSolver ) {
PETScMatrixSolver* self = (PETScMatrixSolver*)matrixSolver;
PetscScalar norm;
+ PetscErrorCode ec;
assert( self && Stg_CheckType( self, PETScMatrixSolver ) );
- KSPGetResidualNorm( self->ksp, &norm );
+ ec = KSPGetResidualNorm( self->ksp, &norm );
+ CheckPETScError( ec );
return (double)norm;
}
@@ -270,18 +302,22 @@
void PETScMatrixSolver_SetKSPType( void* matrixSolver, PETScMatrixSolver_KSPType type ) {
PETScMatrixSolver* self = (PETScMatrixSolver*)matrixSolver;
+ PetscErrorCode ec;
assert( self && Stg_CheckType( self, PETScMatrixSolver ) );
switch( type ) {
case PETScMatrixSolver_KSPType_Richardson:
- KSPSetType( self->ksp, KSPRICHARDSON );
+ ec = KSPSetType( self->ksp, KSPRICHARDSON );
+ CheckPETScError( ec );
break;
case PETScMatrixSolver_KSPType_GMRes:
- KSPSetType( self->ksp, KSPGMRES );
+ ec = KSPSetType( self->ksp, KSPGMRES );
+ CheckPETScError( ec );
break;
case PETScMatrixSolver_KSPType_PreOnly:
- KSPSetType( self->ksp, KSPPREONLY );
+ ec = KSPSetType( self->ksp, KSPPREONLY );
+ CheckPETScError( ec );
break;
default:
assert( 0 );
@@ -291,51 +327,67 @@
void PETScMatrixSolver_SetPCType( void* matrixSolver, PETScMatrixSolver_PCType type ) {
PETScMatrixSolver* self = (PETScMatrixSolver*)matrixSolver;
PC pc;
+ PetscErrorCode ec;
assert( self && Stg_CheckType( self, PETScMatrixSolver ) );
- KSPGetPC( self->ksp, &pc );
+ ec = KSPGetPC( self->ksp, &pc );
+ CheckPETScError( ec );
if( type == PETScMatrixSolver_PCType_Jacobi ) {
- PCSetType( pc, PCJACOBI );
+ ec = PCSetType( pc, PCJACOBI );
+ CheckPETScError( ec );
}
else if( type == PETScMatrixSolver_PCType_BlockJacobi ) {
- PCSetType( pc, PCBJACOBI );
+ ec = PCSetType( pc, PCBJACOBI );
+ CheckPETScError( ec );
}
else if( type == PETScMatrixSolver_PCType_SOR ) {
- PCSetType( pc, PCSOR );
+ ec = PCSetType( pc, PCSOR );
+ CheckPETScError( ec );
}
-#if 0
- else if( type == PETScMatrixSolver_PCType_ParallelSOR ) {
- KSP* subKSPs;
- PC subPC;
-
- PCSetType( pc, BJACOBI );
- PCSetUp( pc );
- PCBJacobiGetSubKSP( pc, PETSC_NULL, PETSC_NULL, &subKSPs );
- KSPGetPC( subKSPs[0], &subPC );
- PCSetType( subPC, PCSOR );
- }
-#endif
else if( type == PETScMatrixSolver_PCType_LU ) {
- PCSetType( pc, PCLU );
+ ec = PCSetType( pc, PCLU );
+ CheckPETScError( ec );
}
else if( type == PETScMatrixSolver_PCType_RedundantLU ) {
- PCSetType( pc, PCREDUNDANT );
- PCRedundantGetPC( pc, &pc );
- PCSetType( pc, PCLU );
+ ec = PCSetType( pc, PCREDUNDANT );
+ CheckPETScError( ec );
+ ec = PCRedundantGetPC( pc, &pc );
+ CheckPETScError( ec );
+ ec = PCSetType( pc, PCLU );
+ CheckPETScError( ec );
}
else if( type == PETScMatrixSolver_PCType_ILU ) {
- PCSetType( pc, PCILU );
+ ec = PCSetType( pc, PCILU );
+ CheckPETScError( ec );
}
+ else if( type == PETScMatrixSolver_PCType_Multigrid ) {
+ ec = PCSetType( pc, PCMG );
+ CheckPETScError( ec );
+ }
else if( type == PETScMatrixSolver_PCType_None ) {
- PCSetType( pc, PCNONE );
+ ec = PCSetType( pc, PCNONE );
+ CheckPETScError( ec );
}
else {
assert( 0 ); /* TODO */
}
}
+void PETScMatrixSolver_EnableShifting( void* matrixSolver, Bool state ) {
+ PETScMatrixSolver* self = (PETScMatrixSolver*)matrixSolver;
+ PC pc;
+ PetscErrorCode ec;
+ assert( self && Stg_CheckType( self, PETScMatrixSolver ) );
+
+ ec = KSPGetPC( self->ksp, &pc );
+ CheckPETScError( ec );
+ ec = PCFactorSetShiftPd( pc, (PetscTruth)state );
+ CheckPETScError( ec );
+}
+
+
/*----------------------------------------------------------------------------------------------------------------------------------
** Private Functions
*/
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrixSolver.h
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrixSolver.h 2007-01-12 05:03:46 UTC (rev 5761)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrixSolver.h 2007-01-12 05:03:52 UTC (rev 5762)
@@ -92,6 +92,7 @@
void PETScMatrixSolver_SetUseInitialSolution( void* matrixSolver, Bool state );
void PETScMatrixSolver_Solve( void* matrixSolver, void* rhs, void* solution );
+ void PETScMatrixSolver_Setup( void* matrixSolver, void* rhs, void* solution );
MatrixSolver_Status PETScMatrixSolver_GetSolveStatus( void* matrixSolver );
unsigned PETScMatrixSolver_GetIterations( void* matrixSolver );
@@ -104,6 +105,7 @@
void PETScMatrixSolver_SetKSPType( void* matrixSolver, PETScMatrixSolver_KSPType type );
void PETScMatrixSolver_SetPCType( void* matrixSolver, PETScMatrixSolver_PCType type );
+ void PETScMatrixSolver_EnableShifting( void* matrixSolver, Bool state );
/*--------------------------------------------------------------------------------------------------------------------------
** Private Member functions
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/SROpGenerator.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/SROpGenerator.c 2007-01-12 05:03:46 UTC (rev 5761)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/SROpGenerator.c 2007-01-12 05:03:52 UTC (rev 5762)
@@ -61,7 +61,7 @@
_SROpGenerator_Destroy,
name,
NON_GLOBAL,
- _MGOpGenerator_SetMultigridSolver,
+ _MGOpGenerator_SetNumLevels,
SROpGenerator_HasExpired,
SROpGenerator_Generate );
}
@@ -151,14 +151,20 @@
return False;
}
-void SROpGenerator_Generate( void* srOpGenerator ) {
+void SROpGenerator_Generate( void* srOpGenerator, Matrix*** pOps, Matrix*** rOps ) {
SROpGenerator* self = (SROpGenerator*)srOpGenerator;
assert( self && Stg_CheckType( self, SROpGenerator ) );
+ assert( pOps && rOps );
+ *pOps = AllocArray( Matrix*, self->nLevels );
+ *rOps = AllocArray( Matrix*, self->nLevels );
+ memset( *pOps, 0, self->nLevels * sizeof(Matrix*) );
+ memset( *rOps, 0, self->nLevels * sizeof(Matrix*) );
+
self->fineEqNum = self->fineVar->eqNum;
SROpGenerator_GenMeshes( self );
- SROpGenerator_GenOps( self );
+ SROpGenerator_GenOps( self, *pOps, *rOps );
SROpGenerator_DestructMeshes( self );
}
@@ -190,7 +196,7 @@
assert( self && Stg_CheckType( self, SROpGenerator ) );
- nLevels = MultigridSolver_GetNumLevels( self->solver );
+ nLevels = self->nLevels;
self->meshes = AllocArray( Mesh*, nLevels );
memset( self->meshes, 0, nLevels * sizeof(Mesh*) );
self->topMaps = AllocArray( unsigned*, nLevels );
@@ -221,7 +227,7 @@
assert( self && Stg_CheckType( self, SROpGenerator ) );
assert( self->meshes );
- assert( level < MultigridSolver_GetNumLevels( self->solver ) );
+ assert( level < self->nLevels );
fMesh = self->meshes[level + 1];
nDims = Mesh_GetDimSize( fMesh );
@@ -266,12 +272,12 @@
assert( self && Stg_CheckType( self, SROpGenerator ) );
assert( self->meshes );
- assert( level < MultigridSolver_GetNumLevels( self->solver ) );
+ assert( level < self->nLevels );
fMesh = self->meshes[level + 1];
cMesh = self->meshes[level];
nDims = Mesh_GetDimSize( fMesh );
- nLevels = MultigridSolver_GetNumLevels( self->solver );
+ nLevels = self->nLevels;
nDomainNodes = Mesh_GetDomainSize( cMesh, MT_VERTEX );
self->topMaps[level] = ReallocArray( self->topMaps[level], unsigned, nDomainNodes );
for( n_i = 0; n_i < nDomainNodes; n_i++ ) {
@@ -319,7 +325,7 @@
assert( self && Stg_CheckType( self, SROpGenerator ) );
assert( self->meshes && self->topMaps && self->eqNums );
- assert( level < MultigridSolver_GetNumLevels( self->solver ) - 1 );
+ assert( level < self->nLevels );
cMesh = self->meshes[level];
nDomainNodes = Mesh_GetDomainSize( cMesh, MT_VERTEX );
@@ -400,16 +406,17 @@
self->eqNumBases[level] = base;
}
-void SROpGenerator_GenOps( SROpGenerator* self ) {
+void SROpGenerator_GenOps( SROpGenerator* self, Matrix** pOps, Matrix** rOps ) {
unsigned nLevels;
Matrix *fineMat, *P;
unsigned nRows, nCols;
unsigned l_i;
assert( self && Stg_CheckType( self, SROpGenerator ) );
+ assert( pOps && rOps );
fineMat = MatrixSolver_GetMatrix( self->solver );
- nLevels = MultigridSolver_GetNumLevels( self->solver );
+ nLevels = self->nLevels;
for( l_i = nLevels - 1; l_i > 0; l_i-- ) {
nRows = self->eqNums[l_i] ? self->nLocalEqNums[l_i] :
@@ -429,8 +436,8 @@
}
SROpGenerator_GenLevelOp( self, l_i, P );
- MultigridSolver_SetRestriction( self->solver, l_i, P );
- MultigridSolver_SetProlongation( self->solver, l_i, P );
+ pOps[l_i] = P;
+ rOps[l_i] = P;
}
}
@@ -449,7 +456,7 @@
assert( self && Stg_CheckType( self, SROpGenerator ) );
assert( self->meshes );
- assert( level < MultigridSolver_GetNumLevels( self->solver ) );
+ assert( level < self->nLevels );
assert( P );
fMesh = self->meshes[level];
@@ -519,7 +526,7 @@
assert( self && Stg_CheckType( self, SROpGenerator ) );
assert( self->meshes );
- assert( level < MultigridSolver_GetNumLevels( self->solver ) );
+ assert( level < self->nLevels );
fMesh = self->meshes[level];
cMesh = self->meshes[level - 1];
@@ -591,7 +598,7 @@
assert( self && Stg_CheckType( self, SROpGenerator ) );
if( self->meshes ) {
- nLevels = MultigridSolver_GetNumLevels( self->solver );
+ nLevels = self->nLevels;
for( l_i = 0; l_i < nLevels - 1; l_i++ ) {
FreeObject( self->meshes[l_i] );
FreeArray( self->topMaps[l_i] );
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/SROpGenerator.h
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/SROpGenerator.h 2007-01-12 05:03:46 UTC (rev 5761)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/SROpGenerator.h 2007-01-12 05:03:52 UTC (rev 5762)
@@ -91,7 +91,7 @@
void _SROpGenerator_Destroy( void* srOpGenerator, void* data );
Bool SROpGenerator_HasExpired( void* srOpGenerator );
- void SROpGenerator_Generate( void* srOpGenerator );
+ void SROpGenerator_Generate( void* srOpGenerator, Matrix*** pOps, Matrix*** rOps );
/*--------------------------------------------------------------------------------------------------------------------------
** Public functions
@@ -108,7 +108,7 @@
void SROpGenerator_GenLevelTopMap( SROpGenerator* self, unsigned level );
void SROpGenerator_GenLevelEqNums( SROpGenerator* self, unsigned level );
- void SROpGenerator_GenOps( SROpGenerator* self );
+ void SROpGenerator_GenOps( SROpGenerator* self, Matrix** pOps, Matrix** rOps );
void SROpGenerator_GenLevelOp( SROpGenerator* self, unsigned level, Matrix* P );
void SROpGenerator_CalcOpNonZeros( SROpGenerator* self, unsigned level,
unsigned** nDiagNonZeros, unsigned** nOffDiagNonZeros );
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/types.h
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/types.h 2007-01-12 05:03:46 UTC (rev 5761)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/types.h 2007-01-12 05:03:52 UTC (rev 5762)
@@ -61,6 +61,7 @@
typedef struct PETScVector PETScVector;
typedef struct PETScMatrix PETScMatrix;
typedef struct PETScMatrixSolver PETScMatrixSolver;
+ typedef struct PETScMGSolver PETScMGSolver;
#endif
typedef enum {
@@ -89,6 +90,7 @@
PETScMatrixSolver_PCType_LU,
PETScMatrixSolver_PCType_RedundantLU,
PETScMatrixSolver_PCType_ILU,
+ PETScMatrixSolver_PCType_Multigrid,
PETScMatrixSolver_PCType_None
} PETScMatrixSolver_PCType;
More information about the cig-commits
mailing list