[cig-commits] r6097 - in long/3D/Gale/trunk/src/StgFEM: .
SLE/LinearAlgebra/src
walter at geodynamics.org
walter at geodynamics.org
Fri Feb 23 10:04:01 PST 2007
Author: walter
Date: 2007-02-23 10:03:59 -0800 (Fri, 23 Feb 2007)
New Revision: 6097
Modified:
long/3D/Gale/trunk/src/StgFEM/
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/PETScMatrix.c
Log:
r983 at earth (orig r725): LukeHodkinson | 2007-01-31 19:45:22 -0800
Cleaning up the multigrid code.
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:724
db209038-57f2-0310-97fa-b160e0ae9d04:/trunk:669
+ 38867592-cf10-0410-9e16-a142ea72ac34:/cig:880
db209038-57f2-0310-97fa-b160e0ae9d04:/branches/decomp3d:725
db209038-57f2-0310-97fa-b160e0ae9d04:/trunk:669
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMGSolver.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMGSolver.c 2007-02-23 18:03:57 UTC (rev 6096)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMGSolver.c 2007-02-23 18:03:59 UTC (rev 6097)
@@ -133,6 +133,7 @@
void _PETScMGSolver_Construct( void* matrixSolver, Stg_ComponentFactory* cf, void* data ) {
PETScMGSolver* self = (PETScMGSolver*)matrixSolver;
unsigned nLevels;
+ unsigned nDownIts, nUpIts;
assert( self && Stg_CheckType( self, PETScMGSolver ) );
assert( cf );
@@ -140,7 +141,13 @@
_PETScMatrixSolver_Construct( self, cf, data );
nLevels = Stg_ComponentFactory_GetUnsignedInt( cf, self->name, "levels", 1 );
+ nDownIts = Stg_ComponentFactory_GetUnsignedInt( cf, self->name, "downIterations", 1 );
+ nUpIts = Stg_ComponentFactory_GetUnsignedInt( cf, self->name, "upIterations", 1 );
+
PETScMGSolver_SetLevels( self, nLevels );
+ PETScMGSolver_SetAllDownIterations( self, nDownIts );
+ PETScMGSolver_SetAllUpIterations( self, nUpIts );
+
self->opGen = Stg_ComponentFactory_ConstructByKey( cf, self->name, "opGenerator", MGOpGenerator,
True, data );
MGOpGenerator_SetMatrixSolver( self->opGen, self );
@@ -242,48 +249,47 @@
level->nCycles = 1;
level->R = NULL;
level->P = NULL;
- level->workRes = PETSC_NULL;
- level->workSol = PETSC_NULL;
- level->workRHS = PETSC_NULL;
+ level->A = NULL;
+ level->workRes = NULL;
+ level->workSol = NULL;
+ level->workRHS = NULL;
}
}
-void PETScMGSolver_SetRestriction( void* matrixSolver, unsigned level, void* _R ) {
-#if 0
- PETScMGSolver* self = (PETScMGSolver*)matrixSolver;
- Matrix* R = (Matrix*)_R;
+void PETScMGSolver_SetRestriction( void* matrixSolver, unsigned levelInd, void* R ) {
+ PETScMGSolver* self = (PETScMGSolver*)matrixSolver;
+ PETScMGSolver_Level* level;
assert( self && Stg_CheckType( self, PETScMGSolver ) );
- assert( level < self->nLevels && level > 0 );
- assert( !R || Stg_CheckType( R, Matrix ) );
+ assert( levelInd < self->nLevels && levelInd > 0 );
+ assert( !R || Stg_CheckType( R, PETScMatrix ) );
- 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;
+ level = self->levels + levelInd;
+ if( level->R != R )
+ self->opsChanged = True;
if( R )
Stg_Class_AddRef( R );
-#endif
+ if( level->R )
+ Stg_Class_RemoveRef( level->R );
+ level->R = R;
}
-void PETScMGSolver_SetProlongation( void* matrixSolver, unsigned level, void* _P ) {
-#if 0
- PETScMGSolver* self = (PETScMGSolver*)matrixSolver;
- Matrix* P = (Matrix*)_P;
+void PETScMGSolver_SetProlongation( void* matrixSolver, unsigned levelInd, void* P ) {
+ PETScMGSolver* self = (PETScMGSolver*)matrixSolver;
+ PETScMGSolver_Level* level;
assert( self && Stg_CheckType( self, PETScMGSolver ) );
- assert( level < self->nLevels && level > 0 );
- assert( !P || Stg_CheckType( P, Matrix ) );
+ assert( levelInd < self->nLevels && levelInd > 0 );
+ assert( !P || Stg_CheckType( P, PETScMatrix ) );
- 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;
+ level = self->levels + levelInd;
+ if( level->P != P )
+ self->opsChanged = True;
if( P )
Stg_Class_AddRef( P );
-#endif
+ if( level->P )
+ Stg_Class_RemoveRef( level->P );
+ level->P = P;
}
void PETScMGSolver_SetLevelDownIterations( void* matrixSolver, unsigned level, unsigned nIts ) {
@@ -293,7 +299,7 @@
assert( level < self->nLevels );
self->levels[level].nDownIts = nIts;
- self->solversChanged = False;
+ self->solversChanged = True;
}
void PETScMGSolver_SetLevelUpIterations( void* matrixSolver, unsigned level, unsigned nIts ) {
@@ -304,7 +310,7 @@
assert( level > 0 );
self->levels[level].nUpIts = nIts;
- self->solversChanged = False;
+ self->solversChanged = True;
}
void PETScMGSolver_SetLevelCycles( void* matrixSolver, unsigned level, unsigned nCycles ) {
@@ -315,12 +321,12 @@
assert( level > 0 );
self->levels[level].nCycles = nCycles;
- self->solversChanged = False;
+ self->solversChanged = True;
}
-void PETScMGSolver_SetAllDownIterations( void* matrixSolver, unsigned level, unsigned nIts ) {
+void PETScMGSolver_SetAllDownIterations( void* matrixSolver, unsigned nIts ) {
PETScMGSolver* self = (PETScMGSolver*)matrixSolver;
- unsigned l_i;
+ unsigned l_i;
assert( self && Stg_CheckType( self, PETScMGSolver ) );
@@ -328,9 +334,9 @@
PETScMGSolver_SetLevelDownIterations( self, l_i, nIts );
}
-void PETScMGSolver_SetAllUpIterations( void* matrixSolver, unsigned level, unsigned nIts ) {
+void PETScMGSolver_SetAllUpIterations( void* matrixSolver, unsigned nIts ) {
PETScMGSolver* self = (PETScMGSolver*)matrixSolver;
- unsigned l_i;
+ unsigned l_i;
assert( self && Stg_CheckType( self, PETScMGSolver ) );
@@ -363,14 +369,16 @@
CheckPETScError( ec );
MGOpGenerator_Generate( self->opGen, (Matrix***)&pOps, (Matrix***)&rOps );
+
for( l_i = 1; l_i < self->nLevels; l_i++ ) {
assert( Stg_CheckType( pOps[l_i], PETScMatrix ) );
assert( Stg_CheckType( rOps[l_i], PETScMatrix ) );
- self->levels[l_i].P = pOps[l_i]->petscMat;
+ PETScMGSolver_SetProlongation( self, l_i, pOps[l_i] );
ec = PCMGSetInterpolate( pc, l_i, pOps[l_i]->petscMat );
CheckPETScError( ec );
- self->levels[l_i].R = rOps[l_i]->petscMat;
+
+ PETScMGSolver_SetRestriction( self, l_i, rOps[l_i] );
ec = PCMGSetRestriction( pc, l_i, rOps[l_i]->petscMat );
CheckPETScError( ec );
}
@@ -382,9 +390,6 @@
void PETScMGSolver_UpdateMatrices( PETScMGSolver* self ) {
PETScMGSolver_Level* level;
PC pc;
- Mat* mats;
- Mat downMat, upMat, pMat;
- MatStructure flag;
KSP levelKSP;
PetscErrorCode ec;
unsigned l_i;
@@ -395,55 +400,36 @@
ec = KSPGetPC( self->ksp, &pc );
CheckPETScError( ec );
- mats = AllocArray( Mat, self->nLevels * sizeof(Mat) );
- memset( mats, 0, self->nLevels * sizeof(Mat) );
- mats[self->nLevels - 1] = ((PETScMatrix*)self->matrix)->petscMat;
-
for( l_i = self->nLevels - 1; l_i < self->nLevels; l_i-- ) {
level = self->levels + l_i;
- if( l_i < self->nLevels - 1 ) {
- ec = MatPtAP( mats[l_i + 1], self->levels[l_i + 1].P, MAT_INITIAL_MATRIX, 1.0, mats + l_i );
- CheckPETScError( ec );
+ if( l_i == self->nLevels - 1 )
+ level->A = self->matrix;
+ else {
+ KillObject( level->A );
+ Matrix_PtAP( self->levels[l_i + 1].A, self->levels[l_i + 1].P, (void**)&level->A );
}
ec = PCMGGetSmootherDown( pc, l_i, &levelKSP );
CheckPETScError( ec );
- ec = KSPGetOperators( levelKSP, &downMat, &pMat, &flag );
+ ec = KSPSetOperators( levelKSP, level->A->petscMat, level->A->petscMat, DIFFERENT_NONZERO_PATTERN );
CheckPETScError( ec );
- ec = KSPSetOperators( levelKSP, mats[l_i], mats[l_i], DIFFERENT_NONZERO_PATTERN );
+ ec = PCMGSetResidual( pc, l_i, PCMGDefaultResidual, level->A->petscMat );
CheckPETScError( ec );
- ec = PCMGSetResidual( pc, l_i, PCMGDefaultResidual, mats[l_i] );
- CheckPETScError( ec );
if( l_i > 0 ) {
PCMGGetSmootherUp( pc, l_i, &levelKSP );
- ec = KSPGetOperators( levelKSP, &upMat, &pMat, &flag );
+ ec = KSPSetOperators( levelKSP, level->A->petscMat, level->A->petscMat, DIFFERENT_NONZERO_PATTERN );
CheckPETScError( ec );
- ec = KSPSetOperators( levelKSP, mats[l_i], mats[l_i], DIFFERENT_NONZERO_PATTERN );
- CheckPETScError( ec );
}
-
- if( l_i < self->nLevels - 1 ) {
- if( downMat ) {
- ec = MatDestroy( downMat );
- CheckPETScError( ec );
- }
- if( l_i > 0 && upMat && upMat != downMat ) {
- ec = MatDestroy( upMat );
- CheckPETScError( ec );
- }
- }
}
-
- FreeArray( mats );
}
void PETScMGSolver_UpdateWorkVectors( PETScMGSolver* self ) {
PETScMGSolver_Level* level;
PC pc;
+ unsigned size;
PetscErrorCode ec;
- PetscInt nr, nc;
unsigned l_i;
assert( self && Stg_CheckType( self, PETScMGSolver ) );
@@ -457,60 +443,34 @@
for( l_i = 0; l_i < self->nLevels; l_i++ ) {
level = self->levels + l_i;
- if( l_i == self->nLevels - 1 ) {
- if( level->workRes ) {
- ec = VecDestroy( level->workRes );
- CheckPETScError( ec );
- }
- ec = VecDuplicate( ((PETScVector*)self->curSolution)->petscVec, &level->workRes );
+ Matrix_GetLocalSize( level->A, &size, NULL );
+
+ if( l_i > 0 && (!level->workRes || Vector_GetLocalSize( level->workRes ) != size) ) {
+ if( level->workRes )
+ FreeObject( level->workRes );
+ Vector_Duplicate( self->curSolution, (void**)&level->workRes );
+ Vector_SetLocalSize( level->workRes, size );
+ ec = PCMGSetR( pc, l_i, level->workRes->petscVec );
CheckPETScError( ec );
- ec = PCMGSetR( pc, l_i, self->levels[l_i].workRes );
- CheckPETScError( ec );
}
- else {
- ec = MatGetLocalSize( self->levels[l_i + 1].P, &nr, &nc );
- CheckPETScError( ec );
- if( l_i > 0 ) {
- if( level->workRes ) {
- ec = VecDestroy( level->workRes );
- CheckPETScError( ec );
- }
- ec = VecCreate( self->comm, &level->workRes );
+ if( l_i < self->nLevels - 1 ) {
+ if( !level->workSol || Vector_GetLocalSize( level->workSol ) != size ) {
+ if( level->workSol )
+ FreeObject( level->workSol );
+ Vector_Duplicate( self->curSolution, (void**)&level->workSol );
+ Vector_SetLocalSize( level->workSol, size );
+ ec = PCMGSetX( pc, l_i, level->workSol->petscVec );
CheckPETScError( ec );
- ec = VecSetSizes( level->workRes, nc, PETSC_DECIDE );
- CheckPETScError( ec );
- ec = VecSetFromOptions( level->workRes );
- CheckPETScError( ec );
- ec = PCMGSetR( pc, l_i, level->workRes );
- CheckPETScError( ec );
}
- else {
- if( level->workSol ) {
- ec = VecDestroy( level->workSol );
- CheckPETScError( ec );
- }
- ec = VecCreate( self->comm, &level->workSol );
- CheckPETScError( ec );
- ec = VecSetSizes( level->workSol, nc, PETSC_DECIDE );
- CheckPETScError( ec );
- ec = VecSetFromOptions( level->workSol );
- CheckPETScError( ec );
- ec = PCMGSetX( pc, l_i, level->workSol );
- CheckPETScError( ec );
- if( level->workRHS ) {
- ec = VecDestroy( level->workRHS );
- CheckPETScError( ec );
- }
- ec = VecCreate( self->comm, &level->workRHS );
+ if( !level->workRHS || Vector_GetLocalSize( level->workRHS ) != size ) {
+ if( level->workRHS )
+ FreeObject( level->workRHS );
+ Vector_Duplicate( self->curSolution, (void**)&level->workRHS );
+ Vector_SetLocalSize( level->workRHS, size );
+ ec = PCMGSetRhs( pc, l_i, level->workRHS->petscVec );
CheckPETScError( ec );
- ec = VecSetSizes( level->workRHS, nc, PETSC_DECIDE );
- CheckPETScError( ec );
- ec = VecSetFromOptions( level->workRHS );
- CheckPETScError( ec );
- ec = PCMGSetRhs( pc, l_i, level->workRHS );
- CheckPETScError( ec );
}
}
}
@@ -590,9 +550,16 @@
Stg_Class_RemoveRef( level->R );
if( level->P )
Stg_Class_RemoveRef( level->P );
+ if( level->A )
+ Stg_Class_RemoveRef( level->A );
+
+ FreeObject( level->workRes );
+ FreeObject( level->workSol );
+ FreeObject( level->workRHS );
}
KillArray( self->levels );
self->nLevels = 0;
self->solversChanged = True;
+ self->opsChanged = True;
}
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMGSolver.h
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMGSolver.h 2007-02-23 18:03:57 UTC (rev 6096)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMGSolver.h 2007-02-23 18:03:59 UTC (rev 6097)
@@ -52,12 +52,13 @@
unsigned nUpIts;
unsigned nCycles;
- Mat R;
- Mat P;
+ PETScMatrix* R;
+ PETScMatrix* P;
+ PETScMatrix* A;
- Vec workRes;
- Vec workSol;
- Vec workRHS;
+ PETScVector* workRes;
+ PETScVector* workSol;
+ PETScVector* workRHS;
} PETScMGSolver_Level;
#define __PETScMGSolver \
@@ -109,19 +110,13 @@
*/
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_SetRestriction( void* matrixSolver, unsigned levelInd, void* R );
+ void PETScMGSolver_SetProlongation( void* matrixSolver, unsigned levelInd, void* P );
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 );
+ void PETScMGSolver_SetAllDownIterations( void* matrixSolver, unsigned nIts );
+ void PETScMGSolver_SetAllUpIterations( void* matrixSolver, unsigned nIts );
unsigned PETScMGSolver_GetNumLevels( void* matrixSolver );
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrix.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrix.c 2007-02-23 18:03:57 UTC (rev 6096)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrix.c 2007-02-23 18:03:59 UTC (rev 6097)
@@ -489,7 +489,7 @@
if( *dstMatrix )
ec = MatPtAP( self->petscMat, P->petscMat, MAT_REUSE_MATRIX, 1.0, &(*dstMatrix)->petscMat );
else {
- Matrix_Duplicate( self, dstMatrix );
+ Matrix_Duplicate( self, (void**)dstMatrix );
if( (*dstMatrix)->petscMat )
MatDestroy( (*dstMatrix)->petscMat );
ec = MatPtAP( self->petscMat, P->petscMat, MAT_INITIAL_MATRIX, 1.0, &(*dstMatrix)->petscMat );
More information about the cig-commits
mailing list