[cig-commits] r6091 - in long/3D/Gale/trunk/src/StgFEM: .
SLE/LinearAlgebra/src
walter at geodynamics.org
walter at geodynamics.org
Fri Feb 23 10:03:43 PST 2007
Author: walter
Date: 2007-02-23 10:03:43 -0800 (Fri, 23 Feb 2007)
New Revision: 6091
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/MatrixSolver.h
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
Log:
r977 at earth (orig r717): LukeHodkinson | 2007-01-21 22:17:10 -0800
Updating 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:716
db209038-57f2-0310-97fa-b160e0ae9d04:/trunk:669
+ 38867592-cf10-0410-9e16-a142ea72ac34:/cig:880
db209038-57f2-0310-97fa-b160e0ae9d04:/branches/decomp3d:717
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:02:37 UTC (rev 6090)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/Matrix.c 2007-02-23 18:03:43 UTC (rev 6091)
@@ -98,6 +98,9 @@
assert( self && Stg_CheckType( self, Matrix ) );
self->comm = MPI_COMM_WORLD;
+ self->hasChanged = True;
+ self->solvers = List_New();
+ List_SetItemSize( self->solvers, sizeof(MatrixSolver*) );
}
@@ -155,12 +158,36 @@
self->comm = comm;
}
+void _Matrix_AssemblyBegin( void* matrix ) {
+ Matrix* self = (Matrix*)matrix;
+ assert( self && Stg_CheckType( self, Matrix ) );
+
+ if( self->hasChanged ) {
+ Matrix_InvalidateSolvers( self );
+ self->hasChanged = False;
+ }
+}
+
+
/*--------------------------------------------------------------------------------------------------------------------------
** Public Functions
*/
+void Matrix_InvalidateSolvers( void* matrix ) {
+ Matrix* self = (Matrix*)matrix;
+ MatrixSolver* solver;
+ unsigned s_i;
+ assert( self && Stg_CheckType( self, Matrix ) );
+
+ for( s_i = 0; s_i < List_GetSize( self->solvers ); s_i++ ) {
+ solver = *List_Get( self->solvers, s_i, MatrixSolver* );
+ solver->matrixChanged = True;
+ }
+}
+
+
/*----------------------------------------------------------------------------------------------------------------------------------
** Private Functions
*/
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:02:37 UTC (rev 6090)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/Matrix.h 2007-02-23 18:03:43 UTC (rev 6091)
@@ -124,7 +124,9 @@
Matrix_DumpFunc* dumpFunc; \
\
/* Matrix info */ \
- MPI_Comm comm;
+ MPI_Comm comm; \
+ Bool hasChanged; \
+ List* solvers;
struct Matrix { __Matrix };
@@ -218,6 +220,7 @@
void _Matrix_Destroy( void* matrix, void* data );
void _Matrix_SetComm( void* matrix, MPI_Comm comm );
+ void _Matrix_AssemblyBegin( void* matrix );
#define Matrix_SetComm( self, comm ) \
VirtualCall( self, setCommFunc, self, comm )
@@ -313,6 +316,8 @@
** Public functions
*/
+ void Matrix_InvalidateSolvers( void* matrix );
+
/*--------------------------------------------------------------------------------------------------------------------------
** Private Member functions
*/
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:02:37 UTC (rev 6090)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MatrixSolver.c 2007-02-23 18:03:43 UTC (rev 6091)
@@ -84,6 +84,7 @@
self->inversion = NULL;
self->residual = NULL;
self->expiredResidual = True;
+ self->matrixChanged = True;
self->curRHS = NULL;
self->curSolution = NULL;
@@ -144,13 +145,22 @@
self->comm = comm;
}
-void _MatrixSolver_SetMatrix( void* matrixSolver, void* matrix ) {
+void _MatrixSolver_SetMatrix( void* matrixSolver, void* _matrix ) {
MatrixSolver* self = (MatrixSolver*)matrixSolver;
+ Matrix* matrix = (Matrix*)_matrix;
assert( self && Stg_CheckType( self, MatrixSolver ) );
assert( !matrix || Stg_CheckType( matrix, Matrix ) );
+ if( matrix == self->matrix )
+ return;
+
+ if( self->matrix )
+ List_Remove( self->matrix->solvers, &self );
self->matrix = matrix;
+ if( matrix )
+ List_Append( matrix->solvers, &self );
+ self->matrixChanged = True;
}
void _MatrixSolver_Setup( void* matrixSolver, void* rhs, void* solution ) {
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MatrixSolver.h
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MatrixSolver.h 2007-02-23 18:02:37 UTC (rev 6090)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MatrixSolver.h 2007-02-23 18:03:43 UTC (rev 6091)
@@ -87,6 +87,7 @@
Matrix* inversion; \
Vector* residual; \
Bool expiredResidual; \
+ Bool matrixChanged; \
\
Vector* curRHS; \
Vector* curSolution;
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:02:37 UTC (rev 6090)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MultigridSolver.c 2007-02-23 18:03:43 UTC (rev 6091)
@@ -96,7 +96,8 @@
self->nLevels = 0;
self->levels = NULL;
self->opGen = NULL;
- self->ready = False;
+ self->matricesReady = False;
+ self->solversReady = False;
self->curIt = 0;
self->maxIts = 500;
@@ -234,19 +235,15 @@
MultigridSolver_LevelCycle( self, self->nLevels - 1, rhs, solution );
self->curIt++;
}
+
+ self->matricesReady = False;
}
-void MultigridSolver_Setup( void* matrixSolver, void* _rhs, void* _solution ) {
+void MultigridSolver_Setup( void* matrixSolver, void* rhs, void* solution ) {
MultigridSolver* self = (MultigridSolver*)matrixSolver;
- Vector* rhs = (Vector*)_rhs;
- Vector* solution = (Vector*)_solution;
Bool rebuildOps;
- Matrix **pOps, **rOps;
- unsigned l_i;
assert( self && Stg_CheckType( self, MultigridSolver ) );
- assert( rhs && Stg_CheckType( rhs, Vector ) );
- assert( solution && Stg_CheckType( solution, Vector ) );
_MatrixSolver_Setup( self, rhs, solution );
@@ -254,68 +251,29 @@
rebuildOps = MGOpGenerator_HasExpired( self->opGen );
else
rebuildOps = False;
- if( self->ready && !rebuildOps )
- return;
+ if( !rebuildOps ) {
+ unsigned l_i;
- 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;
-
- if( MatrixSolver_GetMaxIterations( level->downSolver ) != level->nDownIts )
- MatrixSolver_SetMaxIterations( level->downSolver, level->nDownIts );
- if( level->upSolver != level->downSolver ) {
- if( MatrixSolver_GetMaxIterations( level->upSolver ) != level->nUpIts )
- MatrixSolver_SetMaxIterations( level->upSolver, level->nUpIts );
- }
-
- if( l_i < self->nLevels - 1 ) {
- Matrix* mat;
-
- mat = MatrixSolver_GetMatrix( level->downSolver );
- if( !mat ) {
- Matrix_Duplicate( self->matrix, (void**)&mat );
- Matrix_AssemblyBegin( mat );
- Matrix_AssemblyEnd( mat );
+ for( l_i = 1; l_i < self->nLevels; l_i++ ) {
+ if( !MultigridSolver_GetRestriction( self, l_i ) || !MultigridSolver_GetProlongation( self, l_i ) ) {
+ rebuildOps = True;
+ break;
}
- MultigridSolver_RestrictMatrix( self, self->levels + l_i + 1, mat );
-
- MatrixSolver_SetMatrix( level->downSolver, mat );
- if( level->upSolver != level->downSolver )
- MatrixSolver_SetMatrix( level->upSolver, mat );
}
- else {
- MatrixSolver_SetMatrix( level->downSolver, self->matrix );
- if( level->upSolver != level->downSolver )
- MatrixSolver_SetMatrix( level->upSolver, self->matrix );
- }
-
- MatrixSolver_SetUseInitialSolution( level->downSolver, False );
- if( level->upSolver != level->downSolver )
- MatrixSolver_SetUseInitialSolution( level->upSolver, True );
}
- for( l_i = 1; l_i < self->nLevels; l_i++ ) {
- MultigridSolver_Level* level = self->levels + l_i;
- MultigridSolver_Level* nextLevel = self->levels + l_i - 1;
- unsigned rowSize, colSize;
-
- Matrix_GetLocalSize( MatrixSolver_GetMatrix( nextLevel->downSolver ), &rowSize, &colSize );
- Vector_Duplicate( solution, (void**)&level->nextSol );
- Vector_SetLocalSize( level->nextSol, rowSize );
- Vector_Duplicate( rhs, (void**)&level->nextRHS );
- Vector_SetLocalSize( level->nextRHS, rowSize );
+ if( rebuildOps )
+ MultigridSolver_UpdateOps( self );
+ if( !self->matricesReady ) {
+ MultigridSolver_UpdateMatrices( self );
+ MultigridSolver_UpdateWorkVectors( self );
+ MatrixSolver_SetMatrix( self->outerSolver, self->matrix );
+ self->matricesReady = True;
}
-
- MatrixSolver_SetMatrix( self->outerSolver, self->matrix );
-
- self->ready = True;
+ if( !self->solversReady ) {
+ MultigridSolver_UpdateSolvers( self );
+ self->solversReady = True;
+ }
}
MatrixSolver_Status MultigridSolver_GetSolveStatus( void* matrixSolver ) {
@@ -377,7 +335,7 @@
for( l_i = 0; l_i < nLevels; l_i++ ) {
MultigridSolver_Level* level = self->levels + l_i;
- level->nDownIts = 5;
+ level->nDownIts = 1;
level->nCycles = 1;
if( l_i == 0 ) {
@@ -396,13 +354,12 @@
PETScMatrixSolver_SetPCType( level->downSolver,
PETScMatrixSolver_PCType_RedundantLU );
}
- PETScMatrixSolver_EnableShifting( level->downSolver, True );
#else
self->downSolver = NULL;
#endif
}
else {
- level->nUpIts = 5;
+ level->nUpIts = 1;
/* Create smoother. */
#ifdef HAVE_PETSC
@@ -432,7 +389,7 @@
assert( !R || Stg_CheckType( R, Matrix ) );
if( self->levels[level].R != R )
- self->ready = False;
+ self->matricesReady = False;
if( self->levels[level].R )
Stg_Class_RemoveRef( self->levels[level].R );
self->levels[level].R = R;
@@ -449,7 +406,7 @@
assert( !P || Stg_CheckType( P, Matrix ) );
if( self->levels[level].P != P )
- self->ready = False;
+ self->matricesReady = False;
if( self->levels[level].P )
Stg_Class_RemoveRef( self->levels[level].P );
self->levels[level].P = P;
@@ -471,7 +428,7 @@
if( solver )
Stg_Class_AddRef( solver );
- self->ready = False;
+ self->solversReady = False;
}
void MultigridSolver_SetLevelDownIterations( void* matrixSolver, unsigned level, unsigned nIts ) {
@@ -481,7 +438,7 @@
assert( level < self->nLevels );
self->levels[level].nDownIts = nIts;
- self->ready = False;
+ self->solversReady = False;
}
void MultigridSolver_SetLevelUpSolver( void* matrixSolver, unsigned level, MatrixSolver* solver ) {
@@ -499,7 +456,7 @@
if( solver )
Stg_Class_AddRef( solver );
- self->ready = False;
+ self->solversReady = False;
}
void MultigridSolver_SetLevelUpIterations( void* matrixSolver, unsigned level, unsigned nIts ) {
@@ -510,7 +467,7 @@
assert( level > 0 );
self->levels[level].nUpIts = nIts;
- self->ready = False;
+ self->solversReady = False;
}
void MultigridSolver_SetLevelCycles( void* matrixSolver, unsigned level, unsigned nCycles ) {
@@ -521,7 +478,7 @@
assert( level > 0 );
self->levels[level].nCycles = nCycles;
- self->ready = False;
+ self->solversReady = False;
}
void MultigridSolver_SetAllDownSolver( void* matrixSolver, MatrixSolver* solver ) {
@@ -581,6 +538,34 @@
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;
@@ -589,7 +574,43 @@
return self->nLevels;
}
+Matrix* MultigridSolver_GetRestriction( void* matrixSolver, unsigned level ) {
+ MultigridSolver* self = (MultigridSolver*)matrixSolver;
+ assert( self && Stg_CheckType( self, MultigridSolver ) );
+ assert( level < self->nLevels && level > 0 );
+
+ return self->levels[level].R;
+}
+
+Matrix* MultigridSolver_GetProlongation( void* matrixSolver, unsigned level ) {
+ MultigridSolver* self = (MultigridSolver*)matrixSolver;
+
+ assert( self && Stg_CheckType( self, MultigridSolver ) );
+ assert( level < self->nLevels && level > 0 );
+
+ 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
*/
@@ -627,15 +648,17 @@
}
if( levelInd > 0 ) {
- unsigned c_i;
+ MultigridSolver_Level* nextLevel;
+ unsigned c_i;
+ nextLevel = self->levels - (levelInd - 1);
if( level->R == level->P )
- Matrix_TransposeMultiply( level->R, MatrixSolver_GetResidual( level->downSolver ), level->nextRHS );
+ Matrix_TransposeMultiply( level->R, MatrixSolver_GetResidual( level->downSolver ), nextLevel->workRHS );
else
- Matrix_Multiply( level->R, MatrixSolver_GetResidual( level->downSolver ), level->nextRHS );
+ Matrix_Multiply( level->R, MatrixSolver_GetResidual( level->downSolver ), nextLevel->workRHS );
for( c_i = 0; c_i < level->nCycles; c_i++ )
- MultigridSolver_LevelCycle( self, levelInd - 1, level->nextRHS, level->nextSol );
- Matrix_MultiplyAdd( level->P, level->nextSol, solution, solution );
+ MultigridSolver_LevelCycle( self, levelInd - 1, nextLevel->workRHS, nextLevel->workSol );
+ Matrix_MultiplyAdd( level->P, nextLevel->workSol, solution, solution );
}
if( level->nUpIts ) {
@@ -647,6 +670,102 @@
}
}
+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;
+
+ assert( self && Stg_CheckType( self, MultigridSolver ) );
+
+ 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( MatrixSolver_GetMaxIterations( level->upSolver ) != level->nUpIts )
+ MatrixSolver_SetMaxIterations( level->upSolver, level->nUpIts );
+ }
+
+ MatrixSolver_SetUseInitialSolution( level->downSolver, False );
+ if( level->upSolver != level->downSolver )
+ MatrixSolver_SetUseInitialSolution( level->upSolver, True );
+ }
+}
+
+void MultigridSolver_UpdateMatrices( MultigridSolver* self ) {
+ MultigridSolver_Level* level;
+ unsigned l_i;
+
+ assert( self && Stg_CheckType( self, MultigridSolver ) );
+
+ for( l_i = self->nLevels - 1; l_i < self->nLevels; l_i-- ) {
+ level = self->levels + l_i;
+
+ if( l_i < self->nLevels - 1 ) {
+ Matrix* mat;
+
+ mat = MatrixSolver_GetMatrix( level->downSolver );
+ if( !mat ) {
+ Matrix_Duplicate( self->matrix, (void**)&mat );
+ Matrix_AssemblyBegin( mat );
+ Matrix_AssemblyEnd( mat );
+ }
+ MultigridSolver_RestrictMatrix( self, self->levels + l_i + 1, mat );
+
+ MatrixSolver_SetMatrix( level->downSolver, mat );
+ if( level->upSolver != level->downSolver )
+ MatrixSolver_SetMatrix( level->upSolver, mat );
+ }
+ else {
+ MatrixSolver_SetMatrix( level->downSolver, self->matrix );
+ if( level->upSolver != level->downSolver )
+ MatrixSolver_SetMatrix( level->upSolver, self->matrix );
+ }
+ }
+}
+
+void MultigridSolver_UpdateOps( MultigridSolver* self ) {
+ 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] );
+ }
+ FreeArray( pOps );
+ FreeArray( rOps );
+}
+
void MultigridSolver_Destruct( MultigridSolver* self ) {
assert( self && Stg_CheckType( self, MultigridSolver ) );
@@ -674,5 +793,6 @@
KillArray( self->levels );
self->nLevels = 0;
- self->ready = False;
+ self->matricesReady = False;
+ self->solversReady = False;
}
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:02:37 UTC (rev 6090)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MultigridSolver.h 2007-02-23 18:03:43 UTC (rev 6091)
@@ -57,9 +57,8 @@
Matrix* R;
Matrix* P;
- Vector* nextRHS;
- Vector* nextSol;
- Vector* residual;
+ Vector* workRHS;
+ Vector* workSol;
} MultigridSolver_Level;
#define __MultigridSolver \
@@ -72,7 +71,8 @@
unsigned nLevels; \
MultigridSolver_Level* levels; \
MGOpGenerator* opGen; \
- Bool ready; \
+ Bool matricesReady; \
+ Bool solversReady; \
\
unsigned curIt; \
unsigned maxIts; \
@@ -140,6 +140,10 @@
void MultigridSolver_SetCoarseSolver( void* matrixSolver, MatrixSolver* 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
@@ -147,7 +151,10 @@
void MultigridSolver_RestrictMatrix( MultigridSolver* self, MultigridSolver_Level* level, Matrix* dstMatrix );
void MultigridSolver_LevelCycle( MultigridSolver* self, unsigned levelInd, Vector* rhs, Vector* solution );
- void MultigridSolver_CalcResidual( MultigridSolver* self, Vector* solution, Vector* residual );
+ void MultigridSolver_UpdateWorkVectors( MultigridSolver* self );
+ void MultigridSolver_UpdateSolvers( MultigridSolver* self );
+ void MultigridSolver_UpdateMatrices( MultigridSolver* self );
+ void MultigridSolver_UpdateOps( 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:02:37 UTC (rev 6090)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMGSolver.c 2007-02-23 18:03:43 UTC (rev 6091)
@@ -49,30 +49,30 @@
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_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 ) {
@@ -96,7 +96,7 @@
self->nLevels = 0;
self->levels = NULL;
self->opGen = NULL;
- self->ready = False;
+ self->solversChanged = True;
}
@@ -181,19 +181,11 @@
PETScMatrixSolver_SetPCType( self, PETScMatrixSolver_PCType_Multigrid );
}
-void PETScMGSolver_Setup( void* matrixSolver, void* rhs, void* _solution ) {
+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 );
@@ -201,57 +193,28 @@
rebuildOps = MGOpGenerator_HasExpired( self->opGen );
else
rebuildOps = False;
- if( self->ready && !rebuildOps )
- return;
+ if( !rebuildOps ) {
+ unsigned l_i;
- 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 = 1; l_i < self->nLevels; l_i++ ) {
- PCMGGetSmoother( pc, l_i, &downKSP );
- KSPSetType( downKSP, KSPRICHARDSON );
- KSPGetPC( downKSP, &downPC );
- PCSetType( downPC, PCSOR );
+ for( l_i = 1; l_i < self->nLevels; l_i++ ) {
+ if( !self->levels[l_i].R || !self->levels[l_i].P ) {
+ rebuildOps = True;
+ break;
+ }
+ }
}
- 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 );
+ if( self->solversChanged || rebuildOps )
+ PETScMGSolver_UpdateSolvers( self );
+ if( rebuildOps )
+ PETScMGSolver_UpdateOps( self );
+ if( self->matrixChanged || self->solversChanged || rebuildOps ) {
+ PETScMGSolver_UpdateMatrices( self );
+ PETScMGSolver_UpdateWorkVectors( self );
}
- FreeArray( pOps );
- FreeArray( rOps );
- self->ready = True;
+ self->solversChanged = False;
+ self->matrixChanged = False;
}
@@ -272,16 +235,19 @@
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->nDownIts = 2;
+ level->nUpIts = (l_i == 0) ? 0 : 2;
level->nCycles = 1;
level->R = NULL;
level->P = NULL;
level->workRes = PETSC_NULL;
+ level->workSol = PETSC_NULL;
+ level->workRHS = PETSC_NULL;
}
}
void PETScMGSolver_SetRestriction( void* matrixSolver, unsigned level, void* _R ) {
+#if 0
PETScMGSolver* self = (PETScMGSolver*)matrixSolver;
Matrix* R = (Matrix*)_R;
@@ -296,9 +262,11 @@
self->levels[level].R = R;
if( R )
Stg_Class_AddRef( R );
+#endif
}
void PETScMGSolver_SetProlongation( void* matrixSolver, unsigned level, void* _P ) {
+#if 0
PETScMGSolver* self = (PETScMGSolver*)matrixSolver;
Matrix* P = (Matrix*)_P;
@@ -313,6 +281,7 @@
self->levels[level].P = P;
if( P )
Stg_Class_AddRef( P );
+#endif
}
void PETScMGSolver_SetLevelDownIterations( void* matrixSolver, unsigned level, unsigned nIts ) {
@@ -322,7 +291,7 @@
assert( level < self->nLevels );
self->levels[level].nDownIts = nIts;
- self->ready = False;
+ self->solversChanged = False;
}
void PETScMGSolver_SetLevelUpIterations( void* matrixSolver, unsigned level, unsigned nIts ) {
@@ -333,7 +302,7 @@
assert( level > 0 );
self->levels[level].nUpIts = nIts;
- self->ready = False;
+ self->solversChanged = False;
}
void PETScMGSolver_SetLevelCycles( void* matrixSolver, unsigned level, unsigned nCycles ) {
@@ -344,7 +313,7 @@
assert( level > 0 );
self->levels[level].nCycles = nCycles;
- self->ready = False;
+ self->solversChanged = False;
}
void PETScMGSolver_SetAllDownIterations( void* matrixSolver, unsigned level, unsigned nIts ) {
@@ -380,6 +349,223 @@
** Private Functions
*/
+void PETScMGSolver_UpdateOps( PETScMGSolver* self ) {
+ PC pc;
+ PETScMatrix **pOps, **rOps;
+ PetscErrorCode ec;
+ unsigned l_i;
+
+ assert( self && Stg_CheckType( self, PETScMGSolver ) );
+
+ ec = KSPGetPC( self->ksp, &pc );
+ 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;
+ ec = PCMGSetInterpolate( pc, l_i, pOps[l_i]->petscMat );
+ CheckPETScError( ec );
+ self->levels[l_i].R = rOps[l_i]->petscMat;
+ ec = PCMGSetRestriction( pc, l_i, rOps[l_i]->petscMat );
+ CheckPETScError( ec );
+ }
+
+ FreeArray( pOps );
+ FreeArray( rOps );
+}
+
+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;
+
+ assert( self && Stg_CheckType( self, PETScMGSolver ) );
+ assert( self->matrix && Stg_CheckType( self->matrix, PETScMatrix ) );
+
+ 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 );
+
+ 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 );
+
+ if( l_i > 0 ) {
+ PCMGGetSmootherUp( pc, l_i, &levelKSP );
+ ec = KSPGetOperators( levelKSP, &upMat, &pMat, &flag );
+ 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;
+ PetscErrorCode ec;
+ PetscInt nr, nc;
+ unsigned l_i;
+
+ assert( self && Stg_CheckType( self, PETScMGSolver ) );
+
+ ec = KSPGetPC( self->ksp, &pc );
+ CheckPETScError( ec );
+
+ 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 );
+ 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 );
+ 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 );
+ 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 );
+ }
+ }
+ }
+}
+
+void PETScMGSolver_UpdateSolvers( PETScMGSolver* self ) {
+ PETScMGSolver_Level* level;
+ PC pc;
+ KSP levelKSP;
+ PC levelPC;
+ PetscErrorCode ec;
+ unsigned l_i;
+
+ assert( self && Stg_CheckType( self, PETScMGSolver ) );
+
+ ec = KSPGetPC( self->ksp, &pc );
+ CheckPETScError( ec );
+
+ ec = PCMGSetLevels( pc, self->nLevels, PETSC_NULL );
+ CheckPETScError( ec );
+ ec = PCMGSetType( pc, PC_MG_MULTIPLICATIVE );
+ CheckPETScError( ec );
+
+ for( l_i = 1; l_i < self->nLevels; l_i++ ) {
+ level = self->levels + l_i;
+
+ ec = PCMGGetSmootherDown( pc, l_i, &levelKSP );
+ CheckPETScError( ec );
+ ec = KSPSetType( levelKSP, KSPRICHARDSON );
+ CheckPETScError( ec );
+ ec = KSPGetPC( levelKSP, &levelPC );
+ CheckPETScError( ec );
+ ec = PCSetType( levelPC, PCSOR );
+ CheckPETScError( ec );
+ ec = KSPSetTolerances( levelKSP, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, level->nDownIts );
+ CheckPETScError( ec );
+ if( l_i == self->nLevels - 1 )
+ ec = KSPSetInitialGuessNonzero( levelKSP, PETSC_TRUE );
+ else
+ ec = KSPSetInitialGuessNonzero( levelKSP, PETSC_FALSE );
+ CheckPETScError( ec );
+
+ ec = PCMGGetSmootherUp( pc, l_i, &levelKSP );
+ CheckPETScError( ec );
+ ec = KSPSetType( levelKSP, KSPRICHARDSON );
+ CheckPETScError( ec );
+ ec = KSPGetPC( levelKSP, &levelPC );
+ CheckPETScError( ec );
+ ec = PCSetType( levelPC, PCSOR );
+ CheckPETScError( ec );
+ ec = KSPSetTolerances( levelKSP, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, level->nUpIts );
+ CheckPETScError( ec );
+ ec = KSPSetInitialGuessNonzero( levelKSP, PETSC_TRUE );
+ CheckPETScError( ec );
+
+ ec = PCMGSetCyclesOnLevel( pc, l_i, level->nCycles );
+ CheckPETScError( ec );
+ }
+}
+
void PETScMGSolver_Destruct( PETScMGSolver* self ) {
assert( self && Stg_CheckType( self, PETScMGSolver ) );
@@ -403,5 +589,5 @@
KillArray( self->levels );
self->nLevels = 0;
- self->ready = False;
+ self->solversChanged = 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:02:37 UTC (rev 6090)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMGSolver.h 2007-02-23 18:03:43 UTC (rev 6091)
@@ -52,10 +52,12 @@
unsigned nUpIts;
unsigned nCycles;
- Matrix* R;
- Matrix* P;
+ Mat R;
+ Mat P;
Vec workRes;
+ Vec workSol;
+ Vec workRHS;
} PETScMGSolver_Level;
#define __PETScMGSolver \
@@ -68,7 +70,7 @@
unsigned nLevels; \
PETScMGSolver_Level* levels; \
MGOpGenerator* opGen; \
- Bool ready;
+ Bool solversChanged;
struct PETScMGSolver { __PETScMGSolver };
@@ -126,6 +128,11 @@
** Private Member functions
*/
+ void PETScMGSolver_UpdateOps( PETScMGSolver* self );
+ void PETScMGSolver_UpdateMatrices( PETScMGSolver* self );
+ void PETScMGSolver_UpdateWorkVectors( PETScMGSolver* self );
+ void PETScMGSolver_UpdateSolvers( PETScMGSolver* self );
+
void PETScMGSolver_Destruct( PETScMGSolver* self );
void PETScMGSolver_DestructLevels( PETScMGSolver* self );
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:02:37 UTC (rev 6090)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrix.c 2007-02-23 18:03:43 UTC (rev 6091)
@@ -183,6 +183,8 @@
MatDestroy( self->petscMat );
ec = MatCreate( self->comm, &self->petscMat );
CheckPETScError( ec );
+
+ self->hasChanged = True;
}
void PETScMatrix_SetGlobalSize( void* matrix, unsigned nRows, unsigned nColumns ) {
@@ -195,6 +197,8 @@
CheckPETScError( ec );
ec = MatSetFromOptions( self->petscMat );
CheckPETScError( ec );
+
+ self->hasChanged = True;
}
void PETScMatrix_SetLocalSize( void* matrix, unsigned nRows, unsigned nColumns ) {
@@ -207,6 +211,8 @@
CheckPETScError( ec );
ec = MatSetFromOptions( self->petscMat );
CheckPETScError( ec );
+
+ self->hasChanged = True;
}
void PETScMatrix_AddEntries( void* matrix, unsigned nRows, unsigned* rowIndices,
@@ -226,6 +232,8 @@
(PetscInt)nColumns, (PetscInt*)columnIndices,
(PetscScalar*)values, ADD_VALUES );
CheckPETScError( ec );
+
+ self->hasChanged = True;
}
void PETScMatrix_InsertEntries( void* matrix, unsigned nRows, unsigned* rowIndices,
@@ -244,6 +252,8 @@
(PetscInt)nColumns, (PetscInt*)columnIndices,
(PetscScalar*)values, INSERT_VALUES );
CheckPETScError( ec );
+
+ self->hasChanged = True;
}
void PETScMatrix_DiagonalAddEntries( void* matrix, void* _vector ) {
@@ -256,6 +266,8 @@
ec = MatDiagonalSet( self->petscMat, vector->petscVec, ADD_VALUES );
CheckPETScError( ec );
+
+ self->hasChanged = True;
}
void PETScMatrix_DiagonalInsertEntries( void* matrix, void* _vector ) {
@@ -268,6 +280,8 @@
ec = MatDiagonalSet( self->petscMat, vector->petscVec, INSERT_VALUES );
CheckPETScError( ec );
+
+ self->hasChanged = True;
}
void PETScMatrix_Zero( void* matrix ) {
@@ -278,6 +292,8 @@
ec = MatZeroEntries( self->petscMat );
CheckPETScError( ec );
+
+ self->hasChanged = True;
}
void PETScMatrix_Load( void* matrix, char* filename ) {
@@ -295,6 +311,8 @@
CheckPETScError( ec );
ec = PetscViewerDestroy( viewer );
CheckPETScError( ec );
+
+ self->hasChanged = True;
}
void PETScMatrix_AssemblyBegin( void* matrix ) {
@@ -303,6 +321,8 @@
assert( self && Stg_CheckType( self, PETScMatrix ) );
+ _Matrix_AssemblyBegin( self );
+
ec = MatAssemblyBegin( self->petscMat, MAT_FINAL_ASSEMBLY );
CheckPETScError( ec );
}
@@ -325,6 +345,8 @@
ec = MatScale( self->petscMat, (PetscScalar)factor );
CheckPETScError( ec );
+
+ self->hasChanged = True;
}
void PETScMatrix_AddScaled( void* matrix, double factor, void* _matrix0 ) {
@@ -337,6 +359,8 @@
ec = MatAXPY( self->petscMat, (PetscScalar)factor, matrix0->petscMat, SAME_NONZERO_PATTERN );
CheckPETScError( ec );
+
+ self->hasChanged = True;
}
void PETScMatrix_DiagonalScale( void* matrix, void* _leftVector, void* _rightVector ) {
@@ -352,6 +376,8 @@
ec = MatDiagonalScale( self->petscMat, leftVector ? leftVector->petscVec : PETSC_NULL,
rightVector ? rightVector->petscVec : PETSC_NULL );
CheckPETScError( ec );
+
+ self->hasChanged = True;
}
void PETScMatrix_Multiply( void* matrix, void* _vector, void* _dstVector ) {
@@ -426,6 +452,8 @@
assert( 0 );
/*ec = MatPAPt( self->petscMat, P->petscMat, MAT_REUSE_MATRIX, (PetscScalar)fillRatio, &dstMatrix->petscMat );*/
CheckPETScError( ec );
+
+ dstMatrix->hasChanged = True;
}
void PETScMatrix_PtAP( void* matrix, void* _P, void* _dstMatrix ) {
@@ -462,6 +490,8 @@
CheckPETScError( ec );
ec = MatPtAP( self->petscMat, P->petscMat, MAT_INITIAL_MATRIX, 1.0, &dstMatrix->petscMat );
CheckPETScError( ec );
+
+ dstMatrix->hasChanged = True;
}
void PETScMatrix_MatrixMultiply( void* matrix, void* _matrix0, void* _dstMatrix ) {
@@ -498,6 +528,8 @@
CheckPETScError( ec );
ec = MatPtAP( self->petscMat, matrix0->petscMat, MAT_INITIAL_MATRIX, 1.0, &dstMatrix->petscMat );
CheckPETScError( ec );
+
+ dstMatrix->hasChanged = True;
}
double PETScMatrix_L2Norm( void* matrix ) {
@@ -530,6 +562,8 @@
else
ec = MatTranspose( self->petscMat, PETSC_NULL );
CheckPETScError( ec );
+
+ dstMatrix->hasChanged = True;
}
void PETScMatrix_GetGlobalSize( void* matrix, unsigned* nRows, unsigned* nColumns ) {
@@ -611,16 +645,19 @@
CheckPETScError( ec );
}
-void PETScMatrix_Duplicate( void* matrix, void** dstMatrix ) {
+void PETScMatrix_Duplicate( void* matrix, void** _dstMatrix ) {
PETScMatrix* self = (PETScMatrix*)matrix;
+ PETScMatrix** dstMatrix = (PETScMatrix**)_dstMatrix;
PetscErrorCode ec;
assert( self && Stg_CheckType( self, PETScMatrix ) );
assert( dstMatrix );
*dstMatrix = self->_defaultConstructor( "" );
- ec = MatDuplicate( self->petscMat, MAT_DO_NOT_COPY_VALUES, &((PETScMatrix*)*dstMatrix)->petscMat );
+ ec = MatDuplicate( self->petscMat, MAT_DO_NOT_COPY_VALUES, &(*dstMatrix)->petscMat );
CheckPETScError( ec );
+
+ (*dstMatrix)->hasChanged = True;
}
void PETScMatrix_CopyEntries( void* matrix, void* _dstMatrix ) {
@@ -633,6 +670,8 @@
ec = MatCopy( self->petscMat, dstMatrix->petscMat, DIFFERENT_NONZERO_PATTERN );
CheckPETScError( ec );
+
+ dstMatrix->hasChanged = True;
}
void PETScMatrix_Dump( void* matrix, char* filename ) {
@@ -683,6 +722,8 @@
ec = MatSeqAIJSetPreallocation( self->petscMat, nNonZeros, PETSC_NULL );
}
CheckPETScError( ec );
+
+ self->hasChanged = True;
}
void PETScMatrix_Draw( void* matrix ) {
More information about the cig-commits
mailing list