Memosa-FVM  0.2
Public Types | Public Member Functions | Private Attributes | List of all members
SquareMatrix< T, N > Class Template Reference

#include <MatrixOperation.h>

Inheritance diagram for SquareMatrix< T, N >:
Inheritance graph
[legend]
Collaboration diagram for SquareMatrix< T, N >:
Collaboration graph
[legend]

Public Types

enum  { NSQR = N*N }
 
typedef Array< T > TArray
 
typedef shared_ptr< TArray > TArrPtr
 
typedef Array< int > IntArray
 
typedef NumTypeTraits< T >
::T_Scalar 		T_Scalar
 
- Public Types inherited from MatrixJML< T >
typedef Array< T > TArray
 

Public Member Functions

 SquareMatrix ()
 
 ~SquareMatrix ()
 
 SquareMatrix (const SquareMatrix &o)
 
 SquareMatrix (const T &s)
 
T & operator() (int i, int j)
 
SquareMatrix & operator= (const T &s)
 
SquareMatrix & operator= (const SquareMatrix &o)
 
void print ()
 
 SquareMatrix (const int N)
 
T & getElement (const int i, const int j)
 
T & operator() (const int i, const int j)
 
void zero ()
 
void Solve (TArray &bVec)
 
void makeCopy (SquareMatrix< T > &o)
 
void printMatrix ()
 
T getTraceAbs ()
 
void multiply (const TArray &x, TArray &b)
 
void testSolve ()
 
- Public Member Functions inherited from MatrixJML< T >
 MatrixJML ()
 
virtual ~MatrixJML ()
 
 MatrixJML ()
 
virtual ~MatrixJML ()
 

Private Attributes

T _data [NSQR]
 
const int _order
 
const int _elements
 
bool _sorted
 
IntArray _pivotRows
 
TArray _maxVals
 
TArray _values
 

Detailed Description

template<class T, int N>
class SquareMatrix< T, N >

Definition at line 9 of file MatrixOperation.h.

Member Typedef Documentation

template<class T, int N>
typedef Array<int> SquareMatrix< T, N >::IntArray

Definition at line 18 of file SquareMatrix.h.

template<class T, int N>
typedef NumTypeTraits<T>::T_Scalar SquareMatrix< T, N >::T_Scalar

Definition at line 19 of file SquareMatrix.h.

template<class T, int N>
typedef Array<T> SquareMatrix< T, N >::TArray

Definition at line 16 of file SquareMatrix.h.

template<class T, int N>
typedef shared_ptr<TArray> SquareMatrix< T, N >::TArrPtr

Definition at line 17 of file SquareMatrix.h.

Member Enumeration Documentation

template<class T, int N>
anonymous enum
Enumerator
NSQR 

Definition at line 12 of file MatrixOperation.h.

12 { NSQR = N*N };

Constructor & Destructor Documentation

template<class T, int N>
SquareMatrix< T, N >::SquareMatrix ( )
inline

Definition at line 14 of file MatrixOperation.h.

14 {}
template<class T, int N>
SquareMatrix< T, N >::~SquareMatrix ( )
inline

Definition at line 15 of file MatrixOperation.h.

15 {}
template<class T, int N>
SquareMatrix< T, N >::SquareMatrix ( const SquareMatrix< T, N > &  o)
inline

Definition at line 17 of file MatrixOperation.h.

References SquareMatrix< T, N >::_data, and SquareMatrix< T, N >::NSQR.

17  {
18  for ( int i=0; i<NSQR; i++){
19  _data[i] = o._data[i];
20  }
21  }
template<class T, int N>
SquareMatrix< T, N >::SquareMatrix ( const T &  s)
inline

Definition at line 23 of file MatrixOperation.h.

23  {
24  for ( int i=0; i<N; i++){
25  for ( int j=0; j<N; j++){
26  if (i==j)
27  (*this)(i,j) = s;
28  else
29  (*this)(i,j) = T(0);
30  }
31  }
32  }
template<class T, int N>
SquareMatrix< T, N >::SquareMatrix ( const int  N)
inline

Definition at line 21 of file SquareMatrix.h.

References SquareMatrix< T, N >::_values, and Array< T >::zero().

21  :
22  _order(N),
23  _elements(N*N),
24  _sorted(false),
25  _pivotRows(N),
26  _maxVals(N),
27  _values(_elements)
28  {_values.zero();}
Array::zero
virtual void zero()
Definition: Array.h:281
SquareMatrix::_maxVals
TArray _maxVals
Definition: SquareMatrix.h:249
SquareMatrix::_order
const int _order
Definition: SquareMatrix.h:245
SquareMatrix::_pivotRows
IntArray _pivotRows
Definition: SquareMatrix.h:248
SquareMatrix::_elements
const int _elements
Definition: SquareMatrix.h:246
SquareMatrix::_values
TArray _values
Definition: SquareMatrix.h:250

Member Function Documentation

template<class T, int N>
T& SquareMatrix< T, N >::getElement ( const int  i,
const int  j 
)
inlinevirtual

Implements MatrixJML< T >.

Definition at line 30 of file SquareMatrix.h.

References SquareMatrix< T, N >::_order, and SquareMatrix< T, N >::_values.

30 {return _values[(i-1)*_order+j-1];}
SquareMatrix::_order
const int _order
Definition: SquareMatrix.h:245
SquareMatrix::_values
TArray _values
Definition: SquareMatrix.h:250
template<class T, int N>
T SquareMatrix< T, N >::getTraceAbs ( )
inline

Definition at line 186 of file SquareMatrix.h.

References SquareMatrix< T, N >::_order, and fabs().

Referenced by COMETDiscretizer< T >::findResidCoarse(), COMETDiscretizer< T >::findResidFine(), and COMETDiscretizer< T >::findResidFull().

187  {
188  T trace=0.;
189  for(int i=1;i<_order+1;i++)
190  trace+=fabs((*this)(i,i));
191  return trace;
192  }
SquareMatrix::_order
const int _order
Definition: SquareMatrix.h:245
fabs
Tangent fabs(const Tangent &a)
Definition: Tangent.h:312
template<class T, int N>
void SquareMatrix< T, N >::makeCopy ( SquareMatrix< T > &  o)
inline

Definition at line 164 of file SquareMatrix.h.

References SquareMatrix< T, N >::_maxVals, SquareMatrix< T, N >::_order, SquareMatrix< T, N >::_pivotRows, SquareMatrix< T, N >::_sorted, and SquareMatrix< T, N >::_values.

165  {
166  if(o._order!=this->_order)
167  throw CException("Cannot copy matrices of different sizes!");
168 
169  o._sorted=this->_sorted;
170  o._pivotRows=this->_pivotRows;
171  o._maxVals=this->_maxVals;
172  o._values=this->_values;
173  }
SquareMatrix::_maxVals
TArray _maxVals
Definition: SquareMatrix.h:249
SquareMatrix::_order
const int _order
Definition: SquareMatrix.h:245
SquareMatrix::_pivotRows
IntArray _pivotRows
Definition: SquareMatrix.h:248
SquareMatrix::_values
TArray _values
Definition: SquareMatrix.h:250
template<class T, int N>
void SquareMatrix< T, N >::multiply ( const TArray &  x,
TArray &  b 
)
inlinevirtual

Implements MatrixJML< T >.

Definition at line 194 of file SquareMatrix.h.

References SquareMatrix< T, N >::_order, Array< T >::getLength(), and Array< T >::zero().

Referenced by COMETDiscretizer< T >::findResidCoarse(), COMETDiscretizer< T >::findResidFine(), COMETDiscretizer< T >::findResidFull(), and SquareMatrix< T, N >::testSolve().

195  {
196  int lenx=x.getLength();
197  int lenb=b.getLength();
198  if(lenx==_order && lenb==_order)
199  {
200  b.zero();
201  for(int i=1;i<_order+1;i++)
202  for(int j=1;j<_order+1;j++)
203  b[i-1]+=(*this)(i,j)*x[j-1];
204  }
205  else
206  throw CException("Array length does not match matrix order!");
207  }
SquareMatrix::_order
const int _order
Definition: SquareMatrix.h:245
template<class T, int N>
T& SquareMatrix< T, N >::operator() ( const int  i,
const int  j 
)
inline

Definition at line 31 of file SquareMatrix.h.

References SquareMatrix< T, N >::_order, and SquareMatrix< T, N >::_values.

31 {return _values[(i-1)*_order+j-1];}
SquareMatrix::_order
const int _order
Definition: SquareMatrix.h:245
SquareMatrix::_values
TArray _values
Definition: SquareMatrix.h:250
template<class T, int N>
T& SquareMatrix< T, N >::operator() ( int  i,
int  j 
)
inline

Definition at line 34 of file MatrixOperation.h.

References SquareMatrix< T, N >::_data.

34 {return _data[i*N+j];}
template<class T, int N>
SquareMatrix& SquareMatrix< T, N >::operator= ( const T &  s)
inline

Definition at line 36 of file MatrixOperation.h.

37  {
38  for(int i=0; i<N; i++)
39  for(int j=0; j<N; j++)
40  (*this)(i,j) = (i==j) ? s : T(0.);
41  return *this;
42  }
template<class T, int N>
SquareMatrix& SquareMatrix< T, N >::operator= ( const SquareMatrix< T, N > &  o)
inline

Definition at line 44 of file MatrixOperation.h.

References SquareMatrix< T, N >::_data, and SquareMatrix< T, N >::NSQR.

45  {
46  for(int i=0;i<NSQR;i++)
47  _data[i] = o._data[i];
48  return *this;
49  }
template<class T, int N>
void SquareMatrix< T, N >::print ( )
inline

Definition at line 51 of file MatrixOperation.h.

51  {
52  for(int i=0; i<N; i++){
53  for (int j=0; j<N; j++)
54  cout << (*this)(i,j) << " ";
55  cout << endl;
56  }
57  }
template<class T, int N>
void SquareMatrix< T, N >::printMatrix ( )
inline

Definition at line 175 of file SquareMatrix.h.

References SquareMatrix< T, N >::_order.

176  {
177  for(int i=1;i<_order+1;i++)
178  {
179  for(int j=1;j<_order+1;j++)
180  cout<<(*this)(i,j)<<" ";
181  cout<<endl;
182  }
183  cout<<endl;
184  }
SquareMatrix::_order
const int _order
Definition: SquareMatrix.h:245
template<class T, int N>
void SquareMatrix< T, N >::Solve ( TArray &  bVec)
inline

Definition at line 34 of file SquareMatrix.h.

References SquareMatrix< T, N >::_maxVals, SquareMatrix< T, N >::_order, SquareMatrix< T, N >::_pivotRows, SquareMatrix< T, N >::_sorted, and fabs().

Referenced by COMETDiscretizer< T >::COMETSolveCoarse(), COMETDiscretizer< T >::COMETSolveFine(), COMETDiscretizer< T >::COMETSolveFull(), and SquareMatrix< T, N >::testSolve().

35  {//Gaussian Elimination w/ scaled partial pivoting
36  //replaces bVec with the solution vector.
37 
38  SquareMatrix<T> LU(_order);
39  (*this).makeCopy(LU);
40  TArray bCpy(_order);
41  IntArray l(_order);
42  TArray s(_order);
43 
44  //find max values in each row if not done yet
45  if(!_sorted)
46  {
47  for(int i=1;i<_order+1;i++)
48  {
49  l[i-1]=i;
50  s[i-1]=fabs((*this)(i,1));
51  for(int j=2;j<_order+1;j++)
52  {
53  if(s[i-1]<fabs((*this)(i,j)))
54  s[i-1]=fabs((*this)(i,j));
55  }
56  if(s[i-1]==0)
57  {
58  cout<<"Row: "<<i<<endl;
59  throw CException("Matrix has row of zeros!");
60  }
61  }
62  _maxVals=s;
63  //cout<<"Max vals:"<<endl;
64  //for(int i=0;i<_order;i++)
65  //cout<<s[i]<<endl;
66  }
67 
68  //Forward sweep
69  if(!_sorted)
70  {
71  for(int i=1;i<_order;i++)
72  {
73  T rmax=fabs((*this)(l[i-1],i)/s[l[i-1]-1]);
74  // cout<<"rmax: "<<rmax<<endl;
75  int newMax=-1;
76  for(int j=i+1;j<_order+1;j++)
77  {
78  T r=fabs((*this)(l[j-1],i)/s[l[j-1]-1]);
79  if(r>rmax)
80  {
81  //cout<<"switching,i,j,r:"<<i<<","<<j<<","<<r<<endl;
82  //cout<<"factors "<<(*this)(l[j-1],i)<<" "<<s[l[j-1]-1]<<endl;
83  //cout<<"mapped j "<<l[j-1]<<endl;
84  rmax=r;
85  newMax=j;
86  }
87  }
88 
89  if(newMax!=-1)
90  {
91  int temp=l[i-1];
92  l[i-1]=newMax;
93  l[newMax-1]=temp;
94  }
95 
96  for(int j=i+1;j<_order+1;j++)
97  {
98  T factor=LU(l[j-1],i)/LU(l[i-1],i);
99  LU(l[j-1],i)=factor;
100  bVec[l[j-1]-1]-=factor*bVec[l[i-1]-1];
101  for(int k=i+1;k<_order+1;k++)
102  {
103  LU(l[j-1],k)-=LU(l[i-1],k)*factor;
104  T test=LU(l[j-1],k);
105  if(isnan(test)||isinf(test))
106  {
107  cout<<"Denom: "<<LU(l[i-1],i)<<endl;
108  cout<<"Num: "<<LU(l[j-1],i)<<endl;
109  cout<<"first: "<<LU(l[j-1],k)<<endl;
110  cout<<"second: "<<LU(l[i-1],k)<<endl;
111  cout<<"factor: "<<factor<<endl;
112  cout<<"test: "<<test<<endl;
113  throw CException("test is nan");
114  }
115  }
116  }
117  }
118  _pivotRows=l;
119  _sorted=true;
120  //cout<<"Order"<<endl;
121  //for(int i=0;i<_order;i++)
122  // cout<<_pivotRows[i]<<endl;
123  //cout<<endl;
124  }
125  else
126  {
127  for(int i=1;i<_order;i++)
128  {
129  for(int j=i+1;j<_order+1;j++)
130  {
131  T factor=LU(_pivotRows[j-1],i)/LU(_pivotRows[i-1],i);
132  LU(_pivotRows[j-1],i)=factor;
133  bVec[_pivotRows[j-1]-1]-=factor*bVec[_pivotRows[i-1]-1];
134  for(int k=i+1;k<_order+1;k++)
135  {
136  LU(_pivotRows[j-1],k)=LU(_pivotRows[j-1],k)
137  -LU(_pivotRows[i-1],k)*factor;
138  }
139  }
140  }
141  }
142 
143  //back solve
144  bVec[_pivotRows[_order-1]-1]=
145  bVec[_pivotRows[_order-1]-1]/LU(_pivotRows[_order-1],_order);
146 
147  T sum=0.;
148  for(int i=_order-1;i>0;i--)
149  {
150  sum=0.;
151  for(int j=i+1;j<_order+1;j++)
152  sum-=LU(_pivotRows[i-1],j)*bVec[_pivotRows[j-1]-1];
153  bVec[_pivotRows[i-1]-1]+=sum;
154  bVec[_pivotRows[i-1]-1]=bVec[_pivotRows[i-1]-1]/LU(_pivotRows[i-1],i);
155  }
156 
157  //reorder
158  bCpy=bVec;
159  for(int i=0;i<_order;i++)
160  bVec[i]=bCpy[_pivotRows[i]-1];
161 
162  }
SquareMatrix::IntArray
Array< int > IntArray
Definition: SquareMatrix.h:18
SquareMatrix::_maxVals
TArray _maxVals
Definition: SquareMatrix.h:249
SquareMatrix::_order
const int _order
Definition: SquareMatrix.h:245
SquareMatrix::_pivotRows
IntArray _pivotRows
Definition: SquareMatrix.h:248
SquareMatrix::TArray
Array< T > TArray
Definition: SquareMatrix.h:16
fabs
Tangent fabs(const Tangent &a)
Definition: Tangent.h:312
template<class T, int N>
void SquareMatrix< T, N >::testSolve ( )
inline

Definition at line 209 of file SquareMatrix.h.

References SquareMatrix< T, N >::_order, SquareMatrix< T, N >::multiply(), and SquareMatrix< T, N >::Solve().

210  {
211  TArray x(_order);
212  TArray b(_order);
213 
214  cout<<"Correct Solution:"<<endl;
215  for(int i=0;i<_order;i++)
216  {
217  x[i]=rand()/1.e9;
218  cout<<"x["<<i<<"]="<<x[i]<<endl;
219  }
220 
221  multiply(x,b);
222 
223  cout<<"Before"<<endl;
224  for(int i=0;i<_order;i++)
225  cout<<"b["<<i<<"]="<<b[i]<<endl;
226 
227  Solve(b);
228 
229  cout<<"After"<<endl;
230  for(int i=0;i<_order;i++)
231  cout<<"b["<<i<<"]="<<b[i]<<endl;
232  cout<<endl;
233 
234  cout<<"Error"<<endl;
235  for(int i=0;i<_order;i++)
236  cout<<"b["<<i<<"]="<<b[i]-x[i]<<endl;
237  cout<<endl;
238 
239 
240  throw CException("finished with test");
241 
242  }
SquareMatrix::multiply
void multiply(const TArray &x, TArray &b)
Definition: SquareMatrix.h:194
SquareMatrix::_order
const int _order
Definition: SquareMatrix.h:245
SquareMatrix::TArray
Array< T > TArray
Definition: SquareMatrix.h:16
SquareMatrix::Solve
void Solve(TArray &bVec)
Definition: SquareMatrix.h:34
template<class T, int N>
void SquareMatrix< T, N >::zero ( )
inlinevirtual

Implements MatrixJML< T >.

Definition at line 32 of file SquareMatrix.h.

References SquareMatrix< T, N >::_values, and Array< T >::zero().

Referenced by COMETDiscretizer< T >::COMETSolveCoarse(), COMETDiscretizer< T >::COMETSolveFine(), COMETDiscretizer< T >::COMETSolveFull(), COMETDiscretizer< T >::findResidCoarse(), COMETDiscretizer< T >::findResidFine(), and COMETDiscretizer< T >::findResidFull().

32 {_values.zero();}
Array::zero
virtual void zero()
Definition: Array.h:281
SquareMatrix::_values
TArray _values
Definition: SquareMatrix.h:250

Member Data Documentation

template<class T, int N>
T SquareMatrix< T, N >::_data[NSQR]
private

Definition at line 62 of file MatrixOperation.h.

Referenced by SquareMatrix< T, N >::operator()(), SquareMatrix< T, N >::operator=(), and SquareMatrix< T, N >::SquareMatrix().

template<class T, int N>
const int SquareMatrix< T, N >::_elements
private

Definition at line 246 of file SquareMatrix.h.

template<class T, int N>
TArray SquareMatrix< T, N >::_maxVals
private

Definition at line 249 of file SquareMatrix.h.

Referenced by SquareMatrix< T, N >::makeCopy(), and SquareMatrix< T, N >::Solve().

template<class T, int N>
const int SquareMatrix< T, N >::_order
private

Definition at line 245 of file SquareMatrix.h.

Referenced by SquareMatrix< T, N >::getElement(), SquareMatrix< T, N >::getTraceAbs(), SquareMatrix< T, N >::makeCopy(), SquareMatrix< T, N >::multiply(), SquareMatrix< T, N >::operator()(), SquareMatrix< T, N >::printMatrix(), SquareMatrix< T, N >::Solve(), and SquareMatrix< T, N >::testSolve().

template<class T, int N>
IntArray SquareMatrix< T, N >::_pivotRows
private

Definition at line 248 of file SquareMatrix.h.

Referenced by SquareMatrix< T, N >::makeCopy(), and SquareMatrix< T, N >::Solve().

template<class T, int N>
bool SquareMatrix< T, N >::_sorted
private

Definition at line 247 of file SquareMatrix.h.

Referenced by SquareMatrix< T, N >::makeCopy(), and SquareMatrix< T, N >::Solve().

template<class T, int N>
TArray SquareMatrix< T, N >::_values
private

Definition at line 250 of file SquareMatrix.h.

Referenced by SquareMatrix< T, N >::getElement(), SquareMatrix< T, N >::makeCopy(), SquareMatrix< T, N >::operator()(), SquareMatrix< T, N >::SquareMatrix(), and SquareMatrix< T, N >::zero().

The documentation for this class was generated from the following files: