Skip to content

Commit

Permalink
Eigen matrix efficiency improved, postbifurcation using linesearch im…
Browse files Browse the repository at this point in the history 
…proved

Eigen matrix now setup using triplets to improve efficiency
  • Loading branch information
@msmejkal
msmejkal committed Sep 4, 2024
1 parent 3367c5a commit 10bc90c
Show file tree
Hide file tree
Showing 4 changed files with 164 additions and 208 deletions.
230 changes: 132 additions & 98 deletions src/oofemlib/eigenmtrx.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -68,154 +68,164 @@ int EigenMtrx ::buildInternalStructure( EngngModel *eModel, int di, const Unknow
int neq = eModel->giveNumberOfDomainEquations( di, s );

EigMat.resize( neq, neq ); // Resize the matrix

// allocation map
std ::vector<std ::set<int> > columns( neq );

for ( auto &elem : domain->giveElements() ) {
elem->giveLocationArray( loc, s );

for ( int ii : loc ) {
if ( ii > 0 ) {
for ( int jj : loc ) {
if ( jj > 0 ) {
columns[jj - 1].insert( ii - 1 ); // for each column nonzero rows are stored
}
}
}
}
}


// loop over active boundary conditions
std ::vector<IntArray> r_locs;
std ::vector<IntArray> c_locs;

for ( auto &gbc : domain->giveBcs() ) {
ActiveBoundaryCondition *bc = dynamic_cast<ActiveBoundaryCondition *>( gbc.get() );
if ( bc != NULL ) {
bc->giveLocationArrays( r_locs, c_locs, UnknownCharType, s, s );
for ( std ::size_t k = 0; k < r_locs.size(); k++ ) {
IntArray &krloc = r_locs[k];
IntArray &kcloc = c_locs[k];
for ( int ii : krloc ) {
if ( ii ) {
for ( int jj : kcloc ) {
if ( jj ) {
columns[jj - 1].insert( ii - 1 );
}
}
}
}
}
}
}

std::vector<int> ColReserve( neq ); // Allocate vector of reserved number of nonzeros for each column

for ( int i = 0; i < neq; i++ ) {
ColReserve[i] = columns[i].size();
}

EigMat.reserve( ColReserve ); // Reserve the number of nonzeros
triplets.reserve( 500*neq );


//// allocation map
//std ::vector<std ::set<int> > columns( neq );

//// Eigen::Triplet<double> T;


//for ( auto &elem : domain->giveElements() ) {
// elem->giveLocationArray( loc, s );

// for ( int ii : loc ) {
// if ( ii > 0 ) {
// for ( int jj : loc ) {
// if ( jj > 0 ) {
// columns[jj - 1].insert( ii - 1 ); // for each column nonzero rows are stored

// //Eigen::Triplet<double> T( ii - 1, jj - 1, 0. ); // Create triplet storing (row, col, val)
// //this->triplets.push_back( T ); // Add to vector of triplets
// }
// }
// }
// }
//}


//// loop over active boundary conditions
//std ::vector<IntArray> r_locs;
//std ::vector<IntArray> c_locs;

//for ( auto &gbc : domain->giveBcs() ) {
// ActiveBoundaryCondition *bc = dynamic_cast<ActiveBoundaryCondition *>( gbc.get() );
// if ( bc != NULL ) {
// bc->giveLocationArrays( r_locs, c_locs, UnknownCharType, s, s );
// for ( std ::size_t k = 0; k < r_locs.size(); k++ ) {
// IntArray &krloc = r_locs[k];
// IntArray &kcloc = c_locs[k];
// for ( int ii : krloc ) {
// if ( ii ) {
// for ( int jj : kcloc ) {
// if ( jj ) {
// columns[jj - 1].insert( ii - 1 );

// //Eigen::Triplet<double> T( ii - 0, jj - 0, 0. ); // Create triplet storing (row, col, val)
// //this->triplets.push_back( T ); // Add to vector of triplets
// }
// }
// }
// }
// }
// }
//}

//// std::vector<int> ColReserve( neq ); // Allocate vector of reserved number of nonzeros for each column
//Eigen::VectorXi ColReserve( neq ); // Allocate vector of reserved number of nonzeros for each column

//for ( int i = 0; i < neq; i++ ) {
// ColReserve[i] = columns[i].size() * 2; // overestimate two times
//}


//EigMat.reserve( ColReserve ); // Reserve the number of nonzeros
//this->EigMat.setFromTriplets( this->triplets.begin(), this->triplets.end() );
// this->EigMat.makeCompressed();

this->version++;

return true;
}


int EigenMtrx ::assemble( const IntArray &loc, const FloatMatrix &mat )
{
int dim = mat.giveNumberOfRows();

int ii, jj;
for ( int j = 0; j < dim; j++ ) {
int jj = loc[j];
jj = loc[j];
if ( jj ) {
for ( int i = 0; i < dim; i++ ) {
int ii = loc[i];
ii = loc[i];
if ( ii ) {
EigMat.coeffRef( ii - 1, jj - 1 ) += mat( i, j );
//EigMat.coeffRef( ii - 1, jj - 1 ) += mat( i, j );

// when set from triplets
this->triplets.push_back( Eigen::Triplet<double>( ii - 1, jj - 1, mat( i, j ) ) ); // Add to vector of triplets
}
}
}

}

//std::cout << "EigMat = " << std::endl
// << EigMat << std::endl;

this->version++;
return 1;

}

int EigenMtrx ::assemble( const IntArray &rloc, const IntArray &cloc, const FloatMatrix &mat )
{
int dim1, dim2;

dim1 = mat.giveNumberOfRows();
dim2 = mat.giveNumberOfColumns();

for ( int j = 0; j < dim2; j++ ) {
for ( int j = 0; j < mat.giveNumberOfColumns(); j++ ) {
int jj = cloc[j];
if ( jj ) {
for ( int i = 0; i < dim1; i++ ) {
for ( int i = 0; i < mat.giveNumberOfRows(); i++ ) {
int ii = rloc[i];
if ( ii ) {
EigMat.coeffRef( ii - 1, jj - 1 ) += mat( i, j );
//EigMat.coeffRef( ii - 1, jj - 1 ) += mat( i, j );

// when set from triplets
this->triplets.push_back( Eigen::Triplet<double>( ii - 1, jj - 1, mat( i, j ) ) ); // Add to vector of triplets
}
}
}

}

this->version++;

return 1;
}

void EigenMtrx ::zero()
{
EigMat.setZero();
this->triplets.clear(); // When set from triplets
this->version++;
}

double &EigenMtrx ::at( int i, int j )
{
this->version++;
if ( i > this->giveNumberOfRows() && j > this->giveNumberOfColumns() ) {
OOFEM_ERROR( "Array accessing exception -- (%d,%d) out of bounds", i, j );
} else {
double a = EigMat.coeff( i, j );
return a;
return this->giveMatrix().coeffRef( i, j );
}
}


double EigenMtrx ::at( int i, int j ) const
{
// Create copy
Eigen::SparseMatrix<double> EigMatCopy( this->EigMat );
if ( this->BuiltFromTripletsAtVersion != this->version ) {
EigMatCopy.setFromTriplets( this->triplets.begin(), this->triplets.end() );
}

if ( i > this->giveNumberOfRows() && j > this->giveNumberOfColumns() ) {

OOFEM_ERROR( "Array accessing exception -- (%d,%d) out of bounds", i, j );

} else {

double a = EigMat.coeff( i, j );

return a;
//return EigMat.coeff( i, j );
return EigMatCopy.coeff( i, j );
}
}

Eigen::SparseMatrix<double> EigenMtrx::giveMatrix()
Eigen::SparseMatrix<double>& EigenMtrx::giveMatrix()
{
this->createMatrixFromTriplets(); // When the matrix is created from triplets
return EigMat;
}



template <typename Derived>
Eigen::SparseSolverBase<Derived>& EigenMtrx::giveFactorization( FactorizationType Factorization )
Eigen::SparseSolverBase<Derived> &EigenMtrx::giveFactorization( FactorizationType Factorization )
{
if ( !this->isFactorized( Factorization ) ) {
this->computeFactorization( Factorization )
Expand All @@ -237,11 +247,10 @@ Eigen::SparseSolverBase<Derived>& EigenMtrx::giveFactorization( FactorizationTyp
default:
OOFEM_ERROR( "Unknown factorization type" );
}

}


SimplicialLDLTderived<Eigen::SparseMatrix<double>>& EigenMtrx::giveLDLTFactorization()
SimplicialLDLTderived<Eigen::SparseMatrix<double> > &EigenMtrx::giveLDLTFactorization()
{
if ( !this->isFactorized( FT_LDLT ) ) {
this->computeFactorization( FT_LDLT );
Expand All @@ -251,11 +260,11 @@ SimplicialLDLTderived<Eigen::SparseMatrix<double>>& EigenMtrx::giveLDLTFactoriza
}



bool EigenMtrx::isFactorized(FactorizationType factorizationType) {
bool EigenMtrx::isFactorized( FactorizationType factorizationType )
{
if ( this->versionUpdate != this->version ) {
areFactorized = { false, false, false, false };
}
}
return this->areFactorized[factorizationType];
}

Expand All @@ -264,38 +273,63 @@ void EigenMtrx::computeFactorization( FactorizationType factorizationType )
{
switch ( factorizationType ) {
case FT_LLT:
this->LLT_factorization.compute( EigMat );
this->LLT_factorization.compute( this->giveMatrix() );
break;
case FT_LU:
this->LU_factorization.compute( EigMat );
this->LU_factorization.compute( this->giveMatrix() );
break;
case FT_QR:
this->QR_factorization.compute( EigMat );
this->QR_factorization.compute( this->giveMatrix() );
break;
case FT_LDLT:
this->LDLT_factorization.compute( EigMat );
this->LDLT_factorization.compute( this->giveMatrix() );
break;
default:
OOFEM_ERROR( "Unknown factorization type" );
}
this->versionUpdate = this->version;
this->versionUpdate = this->version;
this->areFactorized[factorizationType] = true;
}


void EigenMtrx::printYourself() const{
// xreate dense matrix
auto denseMat = Eigen::MatrixXd( this->EigMat );
void EigenMtrx::printYourself() const
{
// Create copy
Eigen::SparseMatrix<double> EigMatCopy( this->EigMat );
if ( this->BuiltFromTripletsAtVersion != this->version ) {
EigMatCopy.setFromTriplets( this->triplets.begin(), this->triplets.end() );
}

// create dense matrix
auto denseMat = Eigen::MatrixXd( EigMatCopy );
std::cout << denseMat << std::endl;
}


void EigenMtrx::times( const FloatArray &x, FloatArray &answer ) const
{
FloatArray xcopy = x;
Eigen::VectorXd xeig = Eigen::Map<Eigen::VectorXd, Eigen::Unaligned>( xcopy.givePointer(), x.giveSize() );
Eigen::VectorXd answereig = this->EigMat * xeig;
answer = FloatArray( answereig.begin(), answereig.end() );
// Create a copy
Eigen::SparseMatrix<double> EigMatCopy(this->EigMat);
if ( this->BuiltFromTripletsAtVersion != this->version ) {
EigMatCopy.setFromTriplets( this->triplets.begin(), this->triplets.end() );
}

FloatArray xcopy = x;
Eigen::VectorXd xeig = Eigen::Map<Eigen::VectorXd, Eigen::Unaligned>( xcopy.givePointer(), x.giveSize() );

//Eigen::VectorXd answereig = this->EigMat * xeig;
Eigen::VectorXd answereig = EigMatCopy * xeig; // If set from triplets

answer = FloatArray( answereig.begin(), answereig.end() );
}

void EigenMtrx::createMatrixFromTriplets()
{
if ( this->BuiltFromTripletsAtVersion != this->version ) {
this->EigMat.setFromTriplets( this->triplets.begin(), this->triplets.end() );
this->EigMat.makeCompressed();
BuiltFromTripletsAtVersion = this->version;
}
}

} // end namespace oofem
8 changes: 7 additions & 1 deletion src/oofemlib/eigenmtrx.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -106,6 +106,10 @@ class OOFEM_EXPORT EigenMtrx : public SparseMtrx
std::vector<bool> areFactorized = { false, false, false, false };
int versionUpdate;

std::vector<Eigen::Triplet<double> > triplets; // Allocate vector of triplets
int BuiltFromTripletsAtVersion = 0;


// factorizations
Eigen::SimplicialLLT<Eigen::SparseMatrix<double> > LLT_factorization;
Eigen::SparseLU<Eigen::SparseMatrix<double> > LU_factorization;
Expand Down Expand Up @@ -145,7 +149,7 @@ class OOFEM_EXPORT EigenMtrx : public SparseMtrx
SparseMtrxType giveType() const override { return SMT_EigenMtrx; }
bool isAsymmetric() const override { return true; }

Eigen::SparseMatrix<double> giveMatrix();
Eigen::SparseMatrix<double>& giveMatrix();

//template <typename Derived>
//std::shared_ptr<Eigen::SparseSolverBase<Derived> > giveFactorization( FactorizationType factorizationType );
Expand All @@ -167,6 +171,8 @@ class OOFEM_EXPORT EigenMtrx : public SparseMtrx
bool isFactorized( FactorizationType factorizationType );

void times( const FloatArray &x, FloatArray &answer ) const override;

void createMatrixFromTriplets();
};
} // end namespace oofem
#endif // eigenmtrx_h
Loading

0 comments on commit 10bc90c

Please sign in to comment.