Commit 4fb64953 authored by GILLES Sebastien's avatar GILLES Sebastien
Browse files

#899 Fix a bug in Petsc: reuse case was not properly implemented in MatMatMult...

#899 Fix a bug in Petsc: reuse case was not properly implemented in MatMatMult and similar functions.
parent 75beec38
......@@ -477,7 +477,8 @@ namespace HappyHeart
PtAP(const MatrixT& A,
const MatrixT& P,
MatrixU& out,
const char* invoking_file, int invoking_line);
const char* invoking_file, int invoking_line,
DoReuseMatrix do_reuse_matrix = DoReuseMatrix::no);
......
......@@ -70,18 +70,31 @@ namespace HappyHeart
DoReuseMatrix do_reuse_matrix)
{
Mat result;
if (do_reuse_matrix == DoReuseMatrix::yes)
result = out.Internal();
int error_code = ::MatMatMult(matrix1.Internal(),
matrix2.Internal(),
MAT_INITIAL_MATRIX,
PETSC_DEFAULT,
&result);
if (do_reuse_matrix == DoReuseMatrix::no)
out.SetFromPetscMat(result);
int error_code;
switch(do_reuse_matrix)
{
case DoReuseMatrix::yes:
{
result = out.Internal();
error_code = ::MatMatMult(matrix1.Internal(),
matrix2.Internal(),
MAT_REUSE_MATRIX,
PETSC_DEFAULT,
&result);
break;
}
case DoReuseMatrix::no:
{
error_code = ::MatMatMult(matrix1.Internal(),
matrix2.Internal(),
MAT_INITIAL_MATRIX,
PETSC_DEFAULT,
&result);
out.SetFromPetscMat(result);
break;
}
} // switch
if (error_code)
throw ExceptionNS::Exception(error_code, "MatMatMult", invoking_file, invoking_line);
......@@ -165,19 +178,37 @@ namespace HappyHeart
DoReuseMatrix do_reuse_matrix)
{
Mat result;
if (do_reuse_matrix == DoReuseMatrix::yes)
result = out.Internal();
int error_code = ::MatMatMatMult(matrix1.Internal(),
matrix2.Internal(),
matrix3.Internal(),
MAT_INITIAL_MATRIX,
PETSC_DEFAULT,
&result);
if (do_reuse_matrix == DoReuseMatrix::no)
out.SetFromPetscMat(result);
int error_code;
switch(do_reuse_matrix)
{
case DoReuseMatrix::yes:
{
result = out.Internal();
error_code = ::MatMatMatMult(matrix1.Internal(),
matrix2.Internal(),
matrix3.Internal(),
MAT_REUSE_MATRIX,
PETSC_DEFAULT,
&result);
break;
}
case DoReuseMatrix::no:
{
error_code = ::MatMatMatMult(matrix1.Internal(),
matrix2.Internal(),
matrix3.Internal(),
MAT_INITIAL_MATRIX,
PETSC_DEFAULT,
&result);
out.SetFromPetscMat(result);
break;
}
} // switch
if (error_code)
throw ExceptionNS::Exception(error_code, "MatMatMatMult", invoking_file, invoking_line);
......@@ -340,17 +371,35 @@ namespace HappyHeart
{
Mat result;
if (do_reuse_matrix == DoReuseMatrix::yes)
result = matrix3.Internal();
int error_code = ::MatTransposeMatMult(matrix1.Internal(),
matrix2.Internal(),
MAT_INITIAL_MATRIX,
PETSC_DEFAULT,
&result);
if (do_reuse_matrix == DoReuseMatrix::no)
matrix3.SetFromPetscMat(result);
int error_code;
switch(do_reuse_matrix)
{
case DoReuseMatrix::yes:
{
result = matrix3.Internal();
error_code = ::MatTransposeMatMult(matrix1.Internal(),
matrix2.Internal(),
MAT_REUSE_MATRIX,
PETSC_DEFAULT,
&result);
break;
}
case DoReuseMatrix::no:
{
error_code = ::MatTransposeMatMult(matrix1.Internal(),
matrix2.Internal(),
MAT_INITIAL_MATRIX,
PETSC_DEFAULT,
&result);
matrix3.SetFromPetscMat(result);
break;
}
} // switch
if (error_code)
throw ExceptionNS::Exception(error_code, "MatTransposeMatMult", invoking_file, invoking_line);
......@@ -375,18 +424,33 @@ namespace HappyHeart
DoReuseMatrix do_reuse_matrix)
{
Mat result;
if (do_reuse_matrix == DoReuseMatrix::yes)
result = matrix3.Internal();
int error_code = ::MatMatTransposeMult(matrix1.Internal(),
matrix2.Internal(),
MAT_INITIAL_MATRIX,
PETSC_DEFAULT,
&result);
if (do_reuse_matrix == DoReuseMatrix::no)
matrix3.SetFromPetscMat(result);
int error_code;
switch(do_reuse_matrix)
{
case DoReuseMatrix::yes:
{
result = matrix3.Internal();
error_code = ::MatMatTransposeMult(matrix1.Internal(),
matrix2.Internal(),
MAT_REUSE_MATRIX,
PETSC_DEFAULT,
&result);
break;
}
case DoReuseMatrix::no:
{
error_code = ::MatMatTransposeMult(matrix1.Internal(),
matrix2.Internal(),
MAT_INITIAL_MATRIX,
PETSC_DEFAULT,
&result);
matrix3.SetFromPetscMat(result);
break;
}
} // switch
if (error_code)
throw ExceptionNS::Exception(error_code, "MatMatTransposeMult", invoking_file, invoking_line);
......@@ -408,18 +472,36 @@ namespace HappyHeart
PtAP(const MatrixT& A,
const MatrixT& P,
MatrixU& out,
const char* invoking_file, int invoking_line)
const char* invoking_file, int invoking_line,
DoReuseMatrix do_reuse_matrix)
{
Mat result;
int error_code = ::MatPtAP(A.Internal(),
P.Internal(),
MAT_INITIAL_MATRIX,
PETSC_DEFAULT,
&result);
out.SetFromPetscMat(result);
int error_code;
switch(do_reuse_matrix)
{
case DoReuseMatrix::yes:
{
result = out.Internal();
error_code = ::MatPtAP(A.Internal(),
P.Internal(),
MAT_REUSE_MATRIX,
PETSC_DEFAULT,
&result);
break;
}
case DoReuseMatrix::no:
{
error_code = ::MatPtAP(A.Internal(),
P.Internal(),
MAT_INITIAL_MATRIX,
PETSC_DEFAULT,
&result);
out.SetFromPetscMat(result);
break;
}
} // switch
if (error_code)
throw ExceptionNS::Exception(error_code, "MatPtAP", invoking_file, invoking_line);
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment