[cig-commits] r6093 - in long/3D/Gale/trunk/src/StgFEM: .
SLE/LinearAlgebra/src
walter at geodynamics.org
walter at geodynamics.org
Fri Feb 23 10:03:50 PST 2007
Author: walter
Date: 2007-02-23 10:03:49 -0800 (Fri, 23 Feb 2007)
New Revision: 6093
Modified:
long/3D/Gale/trunk/src/StgFEM/
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/Matrix.c
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/Matrix.h
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MatrixSolver.c
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MultigridSolver.c
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MultigridSolver.h
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
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrix.h
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/Vector.c
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/types.h
Log:
r979 at earth (orig r721): LukeHodkinson | 2007-01-29 14:04:23 -0800
Quite a lot of experimental changes directed at
figuring out why multigrid fails when the geometry
of the system deforms.
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:720
db209038-57f2-0310-97fa-b160e0ae9d04:/trunk:669
+ 38867592-cf10-0410-9e16-a142ea72ac34:/cig:880
db209038-57f2-0310-97fa-b160e0ae9d04:/branches/decomp3d:721
db209038-57f2-0310-97fa-b160e0ae9d04:/trunk:669
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/Matrix.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/Matrix.c 2007-02-23 18:03:47 UTC (rev 6092)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/Matrix.c 2007-02-23 18:03:49 UTC (rev 6093)
@@ -109,10 +109,18 @@
*/
void _Matrix_Delete( void* matrix ) {
- Matrix* self = (Matrix*)matrix;
+ Matrix* self = (Matrix*)matrix;
+ MatrixSolver* solver;
+ unsigned s_i;
assert( self && Stg_CheckType( self, Matrix ) );
+ /* Remove from solvers. */
+ for( s_i = 0; s_i < List_GetSize( self->solvers ); s_i++ ) {
+ solver = *List_Get( self->solvers, s_i, MatrixSolver* );
+ solver->matrix = NULL;
+ }
+
/* Delete the parent. */
_Stg_Component_Delete( self );
}
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/Matrix.h
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/Matrix.h 2007-02-23 18:03:47 UTC (rev 6092)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/Matrix.h 2007-02-23 18:03:49 UTC (rev 6093)
@@ -67,8 +67,8 @@
typedef void (Matrix_MultiplyFunc)( void* matrix, void* vector, void* dstVector );
typedef void (Matrix_TransposeMultiplyFunc)( void* matrix, void* vector, void* dstVector );
typedef void (Matrix_MultiplyAddFunc)( void* matrix, void* vector0, void* vector1, void* dstVector );
- typedef void (Matrix_PAPtFunc)( void* matrix, void* P, void* dstMatrix );
- typedef void (Matrix_PtAPFunc)( void* matrix, void* P, void* dstMatrix );
+ typedef void (Matrix_PAPtFunc)( void* matrix, void* P, void** dstMatrix );
+ typedef void (Matrix_PtAPFunc)( void* matrix, void* P, void** dstMatrix );
typedef void (Matrix_MatrixMultiplyFunc)( void* matrix, void* matrix0, void* dstMatrix );
typedef double (Matrix_L2NormFunc)( void* matrix );
typedef void (Matrix_TransposeFunc)( void* matrix, void* dstMatrix );
@@ -309,7 +309,7 @@
#define Matrix_CopyEntries( self, dstMatrix ) \
VirtualCall( self, copyEntriesFunc, self, dstMatrix )
- #define Matrix_Dump( self, filename ) \
+ #define Matrix_Dump( self, filename ) \
VirtualCall( self, dumpFunc, self, filename )
/*--------------------------------------------------------------------------------------------------------------------------
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MatrixSolver.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MatrixSolver.c 2007-02-23 18:03:47 UTC (rev 6092)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MatrixSolver.c 2007-02-23 18:03:49 UTC (rev 6093)
@@ -155,11 +155,15 @@
if( matrix == self->matrix )
return;
- if( self->matrix )
+ if( self->matrix ) {
List_Remove( self->matrix->solvers, &self );
+ Stg_Class_RemoveRef( self->matrix );
+ }
self->matrix = matrix;
- if( matrix )
+ if( matrix ) {
+ Stg_Class_AddRef( matrix );
List_Append( matrix->solvers, &self );
+ }
self->matrixChanged = True;
}
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MultigridSolver.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MultigridSolver.c 2007-02-23 18:03:47 UTC (rev 6092)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MultigridSolver.c 2007-02-23 18:03:49 UTC (rev 6093)
@@ -96,22 +96,16 @@
self->nLevels = 0;
self->levels = NULL;
self->opGen = NULL;
- self->matricesReady = False;
- self->solversReady = False;
+ self->solversChanged = True;
+ self->opsChanged = True;
self->curIt = 0;
self->maxIts = 500;
+ self->relTol = 1e-5;
+ self->tol = 0.0;
+ self->rnorm = 0.0;
self->useInitial = False;
-
-#ifdef HAVE_PETSC
- self->outerSolver = (MatrixSolver*)PETScMatrixSolver_New( "" );
- PETScMatrixSolver_SetKSPType( self->outerSolver, PETScMatrixSolver_KSPType_Richardson );
- PETScMatrixSolver_SetPCType( self->outerSolver, PETScMatrixSolver_PCType_None );
- PETScMatrixSolver_SetUseInitialSolution( self->outerSolver, True );
- PETScMatrixSolver_SetMaxIterations( self->outerSolver, 1 );
-#else
self->outerSolver = NULL;
-#endif
}
@@ -198,7 +192,9 @@
assert( self && Stg_CheckType( self, MultigridSolver ) );
- MatrixSolver_SetRelativeTolerance( self->outerSolver, tolerance );
+ self->relTol = tolerance;
+ self->solversChanged = True;
+
}
void MultigridSolver_SetAbsoluteTolerance( void* matrixSolver, double tolerance ) {
@@ -235,8 +231,6 @@
MultigridSolver_LevelCycle( self, self->nLevels - 1, rhs, solution );
self->curIt++;
}
-
- self->matricesReady = False;
}
void MultigridSolver_Setup( void* matrixSolver, void* rhs, void* solution ) {
@@ -247,6 +241,7 @@
_MatrixSolver_Setup( self, rhs, solution );
+ /* Need to rebuild the operators? */
if( self->opGen )
rebuildOps = MGOpGenerator_HasExpired( self->opGen );
else
@@ -264,30 +259,31 @@
if( rebuildOps )
MultigridSolver_UpdateOps( self );
- if( !self->matricesReady ) {
+ if( self->solversChanged )
+ MultigridSolver_UpdateSolvers( self );
+ if( self->matrixChanged || rebuildOps || self->opsChanged || self->solversChanged ) {
MultigridSolver_UpdateMatrices( self );
MultigridSolver_UpdateWorkVectors( self );
- MatrixSolver_SetMatrix( self->outerSolver, self->matrix );
- self->matricesReady = True;
}
- if( !self->solversReady ) {
- MultigridSolver_UpdateSolvers( self );
- self->solversReady = True;
- }
+ self->solversChanged = False;
+ self->matrixChanged = False;
+ self->opsChanged = False;
}
MatrixSolver_Status MultigridSolver_GetSolveStatus( void* matrixSolver ) {
MultigridSolver* self = (MultigridSolver*)matrixSolver;
MatrixSolver_Status outerStatus;
+ Vector* residual;
assert( self && Stg_CheckType( self, MultigridSolver ) );
MatrixSolver_Solve( self->outerSolver, self->curRHS, self->curSolution );
outerStatus = MatrixSolver_GetSolveStatus( self->outerSolver );
- if( outerStatus == MatrixSolver_Status_ConvergedIterations )
- return MatrixSolver_Status_ConvergedIterations;
- else if( outerStatus == MatrixSolver_Status_DivergedIterations )
+
+ if( outerStatus == MatrixSolver_Status_ConvergedIterations || outerStatus == MatrixSolver_Status_Iterating )
return MatrixSolver_Status_Iterating;
+ else if( outerStatus == MatrixSolver_Status_ConvergedRelative )
+ return MatrixSolver_Status_ConvergedRelative;
else
return outerStatus;
}
@@ -329,106 +325,81 @@
assert( self && Stg_CheckType( self, MultigridSolver ) );
MultigridSolver_DestructLevels( self );
+
self->nLevels = nLevels;
self->levels = AllocArray( MultigridSolver_Level, nLevels );
MPI_Comm_size( self->comm, (int*)&nProcs );
+
for( l_i = 0; l_i < nLevels; l_i++ ) {
MultigridSolver_Level* level = self->levels + l_i;
+ level->downSolver = NULL;
level->nDownIts = 1;
+ level->upSolver = NULL;
+ level->nUpIts = l_i ? 1 : 0;
level->nCycles = 1;
- if( l_i == 0 ) {
- level->nUpIts = 0;
-
- /* Create direct solve. */
-#ifdef HAVE_PETSC
- level->downSolver = (MatrixSolver*)PETScMatrixSolver_New( "" );
- PETScMatrixSolver_SetKSPType( level->downSolver,
- PETScMatrixSolver_KSPType_PreOnly );
- if( nProcs == 1 ) {
- PETScMatrixSolver_SetPCType( level->downSolver,
- PETScMatrixSolver_PCType_LU );
- }
- else {
- PETScMatrixSolver_SetPCType( level->downSolver,
- PETScMatrixSolver_PCType_RedundantLU );
- }
-#else
- self->downSolver = NULL;
-#endif
- }
- else {
- level->nUpIts = 1;
-
- /* Create smoother. */
-#ifdef HAVE_PETSC
- level->downSolver = (MatrixSolver*)PETScMatrixSolver_New( "" );
- PETScMatrixSolver_SetKSPType( level->downSolver,
- PETScMatrixSolver_KSPType_Richardson );
- PETScMatrixSolver_SetPCType( level->downSolver,
- PETScMatrixSolver_PCType_SOR );
-#else
- level->downSolver = NULL;
-#endif
- }
- level->upSolver = level->downSolver;
- Stg_Class_AddRef( level->upSolver );
-
level->R = NULL;
level->P = NULL;
+
+ level->workRHS = NULL;
+ level->workSol = NULL;
}
}
-void MultigridSolver_SetRestriction( void* matrixSolver, unsigned level, void* _R ) {
+void MultigridSolver_SetRestriction( void* matrixSolver, unsigned levelInd, void* _R ) {
MultigridSolver* self = (MultigridSolver*)matrixSolver;
+ MultigridSolver_Level* level;
Matrix* R = (Matrix*)_R;
assert( self && Stg_CheckType( self, MultigridSolver ) );
- assert( level < self->nLevels && level > 0 );
+ assert( levelInd < self->nLevels && levelInd > 0 );
assert( !R || Stg_CheckType( R, Matrix ) );
- if( self->levels[level].R != R )
- self->matricesReady = 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( level->R )
+ Stg_Class_RemoveRef( level->R );
+ level->R = R;
if( R )
Stg_Class_AddRef( R );
}
-void MultigridSolver_SetProlongation( void* matrixSolver, unsigned level, void* _P ) {
+void MultigridSolver_SetProlongation( void* matrixSolver, unsigned levelInd, void* P ) {
MultigridSolver* self = (MultigridSolver*)matrixSolver;
- Matrix* P = (Matrix*)_P;
+ MultigridSolver_Level* level;
assert( self && Stg_CheckType( self, MultigridSolver ) );
- assert( level < self->nLevels && level > 0 );
+ assert( levelInd < self->nLevels && levelInd > 0 );
assert( !P || Stg_CheckType( P, Matrix ) );
- if( self->levels[level].P != P )
- self->matricesReady = 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( level->P )
+ Stg_Class_RemoveRef( level->P );
+ level->P = P;
if( P )
Stg_Class_AddRef( P );
}
-void MultigridSolver_SetLevelDownSolver( void* matrixSolver, unsigned level, MatrixSolver* solver ) {
+void MultigridSolver_SetLevelDownSolver( void* matrixSolver, unsigned levelInd, void* solver ) {
MultigridSolver* self = (MultigridSolver*)matrixSolver;
- MatrixSolver* curSolver;
+ MultigridSolver_Level* level;
assert( self && Stg_CheckType( self, MultigridSolver ) );
- assert( level < self->nLevels );
+ assert( levelInd < self->nLevels );
+ assert( !solver || Stg_CheckType( solver, MatrixSolver ) );
- curSolver = self->levels[level].downSolver;
- if( curSolver )
- Stg_Class_RemoveRef( curSolver );
- self->levels[level].downSolver = solver;
+ level = self->levels + levelInd;
+ if( level->downSolver != solver )
+ self->solversChanged = True;
+ if( level->downSolver )
+ Stg_Class_RemoveRef( level->downSolver );
+ level->downSolver = solver;
if( solver )
Stg_Class_AddRef( solver );
-
- self->solversReady = False;
}
void MultigridSolver_SetLevelDownIterations( void* matrixSolver, unsigned level, unsigned nIts ) {
@@ -438,25 +409,24 @@
assert( level < self->nLevels );
self->levels[level].nDownIts = nIts;
- self->solversReady = False;
+ self->solversChanged = True;
}
-void MultigridSolver_SetLevelUpSolver( void* matrixSolver, unsigned level, MatrixSolver* solver ) {
+void MultigridSolver_SetLevelUpSolver( void* matrixSolver, unsigned levelInd, void* solver ) {
MultigridSolver* self = (MultigridSolver*)matrixSolver;
- MatrixSolver* curSolver;
+ MultigridSolver_Level* level;
assert( self && Stg_CheckType( self, MultigridSolver ) );
- assert( level < self->nLevels );
- assert( level > 0 );
+ assert( levelInd < self->nLevels && levelInd > 0 );
- curSolver = self->levels[level].upSolver;
- if( curSolver )
- Stg_Class_RemoveRef( curSolver );
- self->levels[level].upSolver = solver;
+ level = self->levels + levelInd;
+ if( level->downSolver != solver )
+ self->solversChanged = True;
+ if( level->downSolver )
+ Stg_Class_RemoveRef( level->downSolver );
+ level->downSolver = solver;
if( solver )
Stg_Class_AddRef( solver );
-
- self->solversReady = False;
}
void MultigridSolver_SetLevelUpIterations( void* matrixSolver, unsigned level, unsigned nIts ) {
@@ -467,7 +437,7 @@
assert( level > 0 );
self->levels[level].nUpIts = nIts;
- self->solversReady = False;
+ self->solversChanged = True;
}
void MultigridSolver_SetLevelCycles( void* matrixSolver, unsigned level, unsigned nCycles ) {
@@ -478,10 +448,10 @@
assert( level > 0 );
self->levels[level].nCycles = nCycles;
- self->solversReady = False;
+ self->solversChanged = True;
}
-void MultigridSolver_SetAllDownSolver( void* matrixSolver, MatrixSolver* solver ) {
+void MultigridSolver_SetAllDownSolver( void* matrixSolver, void* solver ) {
MultigridSolver* self = (MultigridSolver*)matrixSolver;
unsigned l_i;
@@ -501,7 +471,7 @@
MultigridSolver_SetLevelDownIterations( self, l_i, nIts );
}
-void MultigridSolver_SetAllUpSolver( void* matrixSolver, MatrixSolver* solver ) {
+void MultigridSolver_SetAllUpSolver( void* matrixSolver, void* solver ) {
MultigridSolver* self = (MultigridSolver*)matrixSolver;
unsigned l_i;
@@ -521,7 +491,7 @@
MultigridSolver_SetLevelUpIterations( self, l_i, nIts );
}
-void MultigridSolver_SetAllSolver( void* matrixSolver, MatrixSolver* solver ) {
+void MultigridSolver_SetAllSolver( void* matrixSolver, void* solver ) {
MultigridSolver* self = (MultigridSolver*)matrixSolver;
assert( self );
@@ -530,7 +500,7 @@
MultigridSolver_SetAllUpSolver( self, solver );
}
-void MultigridSolver_SetCoarseSolver( void* matrixSolver, MatrixSolver* solver ) {
+void MultigridSolver_SetCoarseSolver( void* matrixSolver, void* solver ) {
MultigridSolver* self = (MultigridSolver*)matrixSolver;
assert( self );
@@ -538,34 +508,6 @@
MultigridSolver_SetLevelDownSolver( self, 0, solver );
}
-void MultigridSolver_SetWorkSolution( void* matrixSolver, unsigned level, void* vector ) {
- MultigridSolver* self = (MultigridSolver*)matrixSolver;
-
- assert( self && Stg_CheckType( self, MultigridSolver ) );
- assert( !vector || Stg_CheckType( vector, Vector ) );
- assert( level < self->nLevels - 1 );
-
- if( self->levels[level].workSol )
- Stg_Class_RemoveRef( self->levels[level].workSol );
- self->levels[level].workSol = vector;
- if( vector )
- Stg_Class_AddRef( vector );
-}
-
-void MultigridSolver_SetWorkRHS( void* matrixSolver, unsigned level, void* vector ) {
- MultigridSolver* self = (MultigridSolver*)matrixSolver;
-
- assert( self && Stg_CheckType( self, MultigridSolver ) );
- assert( !vector || Stg_CheckType( vector, Vector ) );
- assert( level < self->nLevels - 1 );
-
- if( self->levels[level].workRHS )
- Stg_Class_RemoveRef( self->levels[level].workRHS );
- self->levels[level].workRHS = vector;
- if( vector )
- Stg_Class_AddRef( vector );
-}
-
unsigned MultigridSolver_GetNumLevels( void* matrixSolver ) {
MultigridSolver* self = (MultigridSolver*)matrixSolver;
@@ -592,30 +534,12 @@
return self->levels[level].P;
}
-Vector* MultigridSolver_GetWorkSolution( void* matrixSolver, unsigned level ) {
- MultigridSolver* self = (MultigridSolver*)matrixSolver;
- assert( self && Stg_CheckType( self, MultigridSolver ) );
- assert( level < self->nLevels - 1 );
-
- return self->levels[level].workSol;
-}
-
-Vector* MultigridSolver_GetWorkRHS( void* matrixSolver, unsigned level ) {
- MultigridSolver* self = (MultigridSolver*)matrixSolver;
-
- assert( self && Stg_CheckType( self, MultigridSolver ) );
- assert( level < self->nLevels - 1);
-
- return self->levels[level].workRHS;
-}
-
-
/*----------------------------------------------------------------------------------------------------------------------------------
** Private Functions
*/
-void MultigridSolver_RestrictMatrix( MultigridSolver* self, MultigridSolver_Level* level, Matrix* dstMatrix ) {
+void MultigridSolver_RestrictMatrix( MultigridSolver* self, MultigridSolver_Level* level, Matrix** dstMatrix ) {
Matrix *srcMat;
unsigned nSrcRows, nSrcCols;
unsigned nRRows, nRCols;
@@ -628,9 +552,9 @@
Matrix_GetGlobalSize( srcMat, &nSrcRows, &nSrcCols );
Matrix_GetGlobalSize( level->R, &nRRows, &nRCols );
if( nRCols == nSrcRows )
- Matrix_PAPt( srcMat, level->R, dstMatrix );
+ Matrix_PAPt( srcMat, level->R, (void**)dstMatrix );
else
- Matrix_PtAP( srcMat, level->R, dstMatrix );
+ Matrix_PtAP( srcMat, level->R, (void**)dstMatrix );
}
void MultigridSolver_LevelCycle( MultigridSolver* self, unsigned levelInd, Vector* rhs, Vector* solution ) {
@@ -641,11 +565,8 @@
level = self->levels + levelInd;
- if( level->nDownIts ) {
- if( levelInd == self->nLevels - 1 && self->curIt > 0 )
- MatrixSolver_SetUseInitialSolution( level->downSolver, True );
+ if( level->nDownIts )
MatrixSolver_Solve( level->downSolver, rhs, solution );
- }
if( levelInd > 0 ) {
MultigridSolver_Level* nextLevel;
@@ -661,42 +582,10 @@
Matrix_MultiplyAdd( level->P, nextLevel->workSol, solution, solution );
}
- if( level->nUpIts ) {
- if( level->upSolver == level->downSolver )
- MatrixSolver_SetUseInitialSolution( level->upSolver, True );
+ if( level->nUpIts )
MatrixSolver_Solve( level->upSolver, rhs, solution );
- if( level->upSolver == level->downSolver )
- MatrixSolver_SetUseInitialSolution( level->upSolver, False );
- }
}
-void MultigridSolver_UpdateWorkVectors( MultigridSolver* self ) {
- MultigridSolver_Level* level;
- unsigned rowSize, colSize;
- Vector* vec;
- unsigned l_i;
-
- assert( self && Stg_CheckType( self, MultigridSolver ) );
-
- for( l_i = 0; l_i < self->nLevels - 1; l_i++ ) {
- level = self->levels + l_i;
-
- Matrix_GetLocalSize( MatrixSolver_GetMatrix( level->downSolver ), &rowSize, &colSize );
- vec = MultigridSolver_GetWorkSolution( self, l_i );
- if( !vec || Vector_GetLocalSize( vec ) != rowSize ) {
- Vector_Duplicate( self->curSolution, (void**)&vec );
- Vector_SetLocalSize( vec, rowSize );
- MultigridSolver_SetWorkSolution( self, l_i, vec );
- }
- vec = MultigridSolver_GetWorkRHS( self, l_i );
- if( !vec || Vector_GetLocalSize( vec ) != rowSize ) {
- Vector_Duplicate( self->curSolution, (void**)&vec );
- Vector_SetLocalSize( vec, rowSize );
- MultigridSolver_SetWorkRHS( self, l_i, vec );
- }
- }
-}
-
void MultigridSolver_UpdateSolvers( MultigridSolver* self ) {
MultigridSolver_Level* level;
unsigned l_i;
@@ -706,17 +595,40 @@
for( l_i = self->nLevels - 1; l_i < self->nLevels; l_i-- ) {
level = self->levels + l_i;
- if( MatrixSolver_GetMaxIterations( level->downSolver ) != level->nDownIts )
- MatrixSolver_SetMaxIterations( level->downSolver, level->nDownIts );
- if( level->upSolver != level->downSolver ) {
+ if( l_i > 0 ) {
+ if( !level->downSolver )
+ level->downSolver = MultigridSolver_CreateSmoother( self );
+ if( !level->upSolver )
+ level->upSolver = MultigridSolver_CreateSmoother( self );
+
+ if( MatrixSolver_GetMaxIterations( level->downSolver ) != level->nDownIts )
+ MatrixSolver_SetMaxIterations( level->downSolver, level->nDownIts );
if( MatrixSolver_GetMaxIterations( level->upSolver ) != level->nUpIts )
MatrixSolver_SetMaxIterations( level->upSolver, level->nUpIts );
- }
- MatrixSolver_SetUseInitialSolution( level->downSolver, False );
- if( level->upSolver != level->downSolver )
+ if( l_i == self->nLevels - 1 )
+ MatrixSolver_SetUseInitialSolution( level->downSolver, True );
+ else
+ MatrixSolver_SetUseInitialSolution( level->downSolver, False );
MatrixSolver_SetUseInitialSolution( level->upSolver, True );
+ }
+ else {
+ if( !level->downSolver )
+ level->downSolver = MultigridSolver_CreateCoarseSolver( self );
+ }
}
+
+#ifdef HAVE_PETSC
+ self->outerSolver = (MatrixSolver*)PETScMatrixSolver_New( "" );
+ PETScMatrixSolver_SetKSPType( self->outerSolver, PETScMatrixSolver_KSPType_Richardson );
+ PETScMatrixSolver_SetPCType( self->outerSolver, PETScMatrixSolver_PCType_SOR );
+ PETScMatrixSolver_SetMaxIterations( self->outerSolver, 1 );
+ PETScMatrixSolver_SetUseInitialSolution( self->outerSolver, True );
+ PETScMatrixSolver_SetNormType( self->outerSolver, PETScMatrixSolver_NormType_Preconditioned );
+ PETScMatrixSolver_SetRelativeTolerance( self->outerSolver, self->relTol );
+#else
+#error Oi! Where is PETSc?
+#endif
}
void MultigridSolver_UpdateMatrices( MultigridSolver* self ) {
@@ -732,40 +644,123 @@
Matrix* mat;
mat = MatrixSolver_GetMatrix( level->downSolver );
- if( !mat ) {
- Matrix_Duplicate( self->matrix, (void**)&mat );
- Matrix_AssemblyBegin( mat );
- Matrix_AssemblyEnd( mat );
+ if( mat ) {
+ KillObject( mat );
}
- MultigridSolver_RestrictMatrix( self, self->levels + l_i + 1, mat );
-
+ MultigridSolver_RestrictMatrix( self, self->levels + l_i + 1, &mat );
MatrixSolver_SetMatrix( level->downSolver, mat );
- if( level->upSolver != level->downSolver )
+ if( l_i > 0 )
MatrixSolver_SetMatrix( level->upSolver, mat );
}
else {
MatrixSolver_SetMatrix( level->downSolver, self->matrix );
- if( level->upSolver != level->downSolver )
+ if( l_i > 0 )
MatrixSolver_SetMatrix( level->upSolver, self->matrix );
}
}
+
+ MatrixSolver_SetMatrix( self->outerSolver, self->matrix );
}
void MultigridSolver_UpdateOps( MultigridSolver* self ) {
- Matrix **pOps, **rOps;
- unsigned l_i;
+ MultigridSolver_Level* level;
+ Matrix **pOps, **rOps;
+ unsigned l_i;
assert( self && Stg_CheckType( self, MultigridSolver ) );
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] );
+ level = self->levels + l_i;
+
+ if( !level->P ) {
+ level->P = pOps[l_i];
+ Stg_Class_AddRef( level->P );
+ }
+ else
+ Stg_Class_RemoveRef( pOps[l_i] );
+
+ if( !level->R ) {
+ level->R = rOps[l_i];
+ Stg_Class_AddRef( level->R );
+ }
+ else
+ Stg_Class_RemoveRef( rOps[l_i] );
}
+
FreeArray( pOps );
FreeArray( rOps );
}
+void MultigridSolver_UpdateWorkVectors( MultigridSolver* self ) {
+ MultigridSolver_Level* level;
+ unsigned rowSize, colSize;
+ unsigned l_i;
+
+ assert( self && Stg_CheckType( self, MultigridSolver ) );
+
+ for( l_i = 0; l_i < self->nLevels - 1; l_i++ ) {
+ level = self->levels + l_i;
+
+ Matrix_GetLocalSize( MatrixSolver_GetMatrix( level->downSolver ), &rowSize, &colSize );
+
+ if( !level->workSol || Vector_GetLocalSize( level->workSol ) != rowSize ) {
+ if( level->workSol )
+ Stg_Class_RemoveRef( level->workSol );
+ Vector_Duplicate( self->curSolution, (void**)&level->workSol );
+ Vector_SetLocalSize( level->workSol, rowSize );
+ }
+
+ if( !level->workRHS || Vector_GetLocalSize( level->workRHS ) != rowSize ) {
+ if( level->workRHS )
+ Stg_Class_RemoveRef( level->workRHS );
+ Vector_Duplicate( self->curSolution, (void**)&level->workRHS );
+ Vector_SetLocalSize( level->workRHS, rowSize );
+ }
+ }
+}
+
+MatrixSolver* MultigridSolver_CreateSmoother( MultigridSolver* self ) {
+ MatrixSolver* smoother;
+
+#ifdef HAVE_PETSC
+ smoother = (MatrixSolver*)PETScMatrixSolver_New( "" );
+ PETScMatrixSolver_SetKSPType( smoother,
+ PETScMatrixSolver_KSPType_Richardson );
+ PETScMatrixSolver_SetPCType( smoother,
+ PETScMatrixSolver_PCType_SOR );
+#else
+#error *** Oi! Where is PETSc?
+#endif
+
+ return smoother;
+}
+
+MatrixSolver* MultigridSolver_CreateCoarseSolver( MultigridSolver* self ) {
+ MatrixSolver* coarseSolver;
+ unsigned nProcs;
+
+ MPI_Comm_size( self->comm, (int*)&nProcs );
+
+#ifdef HAVE_PETSC
+ coarseSolver = (MatrixSolver*)PETScMatrixSolver_New( "" );
+ PETScMatrixSolver_SetKSPType( coarseSolver,
+ PETScMatrixSolver_KSPType_PreOnly );
+ if( nProcs == 1 ) {
+ PETScMatrixSolver_SetPCType( coarseSolver,
+ PETScMatrixSolver_PCType_LU );
+ }
+ else {
+ PETScMatrixSolver_SetPCType( coarseSolver,
+ PETScMatrixSolver_PCType_RedundantLU );
+ }
+#else
+#error *** Oi! Where is PETSc?
+#endif
+
+ return coarseSolver;
+}
+
void MultigridSolver_Destruct( MultigridSolver* self ) {
assert( self && Stg_CheckType( self, MultigridSolver ) );
@@ -781,18 +776,29 @@
for( l_i = 0; l_i < self->nLevels; l_i++ ) {
MultigridSolver_Level* level = self->levels + l_i;
- if( level->downSolver )
+ if( level->downSolver ) {
+ Stg_Class_RemoveRef( MatrixSolver_GetMatrix( level->downSolver ) );
Stg_Class_RemoveRef( level->downSolver );
- if( level->upSolver )
+ }
+ if( level->upSolver ) {
+ Stg_Class_RemoveRef( MatrixSolver_GetMatrix( level->upSolver ) );
Stg_Class_RemoveRef( level->upSolver );
+ }
if( level->R )
Stg_Class_RemoveRef( level->R );
if( level->P )
Stg_Class_RemoveRef( level->P );
+ if( level->workRHS )
+ Stg_Class_RemoveRef( level->workRHS );
+ if( level->workSol )
+ Stg_Class_RemoveRef( level->workSol );
}
KillArray( self->levels );
self->nLevels = 0;
- self->matricesReady = False;
- self->solversReady = False;
+ self->solversChanged = True;
+ self->opsChanged = True;
+
+ /* Temporary. */
+ KillObject( self->outerSolver );
}
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MultigridSolver.h
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MultigridSolver.h 2007-02-23 18:03:47 UTC (rev 6092)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MultigridSolver.h 2007-02-23 18:03:49 UTC (rev 6093)
@@ -71,11 +71,14 @@
unsigned nLevels; \
MultigridSolver_Level* levels; \
MGOpGenerator* opGen; \
- Bool matricesReady; \
- Bool solversReady; \
+ Bool solversChanged; \
+ Bool opsChanged; \
\
unsigned curIt; \
unsigned maxIts; \
+ double relTol; \
+ double tol; \
+ double rnorm; \
Bool useInitial; \
MatrixSolver* outerSolver;
@@ -127,35 +130,36 @@
void MultigridSolver_SetLevels( void* matrixSolver, unsigned nLevels );
void MultigridSolver_SetRestriction( void* matrixSolver, unsigned level, void* R );
void MultigridSolver_SetProlongation( void* matrixSolver, unsigned level, void* P );
- void MultigridSolver_SetLevelDownSolver( void* matrixSolver, unsigned level, MatrixSolver* solver );
+ void MultigridSolver_SetLevelDownSolver( void* matrixSolver, unsigned level, void* solver );
void MultigridSolver_SetLevelDownIterations( void* matrixSolver, unsigned level, unsigned nIts );
- void MultigridSolver_SetLevelUpSolver( void* matrixSolver, unsigned level, MatrixSolver* solver );
+ void MultigridSolver_SetLevelUpSolver( void* matrixSolver, unsigned level, void* solver );
void MultigridSolver_SetLevelUpIterations( void* matrixSolver, unsigned level, unsigned nIts );
void MultigridSolver_SetLevelCycles( void* matrixSolver, unsigned level, unsigned nCycles );
- void MultigridSolver_SetAllDownSolver( void* matrixSolver, MatrixSolver* solver );
+ void MultigridSolver_SetAllDownSolver( void* matrixSolver, void* solver );
void MultigridSolver_SetAllDownIterations( void* matrixSolver, unsigned level, unsigned nIts );
- void MultigridSolver_SetAllUpSolver( void* matrixSolver, MatrixSolver* solver );
+ void MultigridSolver_SetAllUpSolver( void* matrixSolver, void* solver );
void MultigridSolver_SetAllUpIterations( void* matrixSolver, unsigned level, unsigned nIts );
- void MultigridSolver_SetAllSolver( void* matrixSolver, MatrixSolver* solver );
- void MultigridSolver_SetCoarseSolver( void* matrixSolver, MatrixSolver* solver );
+ void MultigridSolver_SetAllSolver( void* matrixSolver, void* solver );
+ void MultigridSolver_SetCoarseSolver( void* matrixSolver, void* solver );
unsigned MultigridSolver_GetNumLevels( void* matrixSolver );
Matrix* MultigridSolver_GetRestriction( void* matrixSolver, unsigned level );
Matrix* MultigridSolver_GetProlongation( void* matrixSolver, unsigned level );
- Vector* MultigridSolver_GetWorkSolution( void* matrixSolver, unsigned level );
- Vector* MultigridSolver_GetWorkRHS( void* matrixSolver, unsigned level );
/*--------------------------------------------------------------------------------------------------------------------------
** Private Member functions
*/
- void MultigridSolver_RestrictMatrix( MultigridSolver* self, MultigridSolver_Level* level, Matrix* dstMatrix );
+ void MultigridSolver_RestrictMatrix( MultigridSolver* self, MultigridSolver_Level* level, Matrix** dstMatrix );
void MultigridSolver_LevelCycle( MultigridSolver* self, unsigned levelInd, Vector* rhs, Vector* solution );
void MultigridSolver_UpdateWorkVectors( MultigridSolver* self );
void MultigridSolver_UpdateSolvers( MultigridSolver* self );
void MultigridSolver_UpdateMatrices( MultigridSolver* self );
void MultigridSolver_UpdateOps( MultigridSolver* self );
+ MatrixSolver* MultigridSolver_CreateSmoother( MultigridSolver* self );
+ MatrixSolver* MultigridSolver_CreateCoarseSolver( MultigridSolver* self );
+
void MultigridSolver_Destruct( MultigridSolver* self );
void MultigridSolver_DestructLevels( MultigridSolver* self );
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:47 UTC (rev 6092)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMGSolver.c 2007-02-23 18:03:49 UTC (rev 6093)
@@ -97,6 +97,7 @@
self->levels = NULL;
self->opGen = NULL;
self->solversChanged = True;
+ self->opsChanged = True;
}
@@ -204,17 +205,18 @@
}
}
- if( self->solversChanged || rebuildOps )
+ if( self->solversChanged || rebuildOps || self->opsChanged )
PETScMGSolver_UpdateSolvers( self );
if( rebuildOps )
PETScMGSolver_UpdateOps( self );
- if( self->matrixChanged || self->solversChanged || rebuildOps ) {
+ if( self->matrixChanged || self->solversChanged || rebuildOps || self->opsChanged ) {
PETScMGSolver_UpdateMatrices( self );
PETScMGSolver_UpdateWorkVectors( self );
}
self->solversChanged = False;
self->matrixChanged = False;
+ self->opsChanged = False;
}
@@ -403,14 +405,14 @@
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 );
+ }
- ec = PCMGGetSmootherDown( pc, l_i, &levelKSP );
- CheckPETScError( ec );
- ec = KSPGetOperators( levelKSP, &downMat, &pMat, &flag );
- CheckPETScError( ec );
- ec = KSPSetOperators( levelKSP, mats[l_i], mats[l_i], DIFFERENT_NONZERO_PATTERN );
- CheckPETScError( ec );
- }
+ ec = PCMGGetSmootherDown( pc, l_i, &levelKSP );
+ CheckPETScError( ec );
+ ec = KSPGetOperators( levelKSP, &downMat, &pMat, &flag );
+ CheckPETScError( ec );
+ ec = KSPSetOperators( levelKSP, mats[l_i], mats[l_i], DIFFERENT_NONZERO_PATTERN );
+ CheckPETScError( ec );
ec = PCMGSetResidual( pc, l_i, PCMGDefaultResidual, mats[l_i] );
CheckPETScError( ec );
@@ -446,6 +448,9 @@
assert( self && Stg_CheckType( self, PETScMGSolver ) );
+ if( self->nLevels == 1 )
+ return;
+
ec = KSPGetPC( self->ksp, &pc );
CheckPETScError( ec );
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:47 UTC (rev 6092)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMGSolver.h 2007-02-23 18:03:49 UTC (rev 6093)
@@ -70,7 +70,8 @@
unsigned nLevels; \
PETScMGSolver_Level* levels; \
MGOpGenerator* opGen; \
- Bool solversChanged;
+ Bool solversChanged; \
+ Bool opsChanged;
struct PETScMGSolver { __PETScMGSolver };
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:47 UTC (rev 6092)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrix.c 2007-02-23 18:03:49 UTC (rev 6093)
@@ -424,7 +424,8 @@
CheckPETScError( ec );
}
-void PETScMatrix_PAPt( void* matrix, void* _P, void* _dstMatrix ) {
+void PETScMatrix_PAPt( void* matrix, void* _P, void** _dstMatrix ) {
+#if 0
PETScMatrix* self = (PETScMatrix*)matrix;
PETScMatrix* P = (PETScMatrix*)_P;
PETScMatrix* dstMatrix = (PETScMatrix*)_dstMatrix;
@@ -454,12 +455,13 @@
CheckPETScError( ec );
dstMatrix->hasChanged = True;
+#endif
}
-void PETScMatrix_PtAP( void* matrix, void* _P, void* _dstMatrix ) {
+void PETScMatrix_PtAP( void* matrix, void* _P, void** _dstMatrix ) {
PETScMatrix* self = (PETScMatrix*)matrix;
PETScMatrix* P = (PETScMatrix*)_P;
- PETScMatrix* dstMatrix = (PETScMatrix*)_dstMatrix;
+ PETScMatrix** dstMatrix = (PETScMatrix**)_dstMatrix;
PetscErrorCode ec;
#if 0
double fillRatio;
@@ -470,7 +472,7 @@
assert( self && Stg_CheckType( self, PETScMatrix ) );
assert( P && Stg_CheckType( P, PETScMatrix ) );
- assert( dstMatrix && Stg_CheckType( dstMatrix, PETScMatrix ) );
+ assert( dstMatrix && (!(*dstMatrix) || Stg_CheckType( *dstMatrix, PETScMatrix )) );
#if 0
MatGetInfo( self->petscMat, MAT_GLOBAL_SUM, &mInfo );
@@ -485,13 +487,17 @@
fillRatio = nzC / (nzA + nzP);
#endif
- if( dstMatrix->petscMat != PETSC_NULL )
- ec = MatDestroy( dstMatrix->petscMat );
+ if( *dstMatrix )
+ ec = MatPtAP( self->petscMat, P->petscMat, MAT_REUSE_MATRIX, 1.0, &(*dstMatrix)->petscMat );
+ else {
+ Matrix_Duplicate( self, dstMatrix );
+ if( (*dstMatrix)->petscMat )
+ MatDestroy( (*dstMatrix)->petscMat );
+ ec = MatPtAP( self->petscMat, P->petscMat, MAT_INITIAL_MATRIX, 1.0, &(*dstMatrix)->petscMat );
+ }
CheckPETScError( ec );
- ec = MatPtAP( self->petscMat, P->petscMat, MAT_INITIAL_MATRIX, 1.0, &dstMatrix->petscMat );
- CheckPETScError( ec );
- dstMatrix->hasChanged = True;
+ (*dstMatrix)->hasChanged = True;
}
void PETScMatrix_MatrixMultiply( void* matrix, void* _matrix0, void* _dstMatrix ) {
@@ -654,7 +660,7 @@
assert( dstMatrix );
*dstMatrix = self->_defaultConstructor( "" );
- ec = MatDuplicate( self->petscMat, MAT_DO_NOT_COPY_VALUES, &(*dstMatrix)->petscMat );
+ ec = MatCreate( self->comm, &(*dstMatrix)->petscMat );
CheckPETScError( ec );
(*dstMatrix)->hasChanged = True;
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrix.h
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrix.h 2007-02-23 18:03:47 UTC (rev 6092)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrix.h 2007-02-23 18:03:49 UTC (rev 6093)
@@ -106,8 +106,8 @@
void PETScMatrix_Multiply( void* matrix, void* vector, void* dstVector );
void PETScMatrix_TransposeMultiply( void* matrix, void* vector, void* dstVector );
void PETScMatrix_MultiplyAdd( void* matrix, void* vector0, void* vector1, void* dstVector );
- void PETScMatrix_PAPt( void* matrix, void* P, void* dstMatrix );
- void PETScMatrix_PtAP( void* matrix, void* P, void* dstMatrix );
+ void PETScMatrix_PAPt( void* matrix, void* P, void** dstMatrix );
+ void PETScMatrix_PtAP( void* matrix, void* P, void** dstMatrix );
void PETScMatrix_MatrixMultiply( void* matrix, void* matrix0, void* dstMatrix );
double PETScMatrix_L2Norm( void* matrix );
void PETScMatrix_Transpose( void* matrix, void* dstMatrix );
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrixSolver.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrixSolver.c 2007-02-23 18:03:47 UTC (rev 6092)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrixSolver.c 2007-02-23 18:03:49 UTC (rev 6093)
@@ -245,14 +245,36 @@
MatrixSolver_Status PETScMatrixSolver_GetSolveStatus( void* matrixSolver ) {
PETScMatrixSolver* self = (PETScMatrixSolver*)matrixSolver;
+ PC pc;
+ KSPType kspType;
+ PCType pcType;
KSPConvergedReason reason;
PetscErrorCode ec;
assert( self && Stg_CheckType( self, PETScMatrixSolver ) );
- ec = KSPGetConvergedReason( self->ksp, &reason );
+ ec = KSPGetType( self->ksp, &kspType );
CheckPETScError( ec );
+ ec = KSPGetPC( self->ksp, &pc );
+ CheckPETScError( ec );
+ ec = PCGetType( pc, &pcType );
+ CheckPETScError( ec );
+ if( !strcmp( kspType, KSPRICHARDSON ) && !strcmp( pcType, PCSOR ) ) {
+ double rnorm;
+ unsigned curIt;
+
+ rnorm = PETScMatrixSolver_GetResidualNorm( self );
+ curIt = PETScMatrixSolver_GetIterations( self );
+ PETScMatrixSolver_SetNormType( self, PETScMatrixSolver_NormType_Preconditioned );
+ ec = KSPDefaultConverged( self->ksp, (PetscInt)curIt, (PetscScalar)rnorm, &reason, PETSC_NULL );
+ CheckPETScError( ec );
+ }
+ else {
+ ec = KSPGetConvergedReason( self->ksp, &reason );
+ CheckPETScError( ec );
+ }
+
return reason;
}
@@ -284,15 +306,33 @@
double PETScMatrixSolver_GetResidualNorm( void* matrixSolver ) {
PETScMatrixSolver* self = (PETScMatrixSolver*)matrixSolver;
- PetscScalar norm;
+ PC pc;
+ KSPType kspType;
+ PCType pcType;
+ PetscScalar rnorm;
PetscErrorCode ec;
assert( self && Stg_CheckType( self, PETScMatrixSolver ) );
- ec = KSPGetResidualNorm( self->ksp, &norm );
+ ec = KSPGetType( self->ksp, &kspType );
CheckPETScError( ec );
+ ec = KSPGetPC( self->ksp, &pc );
+ CheckPETScError( ec );
+ ec = PCGetType( pc, &pcType );
+ CheckPETScError( ec );
- return (double)norm;
+ if( !strcmp( kspType, KSPRICHARDSON ) && !strcmp( pcType, PCSOR ) ) {
+ Vector* residual;
+
+ residual = MatrixSolver_GetResidual( self );
+ rnorm = (PetscScalar)Vector_L2Norm( residual );
+ }
+ else {
+ ec = KSPGetResidualNorm( self->ksp, &rnorm );
+ CheckPETScError( ec );
+ }
+
+ return (double)rnorm;
}
@@ -387,7 +427,17 @@
CheckPETScError( ec );
}
+void PETScMatrixSolver_SetNormType( void* matrixSolver, PETScMatrixSolver_NormType normType ) {
+ PETScMatrixSolver* self = (PETScMatrixSolver*)matrixSolver;
+ PetscErrorCode ec;
+ assert( self && Stg_CheckType( self, PETScMatrixSolver ) );
+
+ ec = KSPSetNormType( self->ksp, (KSPNormType)normType );
+ 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-02-23 18:03:47 UTC (rev 6092)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrixSolver.h 2007-02-23 18:03:49 UTC (rev 6093)
@@ -106,6 +106,7 @@
void PETScMatrixSolver_SetKSPType( void* matrixSolver, PETScMatrixSolver_KSPType type );
void PETScMatrixSolver_SetPCType( void* matrixSolver, PETScMatrixSolver_PCType type );
void PETScMatrixSolver_EnableShifting( void* matrixSolver, Bool state );
+ void PETScMatrixSolver_SetNormType( void* matrixSolver, PETScMatrixSolver_NormType normType );
/*--------------------------------------------------------------------------------------------------------------------------
** 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-02-23 18:03:47 UTC (rev 6092)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/SROpGenerator.c 2007-02-23 18:03:49 UTC (rev 6093)
@@ -223,6 +223,7 @@
CartesianGenerator *fGen, *cGen;
unsigned nDims;
unsigned* cSize;
+ double crdMin[3], crdMax[3];
unsigned d_i;
assert( self && Stg_CheckType( self, SROpGenerator ) );
@@ -248,7 +249,8 @@
for( d_i = 0; d_i < nDims; d_i++ )
cSize[d_i] = fGen->elGrid->sizes[d_i] / 2;
CartesianGenerator_SetTopologyParams( cGen, cSize, fGen->maxDecompDims, fGen->minDecomp, fGen->maxDecomp );
- CartesianGenerator_SetGeometryParams( cGen, fGen->crdMin, fGen->crdMax );
+ Mesh_GetGlobalCoordRange( fMesh, crdMin, crdMax );
+ CartesianGenerator_SetGeometryParams( cGen, crdMin, crdMax );
CartesianGenerator_SetShadowDepth( cGen, 0 );
FreeArray( cSize );
@@ -438,7 +440,9 @@
SROpGenerator_GenLevelOp( self, l_i, P );
pOps[l_i] = P;
+ Stg_Class_AddRef( P );
rOps[l_i] = P;
+ Stg_Class_AddRef( P );
}
}
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/Vector.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/Vector.c 2007-02-23 18:03:47 UTC (rev 6092)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/Vector.c 2007-02-23 18:03:49 UTC (rev 6093)
@@ -168,7 +168,7 @@
Vector_GetArray( self, &array );
Journal_Printf( stream, "%s = [\n", self->name );
for( entry_i = 0; entry_i < size; entry_i++ )
- Journal_Printf( stream , "\t%u: \t %2.4g\n", entry_i, array[entry_i] );
+ Journal_Printf( stream , "\t%u: \t %.12g\n", entry_i, array[entry_i] );
Journal_Printf( stream, "];\n");
Vector_RestoreArray( self, &array );
}
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/types.h
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/types.h 2007-02-23 18:03:47 UTC (rev 6092)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/types.h 2007-02-23 18:03:49 UTC (rev 6093)
@@ -94,5 +94,14 @@
PETScMatrixSolver_PCType_None
} PETScMatrixSolver_PCType;
+ /* Norm enumerations. */
+ typedef enum {
+ PETScMatrixSolver_NormType_None = 0,
+ PETScMatrixSolver_NormType_Preconditioned = 1,
+ PETScMatrixSolver_NormType_Unpreconditioned = 2,
+ PETScMatrixSolver_NormType_Natural = 3
+ } PETScMatrixSolver_NormType;
+
+
#endif /* __StgFEM_SLE_LinearAlgrebra_types_h__ */
More information about the cig-commits
mailing list