Skip to content

Commit 2c2f013

Browse files
dcwhitedcwhite
authored
Merge pull request #1606 from SCIInstitute/solve-complex-system
Solve complex system--part 1
2 parents b9afe85 + 3b9b017 commit 2c2f013

17 files changed

+670
-139
lines changed

‎src/Core/Algorithms/Math/ConvertMatrixType.cc‎

Lines changed: 22 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -41,7 +41,7 @@ using namespace SCIRun::Core::Datatypes;
4141
ALGORITHM_PARAMETER_DEF(Math, OutputMatrixType);
4242

4343
/// @class ConvertMatrixType
44-
/// @image html ConvertMatrixType.png
44+
/// @image html ConvertMatrixType.png
4545
ConvertMatrixTypeAlgorithm::ConvertMatrixTypeAlgorithm()
4646
{
4747
addOption(Parameters::OutputMatrixType, "passThrough", "passThrough|dense|column|sparse");
@@ -51,7 +51,7 @@ MatrixHandle ConvertMatrixTypeAlgorithm::run(MatrixHandle input_matrix) const
5151
{
5252
if (!input_matrix)
5353
{
54-
THROW_ALGORITHM_INPUT_ERROR("No input matrix");
54+
THROW_ALGORITHM_INPUT_ERROR("No input matrix");
5555
}
5656

5757
MatrixHandle omH;
@@ -63,62 +63,62 @@ MatrixHandle ConvertMatrixTypeAlgorithm::run(MatrixHandle input_matrix) const
6363
if(matrixIs::dense(input_matrix))
6464
{
6565
ostr2 << "Input Matrix Type: DENSE MATRIX";
66-
}
66+
}
6767
else if (matrixIs::column(input_matrix))
6868
{
6969
ostr2 << "Input Matrix Type: COLUMN MATRIX";
70-
}
70+
}
7171
else if (matrixIs::sparse(input_matrix))
7272
{
7373
ostr2 << "Input Matrix Type: SPARSE MATRIX";
74-
}
75-
else
74+
}
75+
else
7676
{
77-
THROW_ALGORITHM_INPUT_ERROR("Unknown input matrix type");
77+
THROW_ALGORITHM_INPUT_ERROR("Unknown input matrix type");
7878
}
7979

8080
remark(ostr2.str());
8181

8282
auto outputType = getOption(Parameters::OutputMatrixType);
8383

84-
if ("passThrough" == outputType)
84+
if ("passThrough" == outputType)
8585
{
8686
return input_matrix;
87-
}
87+
}
8888
else
8989
{
9090
if ("column" == outputType && !matrixIs::column(input_matrix))
9191
{
9292
if (input_matrix->ncols()!=1)
9393
{
94-
THROW_ALGORITHM_INPUT_ERROR("Input matrix needs to have a single column to be converted to column matrix type.");
94+
THROW_ALGORITHM_INPUT_ERROR("Input matrix needs to have a single column to be converted to column matrix type.");
9595
}
96-
DenseColumnMatrixHandle output = convertMatrix::toColumn(input_matrix);
97-
if (!output)
96+
auto output = convertMatrix::toColumn(input_matrix);
97+
if (!output)
9898
{
99-
THROW_ALGORITHM_INPUT_ERROR("Conversion to column matrix failed.");
99+
THROW_ALGORITHM_INPUT_ERROR("Conversion to column matrix failed.");
100100
}
101101
return output;
102-
}
102+
}
103103
else if ("dense" == outputType && !matrixIs::dense(input_matrix))
104104
{
105105
auto output = convertMatrix::toDense(input_matrix);
106-
if (!output)
106+
if (!output)
107107
{
108-
THROW_ALGORITHM_INPUT_ERROR("Conversion to dense matrix failed.");
108+
THROW_ALGORITHM_INPUT_ERROR("Conversion to dense matrix failed.");
109109
}
110110
return output;
111-
}
111+
}
112112
else if ("sparse" == outputType && !matrixIs::sparse(input_matrix))
113113
{
114114
auto output = convertMatrix::toSparse(input_matrix);
115-
if (!output)
115+
if (!output)
116116
{
117-
THROW_ALGORITHM_INPUT_ERROR("Conversion to sparse matrix failed.");
117+
THROW_ALGORITHM_INPUT_ERROR("Conversion to sparse matrix failed.");
118118
}
119119
return output;
120-
}
121-
{
120+
}
121+
{
122122
remark(" Datatype unknown or input and output data type are equal. Passing input matrix through. ");
123123
return input_matrix;
124124
}
@@ -131,7 +131,7 @@ AlgorithmOutput ConvertMatrixTypeAlgorithm::run(const AlgorithmInput& input) con
131131

132132
MatrixHandle output_matrix = run(input_matrix);
133133

134-
AlgorithmOutput output;
134+
AlgorithmOutput output;
135135
output[Variables::ResultMatrix] = output_matrix;
136136

137137
return output;

‎src/Core/Algorithms/Math/LinearSystem/SolveLinearSystemAlgo.cc‎

Lines changed: 4 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -1091,7 +1091,7 @@ bool SolveLinearSystemAlgo::run(SparseRowMatrixHandle A,
10911091
ENSURE_POSITIVE_DOUBLE(tolerance, "Tolerance out of range!");
10921092
ENSURE_POSITIVE_INT(maxIterations, "Max iterations out of range!");
10931093

1094-
if (!matrixIs::sparse(A))
1094+
if (false/*!matrixIs::sparse(A)*/)
10951095
{
10961096
THROW_ALGORITHM_INPUT_ERROR("Matrix A is not sparse");
10971097
#ifdef SCIRUN4_CODE_TO_BE_ENABLED_LATER
@@ -1104,7 +1104,7 @@ bool SolveLinearSystemAlgo::run(SparseRowMatrixHandle A,
11041104
#endif
11051105
}
11061106

1107-
if (!matrixIs::dense(b) && !matrixIs::column(b))
1107+
if (false)
11081108
{
11091109
THROW_ALGORITHM_INPUT_ERROR("Matrix b is not a dense or column matrix");
11101110
#ifdef SCIRUN4_CODE_TO_BE_ENABLED_LATER
@@ -1120,16 +1120,11 @@ bool SolveLinearSystemAlgo::run(SparseRowMatrixHandle A,
11201120
if (!x0)
11211121
{
11221122
// create an x0 matrix
1123-
DenseColumnMatrixHandle temp(boost::make_shared<DenseColumnMatrix>(b->nrows()));
1123+
auto temp(boost::make_shared<DenseColumnMatrix>(b->nrows()));
11241124
temp->setZero();
11251125
x0 = temp;
11261126
}
1127-
1128-
if (!matrixIs::dense(x0) && !matrixIs::column(x0))
1129-
{
1130-
THROW_ALGORITHM_INPUT_ERROR("Matrix x0 is not a dense or column matrix");
1131-
}
1132-
1127+
11331128
if ((x0->ncols() != 1) || (b->ncols() != 1))
11341129
{
11351130
THROW_ALGORITHM_INPUT_ERROR("Matrix x0 and b need to have the same number of rows");

‎src/Core/Algorithms/Math/ParallelAlgebra/ParallelLinearAlgebra.cc‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -89,7 +89,6 @@ bool ParallelLinearAlgebra::add_vector(DenseColumnMatrixHandle mat, ParallelVect
8989
if (!mat) { return (false); }
9090
if (mat->ncols() != 1) { return (false); }
9191
if (mat->nrows() != size_) { return (false); }
92-
if (matrixIs::sparse(mat)) { return (false); }
9392

9493
V.data_ = mat->data();
9594
V.size_ = size_;

‎src/Core/Algorithms/Math/SolveLinearSystemWithEigen.cc‎

Lines changed: 36 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -41,20 +41,23 @@ using namespace SCIRun::Core;
4141

4242
namespace
4343
{
44+
template <class ColumnMatrixType>
4445
class SolveLinearSystemAlgorithmEigenCGImpl
4546
{
4647
public:
47-
SolveLinearSystemAlgorithmEigenCGImpl(const DenseColumnMatrix& rhs, double tolerance, int maxIterations) :
48-
rhs_(rhs), tolerance_(tolerance), maxIterations_(maxIterations) {}
48+
SolveLinearSystemAlgorithmEigenCGImpl(const ColumnMatrixType& rhs, double tolerance, int maxIterations) :
49+
tolerance_(tolerance), maxIterations_(maxIterations), rhs_(rhs) {}
50+
51+
using SolutionType = ColumnMatrixType;
4952

5053
template <class MatrixType>
51-
DenseColumnMatrix::EigenBase solveWithEigen(const MatrixType& lhs)
54+
typename ColumnMatrixType::EigenBase solveWithEigen(const MatrixType& lhs)
5255
{
5356
Eigen::ConjugateGradient<typename MatrixType::EigenBase> cg;
5457
cg.compute(lhs);
5558

5659
if (cg.info() != Eigen::Success)
57-
BOOST_THROW_EXCEPTION(AlgorithmInputException()
60+
BOOST_THROW_EXCEPTION(AlgorithmInputException()
5861
<< LinearAlgebraErrorMessage("Conjugate gradient initialization was unsuccessful")
5962
<< EigenComputationInfo(cg.info()));
6063

@@ -69,42 +72,57 @@ namespace
6972
double tolerance_;
7073
int maxIterations_;
7174
private:
72-
const DenseColumnMatrix& rhs_;
75+
const ColumnMatrixType& rhs_;
7376
};
7477
}
7578

7679
SolveLinearSystemAlgorithm::Outputs SolveLinearSystemAlgorithm::run(const Inputs& input, const Parameters& params) const
7780
{
78-
auto A = input.get<0>();
81+
return runImpl<Inputs, Outputs>(input, params);
82+
}
83+
84+
SolveLinearSystemAlgorithm::ComplexOutputs SolveLinearSystemAlgorithm::run(const ComplexInputs& input, const Parameters& params) const
85+
{
86+
return runImpl<ComplexInputs, ComplexOutputs>(input, params);
87+
}
88+
89+
template <typename In, typename Out>
90+
Out SolveLinearSystemAlgorithm::runImpl(const In& input, const Parameters& params) const
91+
{
92+
auto A = std::get<0>(input);
7993
ENSURE_ALGORITHM_INPUT_NOT_NULL(A, "Null input matrix");
8094

81-
auto b = input.get<1>();
95+
auto b = std::get<1>(input);
8296
ENSURE_ALGORITHM_INPUT_NOT_NULL(b, "Null rhs vector");
83-
84-
double tolerance = params.get<0>();
97+
98+
double tolerance = std::get<0>(params);
8599
ENSURE_POSITIVE_DOUBLE(tolerance, "Tolerance out of range!");
86100

87-
int maxIterations = params.get<1>();
101+
int maxIterations = std::get<1>(params);
88102
ENSURE_POSITIVE_INT(maxIterations, "Max iterations out of range!");
89103

90-
SolveLinearSystemAlgorithmEigenCGImpl impl(*b, tolerance, maxIterations);
91-
DenseColumnMatrix x;
104+
using SolutionType = DenseMatrixGeneric<typename std::tuple_element<0, In>::type::element_type::value_type>;
105+
using SolverType = SolveLinearSystemAlgorithmEigenCGImpl<SolutionType>;
106+
SolverType impl(*b, tolerance, maxIterations);
107+
typename SolverType::SolutionType x;
92108
if (matrixIs::dense(A))
93109
{
94-
x = impl.solveWithEigen(*castMatrix::toDense(A));
110+
auto dense = castMatrix::toDense(A);
111+
x = impl.solveWithEigen(*dense);
95112
}
96113
else if (matrixIs::sparse(A))
97114
{
98-
x = impl.solveWithEigen(*castMatrix::toSparse(A));
115+
auto sparse = castMatrix::toSparse(A);
116+
x = impl.solveWithEigen(*sparse);
99117
}
100118
else
101119
BOOST_THROW_EXCEPTION(AlgorithmProcessingException() << ErrorMessage("solveWithEigen can only handle dense and sparse matrices."));
102-
120+
103121
if (x.size() != 0)
104122
{
105123
/// @todo: move ctor
106-
DenseColumnMatrixHandle solution(boost::make_shared<DenseColumnMatrix>(x));
107-
return SolveLinearSystemAlgorithm::Outputs(solution, impl.tolerance_, impl.maxIterations_);
124+
auto solution(boost::make_shared<SolutionType>(x));
125+
return Out(solution, impl.tolerance_, impl.maxIterations_);
108126
}
109127
else
110128
BOOST_THROW_EXCEPTION(AlgorithmProcessingException() << ErrorMessage("solveWithEigen produced an empty solution."));
@@ -113,4 +131,4 @@ SolveLinearSystemAlgorithm::Outputs SolveLinearSystemAlgorithm::run(const Inputs
113131
AlgorithmOutput SolveLinearSystemAlgorithm::run(const AlgorithmInput& input) const
114132
{
115133
throw 2;
116-
}
134+
}

‎src/Core/Algorithms/Math/SolveLinearSystemWithEigen.h‎

Lines changed: 12 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -38,22 +38,28 @@ namespace SCIRun {
3838
namespace Core {
3939
namespace Algorithms {
4040
namespace Math {
41-
41+
4242
/// @todo: this will be the base class of all the solvers. for now it will just contain the Eigen CG impl.
4343
class SCISHARE SolveLinearSystemAlgorithm : public AlgorithmBase
4444
{
4545
public:
46-
typedef boost::tuple<SCIRun::Core::Datatypes::MatrixHandle, SCIRun::Core::Datatypes::DenseColumnMatrixHandle> Inputs;
47-
typedef boost::tuple<double, int> Parameters;
48-
typedef boost::tuple<SCIRun::Core::Datatypes::DenseColumnMatrixHandle, double, int> Outputs;
46+
typedef std::tuple<SCIRun::Core::Datatypes::MatrixHandle, SCIRun::Core::Datatypes::DenseColumnMatrixHandle> Inputs;
47+
typedef std::tuple<SCIRun::Core::Datatypes::ComplexMatrixHandle, SCIRun::Core::Datatypes::ComplexDenseColumnMatrixHandle> ComplexInputs;
48+
typedef std::tuple<double, int> Parameters;
49+
typedef std::tuple<SCIRun::Core::Datatypes::DenseMatrixHandle, double, int> Outputs;
50+
typedef std::tuple<SCIRun::Core::Datatypes::ComplexDenseMatrixHandle, double, int> ComplexOutputs;
4951

5052
Outputs run(const Inputs& input, const Parameters& params) const;
53+
ComplexOutputs run(const ComplexInputs& input, const Parameters& params) const;
5154

52-
AlgorithmOutput run(const AlgorithmInput& input) const;
55+
AlgorithmOutput run(const AlgorithmInput& input) const override;
56+
private:
57+
template <typename In, typename Out>
58+
Out runImpl(const In& input, const Parameters& params) const;
5359
};
5460

5561
typedef boost::error_info<struct tag_eigen_computation, Eigen::ComputationInfo> EigenComputationInfo;
5662

5763
}}}}
5864

59-
#endif
65+
#endif

‎src/Core/Algorithms/Math/Tests/SolveLinearSystemWithEigenTests.cc‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -344,19 +344,19 @@ TEST(EigenSparseSolverTest, DISABLED_CanSolveBigSystem)
344344
ScopedTimer t("using algorithm object");
345345
SolveLinearSystemAlgorithm algo;
346346

347-
x = algo.run(boost::make_tuple(A, bCol), boost::make_tuple(1e-20, 4000));
348-
MatrixHandle solution = x.get<0>();
347+
x = algo.run(std::make_tuple(A, bCol), std::make_tuple(1e-20, 4000));
348+
MatrixHandle solution = std::get<0>(x);
349349

350350
ASSERT_TRUE(solution.get() != nullptr);
351-
std::cout << "error: " << x.get<1>() << std::endl;
352-
std::cout << "iterations: " << x.get<2>() << std::endl;
351+
std::cout << "error: " << std::get<1>(x) << std::endl;
352+
std::cout << "iterations: " << std::get<2>(x) << std::endl;
353353
}
354354

355355
{
356356
ScopedTimer t("comparing solutions.");
357357
auto xFileEigen = TestResources::rootDir() / "CGDarrell" / "xEigenNEW.txt";
358358
std::ofstream output(xFileEigen.string());
359-
auto solution = *x.get<0>();
359+
auto solution = *std::get<0>(x);
360360
output << std::setprecision(15) << solution << std::endl;
361361

362362
auto xFileScirun = TestResources::rootDir() / "CGDarrell" / "xScirunColumn.mat";

‎src/Core/Datatypes/DenseMatrix.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -35,6 +35,7 @@
3535
#define register
3636
#include <Eigen/Dense>
3737
#undef register
38+
#include <Core/Datatypes/share.h>
3839

3940
namespace SCIRun {
4041
namespace Core {

‎src/Core/Datatypes/Matrix.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -67,6 +67,7 @@ namespace Datatypes {
6767
class MatrixBase : public MatrixIOBase, public HasPropertyManager
6868
{
6969
public:
70+
typedef T value_type;
7071
virtual size_t nrows() const = 0;
7172
virtual size_t ncols() const = 0;
7273
size_t get_dense_size() const { return nrows() * ncols(); }

‎src/Core/Datatypes/MatrixFwd.h‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -66,10 +66,13 @@ namespace Datatypes {
6666
class DenseColumnMatrixGeneric;
6767

6868
typedef DenseColumnMatrixGeneric<double> DenseColumnMatrix;
69+
using ComplexDenseColumnMatrix = DenseColumnMatrixGeneric<complex>;
6970

7071
typedef SharedPointer<DenseColumnMatrix> DenseColumnMatrixHandle;
7172
typedef SharedPointer<const DenseColumnMatrix> DenseColumnMatrixConstHandle;
7273

74+
typedef SharedPointer<ComplexDenseColumnMatrix> ComplexDenseColumnMatrixHandle;
75+
7376
template <typename T>
7477
class SparseRowMatrixGeneric;
7578

0 commit comments

Comments
 (0)