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
|
/*=========================================================================
Program: Insight Segmentation & Registration Toolkit
Module: itkFEMLinearSystemWrapperVNL.h
Language: C++
Date: $Date$
Version: $Revision$
Copyright (c) Insight Software Consortium. All rights reserved.
See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#ifndef __itkFEMLinearSystemWrapperVNL_h
#define __itkFEMLinearSystemWrapperVNL_h
#include "itkFEMLinearSystemWrapper.h"
#include "vnl/vnl_sparse_matrix.h"
#include "vnl/vnl_vector.h"
#include <vnl/vnl_sparse_matrix_linear_system.h>
#include <vnl/algo/vnl_lsqr.h>
#include <vector>
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
*/
class LinearSystemWrapperVNL : public LinearSystemWrapper
{
public:
/* values stored in matrices & vectors */
typedef LinearSystemWrapper::Float Float;
/* superclass */
typedef LinearSystemWrapper SuperClass;
/* matrix typedef */
typedef vnl_sparse_matrix<Float> MatrixRepresentation;
/* matrix holder typedef */
typedef std::vector< MatrixRepresentation* > MatrixHolder;
/* constructor & destructor */
LinearSystemWrapperVNL() : LinearSystemWrapper(), m_Matrices(0), m_Vectors(0), m_Solutions(0) {}
virtual ~LinearSystemWrapperVNL();
/* memory management routines */
virtual void InitializeMatrix(unsigned int matrixIndex);
virtual bool IsMatrixInitialized(unsigned int matrixIndex);
virtual void DestroyMatrix(unsigned int matrixIndex);
virtual void InitializeVector(unsigned int vectorIndex);
virtual bool IsVectorInitialized(unsigned int vectorIndex);
virtual void DestroyVector(unsigned int vectorIndex);
virtual void InitializeSolution(unsigned int solutionIndex);
virtual bool IsSolutionInitialized(unsigned int solutionIndex);
virtual void DestroySolution(unsigned int solutionIndex);
virtual void SetMaximumNonZeroValuesInMatrix(unsigned int, unsigned int) {}
/* assembly & solving routines */
virtual Float GetMatrixValue(unsigned int i, unsigned int j, unsigned int matrixIndex) const { return (*((*m_Matrices)[matrixIndex]))(i,j); }
virtual void SetMatrixValue(unsigned int i, unsigned int j, Float value, unsigned int matrixIndex) { (*((*m_Matrices)[matrixIndex]))(i,j) = value; }
virtual void AddMatrixValue(unsigned int i, unsigned int j, Float value, unsigned int matrixIndex) { (*((*m_Matrices)[matrixIndex]))(i,j) += value; }
virtual Float GetVectorValue(unsigned int i, unsigned int vectorIndex) const { return (* ( (*m_Vectors)[vectorIndex] ) )[i]; }
virtual void SetVectorValue(unsigned int i, Float value, unsigned int vectorIndex) { (*((*m_Vectors)[vectorIndex]))(i) = value; }
virtual void AddVectorValue(unsigned int i, Float value, unsigned int vectorIndex) { (*((*m_Vectors)[vectorIndex]))(i) += value; }
virtual Float GetSolutionValue(unsigned int i, unsigned int solutionIndex) const;
virtual void SetSolutionValue(unsigned int i, Float value, unsigned int solutionIndex) { (*((*m_Solutions)[solutionIndex]))(i) = value; }
virtual void AddSolutionValue(unsigned int i, Float value, unsigned int solutionIndex) { (*((*m_Solutions)[solutionIndex]))(i) += value; }
virtual void Solve(void);
/* matrix & vector manipulation routines */
virtual void ScaleMatrix(Float scale, unsigned int matrixIndex);
virtual void SwapMatrices(unsigned int matrixIndex1, unsigned int matrixIndex2);
virtual void SwapVectors(unsigned int vectorIndex1, unsigned int vectorIndex2);
virtual void SwapSolutions(unsigned int solutionIndex1, unsigned int solutionIndex2);
virtual void CopySolution2Vector(unsigned solutionIndex, unsigned int vectorIndex);
virtual void CopyVector2Solution(unsigned int vectorIndex, unsigned int solutionIndex);
virtual void MultiplyMatrixMatrix(unsigned int resultMatrixIndex, unsigned int leftMatrixIndex, unsigned int rightMatrixIndex);
virtual void MultiplyMatrixVector(unsigned int resultVectorIndex, unsigned int matrixIndex, unsigned int vectorIndex);
private:
/** vector of pointers to VNL sparse matrices */
//std::vector< vnl_sparse_matrix<Float>* > *m_Matrices;
MatrixHolder *m_Matrices;
/** vector of pointers to VNL vectors */
std::vector< vnl_vector<Float>* > *m_Vectors;
/** vector of pointers to VNL vectors */
std::vector< vnl_vector<Float>* > *m_Solutions;
};
}} // end namespace itk::fem
#endif
|
