Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

FGMRES refinement for Primal Dual System #773

Open
wants to merge 21 commits into
base: stable/3.14
Choose a base branch
from
Open
Changes from all commits
Commits
Show all changes
21 commits
Select commit Hold shift + click to select a range
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
598 changes: 576 additions & 22 deletions src/Algorithm/IpPDFullSpaceSolver.cpp

Large diffs are not rendered by default.

62 changes: 61 additions & 1 deletion src/Algorithm/IpPDFullSpaceSolver.hpp
Original file line number Diff line number Diff line change
@@ -132,6 +132,9 @@ class PDFullSpaceSolver: public PDSystemSolver
*/
Number residual_improvement_factor_;

/** Use (flexible) GMRES refinement instead of iterative refinement */
bool gmres_refinement_;

/** Tolerance for heuristic to ignore wrong inertia */
Number neg_curv_test_tol_;

@@ -208,7 +211,8 @@ class PDFullSpaceSolver: public PDSystemSolver
Number ComputeResidualRatio(
const IteratesVector& rhs,
const IteratesVector& res,
const IteratesVector& resid
const IteratesVector& resid,
Number AInfNrm
);

/** @name Auxiliary functions */
@@ -224,6 +228,62 @@ class PDFullSpaceSolver: public PDSystemSolver
Vector& X
);
///@}

bool SolveGMRES(
Number alpha,
Number beta,
const IteratesVector& rhs,
IteratesVector& res,
bool allow_inexact = false,
bool improve_solution = false
);

bool GMRES(
const SymMatrix& W,
const Matrix& J_c,
const Matrix& J_d,
const Matrix& Px_L,
const Matrix& Px_U,
const Matrix& Pd_L,
const Matrix& Pd_U,
const Vector& z_L,
const Vector& z_U,
const Vector& v_L,
const Vector& v_U,
const Vector& slack_x_L,
const Vector& slack_x_U,
const Vector& slack_s_L,
const Vector& slack_s_U,
const Vector& sigma_x,
const Vector& sigma_s,
const IteratesVector& rhs,
IteratesVector& res,
IteratesVector& resid,
bool improve_solution,
bool resolve_with_better_quality,
bool pretend_singular,
bool& solve_retval
);

Number NrmInf(
const SymMatrix& W,
const Matrix& J_c,
const Matrix& J_d,
const Matrix& Px_L,
const Matrix& Px_U,
const Matrix& Pd_L,
const Matrix& Pd_U,
const Vector& z_L,
const Vector& z_U,
const Vector& v_L,
const Vector& v_U,
const Vector& slack_x_L,
const Vector& slack_x_U,
const Vector& slack_s_L,
const Vector& slack_s_U,
IteratesVector& tmp
);

};

} // namespace Ipopt
112 changes: 112 additions & 0 deletions src/LinAlg/IpCompoundMatrix.cpp
Original file line number Diff line number Diff line change
@@ -715,6 +715,118 @@ void CompoundMatrix::ComputeColAMaxImpl(
}
}

void CompoundMatrix::ComputeRowA1Impl(
Vector& rows_norms,
bool /*init*/
) const
{
if( !matrices_valid_ )
{
matrices_valid_ = MatricesValid();
}
DBG_ASSERT(matrices_valid_);

// The vector is assumed to be compound Vectors as well except if
// there is only one component
CompoundVector* comp_vec = dynamic_cast<CompoundVector*>(&rows_norms);

#ifndef ALLOW_NESTED
// A few sanity checks
if (comp_vec)
{
DBG_ASSERT(NComps_Rows() == comp_vec->NComps());
}
else
{
DBG_ASSERT(NComps_Rows() == 1);
}
#endif
if( comp_vec )
{
if( NComps_Rows() != comp_vec->NComps() )
{
comp_vec = NULL;
}
}

for( Index jcol = 0; jcol < NComps_Cols(); jcol++ )
{
for( Index irow = 0; irow < NComps_Rows(); irow++ )
{
if( ConstComp(irow, jcol) )
{
SmartPtr<Vector> vec_i;
if( comp_vec )
{
vec_i = comp_vec->GetCompNonConst(irow);
}
else
{
vec_i = &rows_norms;
}
DBG_ASSERT(IsValid(vec_i));
ConstComp(irow, jcol)->ComputeRowA1(*vec_i, false);
}
}
}
}

void CompoundMatrix::ComputeColA1Impl(
Vector& cols_norms,
bool /*init*/
) const
{
if( !matrices_valid_ )
{
matrices_valid_ = MatricesValid();
}
DBG_ASSERT(matrices_valid_);

// The vector is assumed to be compound Vectors as well except if
// there is only one component
CompoundVector* comp_vec = dynamic_cast<CompoundVector*>(&cols_norms);

#ifndef ALLOW_NESTED
// A few sanity checks
if (comp_vec)
{
DBG_ASSERT(NComps_Cols() == comp_vec->NComps());
}
else
{
DBG_ASSERT(NComps_Cols() == 1);
}
#endif
if( comp_vec )
{
if( NComps_Cols() != comp_vec->NComps() )
{
comp_vec = NULL;
}
}

for( Index irow = 0; irow < NComps_Rows(); irow++ )
{
for( Index jcol = 0; jcol < NComps_Cols(); jcol++ )
{
if( ConstComp(irow, jcol) )
{
SmartPtr<Vector> vec_j;
if( comp_vec )
{
vec_j = comp_vec->GetCompNonConst(jcol);
}
else
{
vec_j = &cols_norms;
}
DBG_ASSERT(IsValid(vec_j));
ConstComp(irow, jcol)->ComputeColA1(*vec_j, false);
}
}
}
}

void CompoundMatrix::PrintImpl(
const Journalist& jnlst,
EJournalLevel level,
10 changes: 10 additions & 0 deletions src/LinAlg/IpCompoundMatrix.hpp
Original file line number Diff line number Diff line change
@@ -155,6 +155,16 @@ class IPOPTLIB_EXPORT CompoundMatrix: public Matrix
bool init
) const;

virtual void ComputeRowA1Impl(
Vector& rows_norms,
bool init
) const;

virtual void ComputeColA1Impl(
Vector& cols_norms,
bool init
) const;

virtual void PrintImpl(
const Journalist& jnlst,
EJournalLevel level,
51 changes: 51 additions & 0 deletions src/LinAlg/IpCompoundSymMatrix.cpp
Original file line number Diff line number Diff line change
@@ -237,6 +237,57 @@ void CompoundSymMatrix::ComputeRowAMaxImpl(
}
}

void CompoundSymMatrix::ComputeRowA1Impl(
Vector& rows_norms,
bool /*init*/
) const
{
if( !matrices_valid_ )
{
matrices_valid_ = MatricesValid();
}
DBG_ASSERT(matrices_valid_);

// The vector is assumed to be compound Vectors as well except if
// there is only one component
CompoundVector* comp_vec = dynamic_cast<CompoundVector*>(&rows_norms);

// A few sanity checks
if( comp_vec )
{
DBG_ASSERT(NComps_Dim() == comp_vec->NComps());
}
else
{
DBG_ASSERT(NComps_Dim() == 1);
}

for( Index jcol = 0; jcol < NComps_Dim(); jcol++ )
{
for( Index irow = 0; irow < NComps_Dim(); irow++ )
{
SmartPtr<Vector> vec_i;
if( comp_vec )
{
vec_i = comp_vec->GetCompNonConst(irow);
}
else
{
vec_i = &rows_norms;
}
DBG_ASSERT(IsValid(vec_i));
if( jcol <= irow && ConstComp(irow, jcol) )
{
ConstComp(irow, jcol)->ComputeRowA1(*vec_i, false);
}
else if( jcol > irow && ConstComp(jcol, irow) )
{
ConstComp(jcol, irow)->ComputeRowA1(*vec_i, false);
}
}
}
}

void CompoundSymMatrix::PrintImpl(
const Journalist& jnlst,
EJournalLevel level,
5 changes: 5 additions & 0 deletions src/LinAlg/IpCompoundSymMatrix.hpp
Original file line number Diff line number Diff line change
@@ -117,6 +117,11 @@ class IPOPTLIB_EXPORT CompoundSymMatrix: public SymMatrix
bool init
) const;

virtual void ComputeRowA1Impl(
Vector& rows_norms,
bool init
) const;

virtual void PrintImpl(
const Journalist& jnlst,
EJournalLevel level,
16 changes: 16 additions & 0 deletions src/LinAlg/IpDenseGenMatrix.cpp
Original file line number Diff line number Diff line change
@@ -415,6 +415,22 @@ void DenseGenMatrix::ComputeColAMaxImpl(
}
}

void DenseGenMatrix::ComputeRowA1Impl(
Vector& /*rows_norms*/,
bool /*init*/
) const
{
THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "DenseGenMatrix::ComputeRowA1Impl not implemented");
}

void DenseGenMatrix::ComputeColA1Impl(
Vector& /*cols_norms*/,
bool /*init*/
) const
{
THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "DenseGenMatrix::ComputeColA1Impl not implemented");
}

bool DenseGenMatrix::HasValidNumbersImpl() const
{
DBG_ASSERT(initialized_);
10 changes: 10 additions & 0 deletions src/LinAlg/IpDenseGenMatrix.hpp
Original file line number Diff line number Diff line change
@@ -217,6 +217,16 @@ class IPOPTLIB_EXPORT DenseGenMatrix: public Matrix
bool init
) const;

virtual void ComputeRowA1Impl(
Vector& rows_norms,
bool init
) const;

virtual void ComputeColA1Impl(
Vector& cols_norms,
bool init
) const;

virtual void PrintImpl(
const Journalist& jnlst,
EJournalLevel level,
9 changes: 9 additions & 0 deletions src/LinAlg/IpDenseSymMatrix.cpp
Original file line number Diff line number Diff line change
@@ -254,6 +254,15 @@ void DenseSymMatrix::ComputeRowAMaxImpl(
}
}

void DenseSymMatrix::ComputeRowA1Impl(
Vector& /*rows_norms*/,
bool /*init*/
) const
{
THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "DenseSymMatrix::ComputeRowA1Impl not implemented");
}


void DenseSymMatrix::PrintImpl(
const Journalist& jnlst,
EJournalLevel level,
5 changes: 5 additions & 0 deletions src/LinAlg/IpDenseSymMatrix.hpp
Original file line number Diff line number Diff line change
@@ -133,6 +133,11 @@ class IPOPTLIB_EXPORT DenseSymMatrix: public SymMatrix
bool init
) const;

virtual void ComputeRowA1Impl(
Vector& rows_norms,
bool init
) const;

virtual void PrintImpl(
const Journalist& jnlst,
EJournalLevel level,
19 changes: 19 additions & 0 deletions src/LinAlg/IpDiagMatrix.cpp
Original file line number Diff line number Diff line change
@@ -71,6 +71,25 @@ void DiagMatrix::ComputeRowAMaxImpl(
}
}

void DiagMatrix::ComputeRowA1Impl(
Vector& rows_norms,
bool init
) const
{
DBG_ASSERT(IsValid(diag_));
if( init )
{
rows_norms.Copy(*diag_);
rows_norms.ElementWiseAbs();
}
else
{
SmartPtr<Vector> v = diag_->MakeNewCopy();
v->ElementWiseAbs();
rows_norms.Axpy(Number(1.), *v);
}
}

void DiagMatrix::PrintImpl(
const Journalist& jnlst,
EJournalLevel level,
5 changes: 5 additions & 0 deletions src/LinAlg/IpDiagMatrix.hpp
Original file line number Diff line number Diff line change
@@ -63,6 +63,11 @@ class IPOPTLIB_EXPORT DiagMatrix: public SymMatrix
bool init
) const;

virtual void ComputeRowA1Impl(
Vector& rows_norms,
bool init
) const;

virtual void PrintImpl(
const Journalist& jnlst,
EJournalLevel level,
16 changes: 16 additions & 0 deletions src/LinAlg/IpExpandedMultiVectorMatrix.cpp
Original file line number Diff line number Diff line change
@@ -192,6 +192,22 @@ void ExpandedMultiVectorMatrix::ComputeColAMaxImpl(
THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "ExpandedMultiVectorMatrix::ComputeColAMaxImpl not implemented");
}

void ExpandedMultiVectorMatrix::ComputeRowA1Impl(
Vector& /*rows_norms*/,
bool /*init*/
) const
{
THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "ExpandedMultiVectorMatrix::ComputeRowA1Impl not implemented");
}

void ExpandedMultiVectorMatrix::ComputeColA1Impl(
Vector& /*cols_norms*/,
bool /*init*/
) const
{
THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "ExpandedMultiVectorMatrix::ComputeColA1Impl not implemented");
}

void ExpandedMultiVectorMatrix::PrintImpl(
const Journalist& jnlst,
EJournalLevel level,
12 changes: 11 additions & 1 deletion src/LinAlg/IpExpandedMultiVectorMatrix.hpp
Original file line number Diff line number Diff line change
@@ -101,7 +101,17 @@ class ExpandedMultiVectorMatrix: public Matrix
bool init
) const;

virtual void PrintImpl(
virtual void ComputeRowA1Impl(
Vector& rows_norms,
bool init
) const;

virtual void ComputeColA1Impl(
Vector& cols_norms,
bool init
) const;

virtual void PrintImpl(
const Journalist& jnlst,
EJournalLevel level,
EJournalCategory category,
32 changes: 32 additions & 0 deletions src/LinAlg/IpExpansionMatrix.cpp
Original file line number Diff line number Diff line change
@@ -405,6 +405,38 @@ void ExpansionMatrix::ComputeColAMaxImpl(
}
}

void ExpansionMatrix::ComputeRowA1Impl(
Vector& rows_norms,
bool /*init*/
) const
{
DenseVector* dense_vec = static_cast<DenseVector*>(&rows_norms);
DBG_ASSERT(dynamic_cast<DenseVector*>(&rows_norms));
Number* vec_vals = dense_vec->Values();

const Index* exp_pos = ExpandedPosIndices();

for( Index i = 0; i < NCols(); i++ )
{
vec_vals[exp_pos[i]] += Number(1.);
}
}

void ExpansionMatrix::ComputeColA1Impl(
Vector& cols_norms,
bool init
) const
{
if( init )
{
cols_norms.Set(1.);
}
else
{
cols_norms.AddScalar(Number(1.));
}
}

void ExpansionMatrix::PrintImplOffset(
const Journalist& jnlst,
EJournalLevel level,
10 changes: 10 additions & 0 deletions src/LinAlg/IpExpansionMatrix.hpp
Original file line number Diff line number Diff line change
@@ -103,6 +103,16 @@ class IPOPTLIB_EXPORT ExpansionMatrix: public Matrix
bool init
) const;

virtual void ComputeRowA1Impl(
Vector& rows_norms,
bool init
) const;

virtual void ComputeColA1Impl(
Vector& cols_norms,
bool init
) const;

virtual void PrintImpl(
const Journalist& jnlst,
EJournalLevel level,
15 changes: 15 additions & 0 deletions src/LinAlg/IpIdentityMatrix.cpp
Original file line number Diff line number Diff line change
@@ -62,6 +62,21 @@ void IdentityMatrix::ComputeRowAMaxImpl(
}
}

void IdentityMatrix::ComputeRowA1Impl(
Vector& rows_norms,
bool init
) const
{
if( init )
{
rows_norms.Set(1.);
}
else
{
rows_norms.AddScalar(Number(1.));
}
}

void IdentityMatrix::PrintImpl(
const Journalist& jnlst,
EJournalLevel level,
5 changes: 5 additions & 0 deletions src/LinAlg/IpIdentityMatrix.hpp
Original file line number Diff line number Diff line change
@@ -71,6 +71,11 @@ class IPOPTLIB_EXPORT IdentityMatrix: public SymMatrix
bool init
) const;

virtual void ComputeRowA1Impl(
Vector& rows_norms,
bool init
) const;

virtual void PrintImpl(
const Journalist& jnlst,
EJournalLevel level,
16 changes: 16 additions & 0 deletions src/LinAlg/IpLowRankUpdateSymMatrix.cpp
Original file line number Diff line number Diff line change
@@ -170,6 +170,22 @@ void LowRankUpdateSymMatrix::ComputeColAMaxImpl(
THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "LowRankUpdateSymMatrix::ComputeColAMaxImpl not implemented");
}

void LowRankUpdateSymMatrix::ComputeRowA1Impl(
Vector& /*rows_norms*/,
bool /*init*/
) const
{
THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "LowRankUpdateSymMatrix::ComputeRowA1Impl not implemented");
}

void LowRankUpdateSymMatrix::ComputeColA1Impl(
Vector& /*cols_norms*/,
bool /*init*/
) const
{
THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "LowRankUpdateSymMatrix::ComputeColA1Impl not implemented");
}

void LowRankUpdateSymMatrix::PrintImpl(
const Journalist& jnlst,
EJournalLevel level,
12 changes: 11 additions & 1 deletion src/LinAlg/IpLowRankUpdateSymMatrix.hpp
Original file line number Diff line number Diff line change
@@ -122,7 +122,17 @@ class LowRankUpdateSymMatrix: public SymMatrix
bool init
) const;

virtual void PrintImpl(
virtual void ComputeRowA1Impl(
Vector& rows_norms,
bool init
) const;

virtual void ComputeColA1Impl(
Vector& cols_norms,
bool init
) const;

virtual void PrintImpl(
const Journalist& jnlst,
EJournalLevel level,
EJournalCategory category,
56 changes: 56 additions & 0 deletions src/LinAlg/IpMatrix.hpp
Original file line number Diff line number Diff line change
@@ -163,6 +163,44 @@ class IPOPTLIB_EXPORT Matrix: public TaggedObject
}
///@}

/** Compute the 1-norm of the rows in the matrix.
*
* The result is stored in rows_norms.
* The vector is assumed to be initialized if init is false.
*/
void ComputeRowA1(
Vector& rows_norms,
bool init = true
) const
{
DBG_ASSERT(NRows() == rows_norms.Dim());
if( init )
{
rows_norms.Set(0.);
}
ComputeRowA1Impl(rows_norms, init);
}
///@}

/** Compute the 1-norm of the columns in the matrix.
*
* The result is stored in cols_norms.
* The vector is assumed to be initialized if init is false.
*/
void ComputeColA1(
Vector& cols_norms,
bool init = true
) const
{
DBG_ASSERT(NCols() == cols_norms.Dim());
if( init )
{
cols_norms.Set(0.);
}
ComputeColA1Impl(cols_norms, init);
}
///@}

/** Print detailed information about the matrix.
*
* @attention Do not overload. Overload PrintImpl instead.
@@ -270,6 +308,24 @@ class IPOPTLIB_EXPORT Matrix: public TaggedObject
bool init
) const = 0;

/** Compute the 1-norm of the rows in the matrix.
*
* The result is stored in rows_norms.
* The vector is assumed to be initialized if init is false. */
virtual void ComputeRowA1Impl(
Vector& rows_norms,
bool init
) const = 0;

/** Compute the 1-norm of the cols in the matrix.
*
* The result is stored in cols_norms.
* The vector is assumed to be initialized if init is false. */
virtual void ComputeColA1Impl(
Vector& cols_norms,
bool init
) const = 0;

/** Print detailed information about the matrix. */
virtual void PrintImpl(
const Journalist& jnlst,
16 changes: 16 additions & 0 deletions src/LinAlg/IpMultiVectorMatrix.cpp
Original file line number Diff line number Diff line change
@@ -301,6 +301,22 @@ void MultiVectorMatrix::ComputeColAMaxImpl(
THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "MultiVectorMatrix::ComputeColAMaxImpl not implemented");
}

void MultiVectorMatrix::ComputeRowA1Impl(
Vector& /*rows_norms*/,
bool /*init*/
) const
{
THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "MultiVectorMatrix::ComputeRowA1Impl not implemented");
}

void MultiVectorMatrix::ComputeColA1Impl(
Vector& /*cols_norms*/,
bool /*init*/
) const
{
THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "MultiVectorMatrix::ComputeColA1Impl not implemented");
}

void MultiVectorMatrix::PrintImpl(
const Journalist& jnlst,
EJournalLevel level,
10 changes: 10 additions & 0 deletions src/LinAlg/IpMultiVectorMatrix.hpp
Original file line number Diff line number Diff line change
@@ -165,6 +165,16 @@ class IPOPTLIB_EXPORT MultiVectorMatrix: public Matrix
bool init
) const;

virtual void ComputeRowA1Impl(
Vector& rows_norms,
bool init
) const;

virtual void ComputeColA1Impl(
Vector& cols_norms,
bool init
) const;

virtual void PrintImpl(
const Journalist& jnlst,
EJournalLevel level,
79 changes: 79 additions & 0 deletions src/LinAlg/IpScaledMatrix.cpp
Original file line number Diff line number Diff line change
@@ -5,6 +5,11 @@
// Authors: Carl Laird, Andreas Waechter IBM 2004-08-13

#include "IpScaledMatrix.hpp"
#include "IpGenTMatrix.hpp"
#include "IpDenseVector.hpp"

#include <iostream>
#include <typeinfo>

namespace Ipopt
{
@@ -117,6 +122,80 @@ void ScaledMatrix::ComputeColAMaxImpl(
THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "ScaledMatrix::ComputeColAMaxImpl not implemented");
}

void ScaledMatrix::ComputeRowA1Impl(
Vector& rows_norms,
bool init
) const
{
DBG_ASSERT(IsValid(matrix_));

if( IsValid(owner_space_->ColumnScaling()) )
{
const GenTMatrix* mt = dynamic_cast<const GenTMatrix*>(GetRawPtr(matrix_));
if ( mt )
{
const DenseVector* Dc = dynamic_cast<const DenseVector*>(GetRawPtr(owner_space_->ColumnScaling()));
DenseVector* rn = dynamic_cast<DenseVector*>(&rows_norms);
DBG_ASSERT(Dc);
DBG_ASSERT(rn);
for ( Index i = 0; i < mt->Nonzeros(); i++ )
{
rn->Values()[mt->Irows()[i] - 1] += std::abs(mt->Values()[i] * Dc->Values()[mt->Jcols()[i] - 1]);
}
}
else
{
THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "ScaledMatrix::ComputeRowA1Impl not implemented");
}
}
else
{
matrix_->ComputeRowA1(rows_norms, init);
}
if( IsValid(owner_space_->RowScaling()) )
{
rows_norms.ElementWiseMultiply(*owner_space_->RowScaling());
rows_norms.ElementWiseAbs();
}
}

void ScaledMatrix::ComputeColA1Impl(
Vector& cols_norms,
bool init
) const
{
DBG_ASSERT(IsValid(matrix_));

if( IsValid(owner_space_->RowScaling()) )
{
const GenTMatrix* mt = dynamic_cast<const GenTMatrix*>(GetRawPtr(matrix_));
if ( mt )
{
const DenseVector* Dr = dynamic_cast<const DenseVector*>(GetRawPtr(owner_space_->RowScaling()));
DenseVector* cn = dynamic_cast<DenseVector*>(&cols_norms);
DBG_ASSERT(Dr);
DBG_ASSERT(cn);
for ( Index i = 0; i < mt->Nonzeros(); i++ )
{
cn->Values()[mt->Jcols()[i] - 1] += std::abs(mt->Values()[i] * Dr->Values()[mt->Irows()[i] - 1]);
}
}
else
{
THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "ScaledMatrix::ComputeColA1Impl not implemented");
}
}
else
{
matrix_->ComputeColA1(cols_norms, init);
}
if( IsValid(owner_space_->ColumnScaling()) )
{
cols_norms.ElementWiseMultiply(*owner_space_->ColumnScaling());
cols_norms.ElementWiseAbs();
}
}

void ScaledMatrix::PrintImpl(
const Journalist& jnlst,
EJournalLevel level,
10 changes: 10 additions & 0 deletions src/LinAlg/IpScaledMatrix.hpp
Original file line number Diff line number Diff line change
@@ -89,6 +89,16 @@ class IPOPTLIB_EXPORT ScaledMatrix: public Matrix
bool init
) const;

virtual void ComputeRowA1Impl(
Vector& rows_norms,
bool init
) const;

virtual void ComputeColA1Impl(
Vector& cols_norms,
bool init
) const;

virtual void PrintImpl(
const Journalist& jnlst,
EJournalLevel level,
16 changes: 16 additions & 0 deletions src/LinAlg/IpSumMatrix.cpp
Original file line number Diff line number Diff line change
@@ -136,6 +136,22 @@ void SumMatrix::ComputeColAMaxImpl(
THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "SumMatrix::ComputeColAMaxImpl not implemented");
}

void SumMatrix::ComputeRowA1Impl(
Vector& /*rows_norms*/,
bool /*init*/
) const
{
THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "SumMatrix::ComputeRowA1Impl not implemented");
}

void SumMatrix::ComputeColA1Impl(
Vector& /*cols_norms*/,
bool /*init*/
) const
{
THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "SumMatrix::ComputeColA1Impl not implemented");
}

void SumMatrix::PrintImpl(
const Journalist& jnlst,
EJournalLevel level,
10 changes: 10 additions & 0 deletions src/LinAlg/IpSumMatrix.hpp
Original file line number Diff line number Diff line change
@@ -83,6 +83,16 @@ class SumMatrix: public Matrix
bool init
) const;

virtual void ComputeRowA1Impl(
Vector& rows_norms,
bool init
) const;

virtual void ComputeColA1Impl(
Vector& cols_norms,
bool init
) const;

virtual void PrintImpl(
const Journalist& jnlst,
EJournalLevel level,
24 changes: 24 additions & 0 deletions src/LinAlg/IpSumSymMatrix.cpp
Original file line number Diff line number Diff line change
@@ -108,6 +108,30 @@ void SumSymMatrix::ComputeColAMaxImpl(
THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "SumSymMatrix::ComputeColAMaxImpl not implemented");
}

void SumSymMatrix::ComputeRowA1Impl(
Vector& rows_norms,
bool /*init*/
) const
{
for( Index iterm = 0; iterm < NTerms(); iterm++ )
{
DBG_ASSERT(IsValid(matrices_[iterm]));
matrices_[iterm]->ComputeRowA1(rows_norms, false);
}
}

void SumSymMatrix::ComputeColA1Impl(
Vector& cols_norms,
bool /*init*/
) const
{
for( Index iterm = 0; iterm < NTerms(); iterm++ )
{
DBG_ASSERT(IsValid(matrices_[iterm]));
matrices_[iterm]->ComputeColA1(cols_norms, false);
}
}

void SumSymMatrix::PrintImpl(
const Journalist& jnlst,
EJournalLevel level,
10 changes: 10 additions & 0 deletions src/LinAlg/IpSumSymMatrix.hpp
Original file line number Diff line number Diff line change
@@ -81,6 +81,16 @@ class IPOPTLIB_EXPORT SumSymMatrix: public SymMatrix
bool init
) const;

virtual void ComputeRowA1Impl(
Vector& rows_norms,
bool init
) const;

virtual void ComputeColA1Impl(
Vector& cols_norms,
bool init
) const;

virtual void PrintImpl(
const Journalist& jnlst,
EJournalLevel level,
13 changes: 13 additions & 0 deletions src/LinAlg/IpSymMatrix.hpp
Original file line number Diff line number Diff line change
@@ -74,6 +74,19 @@ class IPOPTLIB_EXPORT SymMatrix: public Matrix
}
///@}

/** Implementation of ComputeColA1Impl, which calls ComputeRowA1Impl.
*
* Since the matrix is symmetric, the row and column 1 norms are identical.
*/
virtual void ComputeColA1Impl(
Vector& cols_norms,
bool init
) const
{
ComputeRowA1Impl(cols_norms, init);
}
///@}

private:
/** Copy of the owner space ptr as a SymMatrixSpace instead
* of a MatrixSpace
38 changes: 38 additions & 0 deletions src/LinAlg/IpSymScaledMatrix.cpp
Original file line number Diff line number Diff line change
@@ -5,6 +5,8 @@
// Authors: Carl Laird, Andreas Waechter IBM 2004-08-13

#include "IpSymScaledMatrix.hpp"
#include "IpSymTMatrix.hpp"
#include "IpDenseVector.hpp"

namespace Ipopt
{
@@ -70,6 +72,42 @@ void SymScaledMatrix::ComputeRowAMaxImpl(
THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "SymScaledMatrix::ComputeRowAMaxImpl not implemented");
}

void SymScaledMatrix::ComputeRowA1Impl(
Vector& rows_norms,
bool init
) const
{
DBG_ASSERT(IsValid(matrix_));

if( IsValid(owner_space_->RowColScaling()) )
{
const SymTMatrix* mt = dynamic_cast<const SymTMatrix*>(GetRawPtr(matrix_));
if( mt )
{
const DenseVector* D = dynamic_cast<const DenseVector*>(GetRawPtr(owner_space_->RowColScaling()));
DenseVector* rn = dynamic_cast<DenseVector*>(&rows_norms);
DBG_ASSERT(D);
DBG_ASSERT(rn);
for ( Index i = 0; i < mt->Nonzeros(); i++ )
{
rn->Values()[mt->Irows()[i] - 1] += std::abs(mt->Values()[i] * D->Values()[mt->Jcols()[i] - 1]);
if( mt->Irows()[i] != mt->Jcols()[i] )
{
rn->Values()[mt->Jcols()[i] - 1] += std::abs(mt->Values()[i] * D->Values()[mt->Irows()[i] - 1]);
}
}
}
else
{
THROW_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED, "ScaledMatrix::ComputeRowA1Impl not implemented");
}
}
else
{
matrix_->ComputeRowA1(rows_norms, init);
}
}

void SymScaledMatrix::PrintImpl(
const Journalist& jnlst,
EJournalLevel level,
5 changes: 5 additions & 0 deletions src/LinAlg/IpSymScaledMatrix.hpp
Original file line number Diff line number Diff line change
@@ -74,6 +74,11 @@ class IPOPTLIB_EXPORT SymScaledMatrix: public SymMatrix
bool init
) const;

virtual void ComputeRowA1Impl(
Vector& rows_norms,
bool init
) const;

virtual void PrintImpl(
const Journalist& jnlst,
EJournalLevel level,
16 changes: 16 additions & 0 deletions src/LinAlg/IpTransposeMatrix.hpp
Original file line number Diff line number Diff line change
@@ -85,6 +85,22 @@ class TransposeMatrix: public Matrix
orig_matrix_->ComputeRowAMax(rows_norms, init);
}

virtual void ComputeRowA1Impl(
Vector& rows_norms,
bool init
) const {
DBG_ASSERT(IsValid(orig_matrix_));
orig_matrix_->ComputeColA1(rows_norms, init);
}

virtual void ComputeColA1Impl(
Vector& cols_norms,
bool init
) const {
DBG_ASSERT(IsValid(orig_matrix_));
orig_matrix_->ComputeRowA1(cols_norms, init);
}

virtual void PrintImpl(
const Journalist& jnlst,
EJournalLevel level,
12 changes: 12 additions & 0 deletions src/LinAlg/IpZeroMatrix.hpp
Original file line number Diff line number Diff line change
@@ -58,6 +58,18 @@ class ZeroMatrix: public Matrix
) const
{ }

virtual void ComputeRowA1Impl(
Vector& /*rows_norms*/,
bool /*init*/
) const
{ }

virtual void ComputeColA1Impl(
Vector& /*cols_norms*/,
bool /*init*/
) const
{ }

virtual void PrintImpl(
const Journalist& jnlst,
EJournalLevel level,
12 changes: 12 additions & 0 deletions src/LinAlg/IpZeroSymMatrix.hpp
Original file line number Diff line number Diff line change
@@ -57,6 +57,18 @@ class IPOPTLIB_EXPORT ZeroSymMatrix: public SymMatrix
) const
{ }

virtual void ComputeRowA1Impl(
Vector& /*rows_norms*/,
bool /*init*/
) const
{ }

virtual void ComputeColA1Impl(
Vector& /*cols_norms*/,
bool /*init*/
) const
{ }

virtual void PrintImpl(
const Journalist& jnlst,
EJournalLevel level,
54 changes: 54 additions & 0 deletions src/LinAlg/TMatrices/IpGenTMatrix.cpp
Original file line number Diff line number Diff line change
@@ -238,6 +238,60 @@ void GenTMatrix::ComputeColAMaxImpl(
}
}

void GenTMatrix::ComputeRowA1Impl(
Vector& rows_norms,
bool /*init*/
) const
{
DBG_ASSERT(initialized_);

if( NRows() == 0 )
{
return;
}

DenseVector* dense_vec = static_cast<DenseVector*>(&rows_norms);
DBG_ASSERT(dynamic_cast<DenseVector*>(&rows_norms));

const Index* irows = Irows();
const Number* val = values_;
Number* vec_vals = dense_vec->Values();
DBG_ASSERT(vec_vals != NULL);

vec_vals--;
for( Index i = 0; i < Nonzeros(); i++ )
{
vec_vals[irows[i]] += std::abs(val[i]);
}
}

void GenTMatrix::ComputeColA1Impl(
Vector& cols_norms,
bool /*init*/
) const
{
DBG_ASSERT(initialized_);

if( NCols() == 0 )
{
return;
}

DenseVector* dense_vec = static_cast<DenseVector*>(&cols_norms);
DBG_ASSERT(dynamic_cast<DenseVector*>(&cols_norms));

const Index* jcols = Jcols();
const Number* val = values_;
Number* vec_vals = dense_vec->Values();
DBG_ASSERT(vec_vals != NULL);

vec_vals--; // to deal with 1-based indexing in jcols, I believe
for( Index i = 0; i < Nonzeros(); i++ )
{
vec_vals[jcols[i]] += std::abs(val[i]);
}
}

void GenTMatrix::PrintImplOffset(
const Journalist& jnlst,
EJournalLevel level,
10 changes: 10 additions & 0 deletions src/LinAlg/TMatrices/IpGenTMatrix.hpp
Original file line number Diff line number Diff line change
@@ -119,6 +119,16 @@ class IPOPTLIB_EXPORT GenTMatrix: public Matrix
bool init
) const;

virtual void ComputeRowA1Impl(
Vector& rows_norms,
bool init
) const;

virtual void ComputeColA1Impl(
Vector& cols_norms,
bool init
) const;

virtual void PrintImpl(
const Journalist& jnlst,
EJournalLevel level,
36 changes: 36 additions & 0 deletions src/LinAlg/TMatrices/IpSymTMatrix.cpp
Original file line number Diff line number Diff line change
@@ -197,6 +197,42 @@ void SymTMatrix::ComputeRowAMaxImpl(
}
}

void SymTMatrix::ComputeRowA1Impl(
Vector& rows_norms,
bool /*init*/
) const
{
DBG_ASSERT(initialized_);

if( NRows() == 0 )
{
return;
}

DenseVector* dense_vec = static_cast<DenseVector*>(&rows_norms);
DBG_ASSERT(dynamic_cast<DenseVector*>(&rows_norms));

const Index* irn = Irows();
const Index* jcn = Jcols();
const Number* val = values_;
Number* vec_vals = dense_vec->Values();
DBG_ASSERT(vec_vals != NULL);

// const Number zero = 0.;
// IpBlasCopy(NRows(), &zero, 0, vec_vals, 1);

vec_vals--; // to deal with 1-based indexing in irn and jcn (I believe)
for( Index i = 0; i < Nonzeros(); i++ )
{
const Number f = std::abs(*val);
vec_vals[*irn] += f;
vec_vals[*jcn] += f;
val++;
irn++;
jcn++;
}
}

void SymTMatrix::PrintImpl(
const Journalist& jnlst,
EJournalLevel level,
5 changes: 5 additions & 0 deletions src/LinAlg/TMatrices/IpSymTMatrix.hpp
Original file line number Diff line number Diff line change
@@ -129,6 +129,11 @@ class IPOPTLIB_EXPORT SymTMatrix: public SymMatrix
bool init
) const;

virtual void ComputeRowA1Impl(
Vector& rows_norms,
bool init
) const;

virtual void PrintImpl(
const Journalist& jnlst,
EJournalLevel level,