1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
|
/*=========================================================================
*
* Copyright NumFOCUS
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* https://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
#ifndef itkFEMLinearSystemWrapperVNL_h
#define itkFEMLinearSystemWrapperVNL_h
#include "itkFEMLinearSystemWrapper.h"
#include "vnl/vnl_sparse_matrix.h"
#include "vnl/vnl_vector.h"
#include <vector>
#include "ITKFEMExport.h"
namespace itk
{
namespace fem
{
/**
* \class LinearSystemWrapperVNL
* \brief LinearSystemWrapper class that uses VNL numeric library functions
* to define a sparse linear system of equations.
* \sa LinearSystemWrapper
* \ingroup ITKFEM
*/
class ITKFEM_EXPORT LinearSystemWrapperVNL : public LinearSystemWrapper
{
public:
/* values stored in matrices & vectors */
using Float = LinearSystemWrapper::Float;
/* superclass */
using Superclass = LinearSystemWrapper;
/* matrix type alias */
using MatrixRepresentation = vnl_sparse_matrix<Float>;
/* matrix holder type alias */
using MatrixHolder = std::vector<MatrixRepresentation *>;
/* constructor & destructor */
LinearSystemWrapperVNL()
: LinearSystemWrapper()
{}
~LinearSystemWrapperVNL() override;
/* memory management routines */
void
InitializeMatrix(unsigned int matrixIndex) override;
bool
IsMatrixInitialized(unsigned int matrixIndex) override;
void
DestroyMatrix(unsigned int matrixIndex) override;
void
InitializeVector(unsigned int vectorIndex) override;
bool
IsVectorInitialized(unsigned int vectorIndex) override;
void
DestroyVector(unsigned int vectorIndex) override;
void
InitializeSolution(unsigned int solutionIndex) override;
bool
IsSolutionInitialized(unsigned int solutionIndex) override;
void
DestroySolution(unsigned int solutionIndex) override;
virtual void
SetMaximumNonZeroValuesInMatrix(unsigned int, unsigned int)
{}
/* assembly & solving routines */
Float
GetMatrixValue(unsigned int i, unsigned int j, unsigned int matrixIndex) const override
{
return (*((*m_Matrices)[matrixIndex]))(i, j);
}
void
SetMatrixValue(unsigned int i, unsigned int j, Float value, unsigned int matrixIndex) override
{
(*((*m_Matrices)[matrixIndex]))(i, j) = value;
}
void
AddMatrixValue(unsigned int i, unsigned int j, Float value, unsigned int matrixIndex) override
{
(*((*m_Matrices)[matrixIndex]))(i, j) += value;
}
Float
GetVectorValue(unsigned int i, unsigned int vectorIndex) const override
{
return (*((*m_Vectors)[vectorIndex]))[i];
}
void
SetVectorValue(unsigned int i, Float value, unsigned int vectorIndex) override
{
(*((*m_Vectors)[vectorIndex]))(i) = value;
}
void
AddVectorValue(unsigned int i, Float value, unsigned int vectorIndex) override
{
(*((*m_Vectors)[vectorIndex]))(i) += value;
}
Float
GetSolutionValue(unsigned int i, unsigned int solutionIndex) const override;
void
SetSolutionValue(unsigned int i, Float value, unsigned int solutionIndex) override
{
(*((*m_Solutions)[solutionIndex]))(i) = value;
}
void
AddSolutionValue(unsigned int i, Float value, unsigned int solutionIndex) override
{
(*((*m_Solutions)[solutionIndex]))(i) += value;
}
void
Solve() override;
/* matrix & vector manipulation routines */
void
ScaleMatrix(Float scale, unsigned int matrixIndex) override;
void
SwapMatrices(unsigned int matrixIndex1, unsigned int matrixIndex2) override;
void
SwapVectors(unsigned int vectorIndex1, unsigned int vectorIndex2) override;
void
SwapSolutions(unsigned int solutionIndex1, unsigned int solutionIndex2) override;
void
CopySolution2Vector(unsigned int solutionIndex, unsigned int vectorIndex) override;
void
CopyVector2Solution(unsigned int vectorIndex, unsigned int solutionIndex) override;
void
MultiplyMatrixMatrix(unsigned int resultMatrixIndex,
unsigned int leftMatrixIndex,
unsigned int rightMatrixIndex) override;
void
MultiplyMatrixVector(unsigned int resultVectorIndex, unsigned int matrixIndex, unsigned int vectorIndex) override;
private:
/** vector of pointers to VNL sparse matrices */
// std::vector< vnl_sparse_matrix<Float>* > *m_Matrices;
MatrixHolder * m_Matrices{ nullptr };
/** vector of pointers to VNL vectors */
std::vector<vnl_vector<Float> *> * m_Vectors{ nullptr };
/** vector of pointers to VNL vectors */
std::vector<vnl_vector<Float> *> * m_Solutions{ nullptr };
};
} // end namespace fem
} // end namespace itk
#endif // itkFEMLinearSystemWrapperVNL_h
|
