NGSolve
4.9
|
00001 #ifndef FILE_SYMMETRICMATRIX 00002 #define FILE_SYMMETRICMATRIX 00003 00004 /****************************************************************************/ 00005 /* File: symmetricmatrix.hpp */ 00006 /* Author: Joachim Schoeberl */ 00007 /* Date: 12. Dec. 2004 */ 00008 /****************************************************************************/ 00009 00010 00011 namespace ngbla 00012 { 00013 00014 00015 00017 template <class T> 00018 class FlatSymmetricMatrix 00019 { 00020 protected: 00022 int n; 00024 T *data; 00025 public: 00026 typedef typename mat_traits<T>::TV_COL TV; 00027 00029 FlatSymmetricMatrix () { ; } 00030 00032 FlatSymmetricMatrix (int an, T * adata) 00033 : n(an), data(adata) 00034 { ; } 00035 00037 void AssignMemory (int an, T * mem) throw() 00038 { 00039 n = an; 00040 data = mem; 00041 } 00042 00043 00044 00046 void Mult (const FlatVector<TV> & x, FlatVector<TV> & y) const 00047 { 00048 for (int i = 0; i < n; i++) 00049 { 00050 TV xi = x(i); 00051 TV sum = (*this)(i,i) * xi; 00052 for (int j = 0; j < i; j++) 00053 { 00054 sum += (*this)(i,j) * x(j); 00055 y(j) += Trans((*this)(i,j)) * xi; 00056 } 00057 y(i) = sum; 00058 } 00059 } 00060 00062 ostream & Print (ostream & ost) const 00063 { 00064 for (int i = 0; i < n; i++) 00065 { 00066 for (int j = 0; j < n; j++) 00067 if (j <= i) 00068 ost << setw(8) << (*this)(i,j) << " "; 00069 else 00070 ost << setw(8) << "sym" << " "; 00071 ost << endl; 00072 } 00073 return *this; 00074 } 00075 00077 int Height() const { return n; } 00078 00079 00081 const T & operator() (int i, int j) const 00082 { return data[ i * (i+1) / 2 + j ]; } 00083 00085 T & operator() (int i, int j) 00086 { return data[ i * (i+1) / 2 + j ]; } 00087 00088 00090 FlatSymmetricMatrix & operator= (const T & val) 00091 { 00092 int nel = n * (n+1) / 2; 00093 for (int i = 0; i < nel; i++) 00094 data[i] = val; 00095 return *this; 00096 } 00097 }; 00098 00099 00100 template<typename T> 00101 inline std::ostream & operator<< (std::ostream & s, const FlatSymmetricMatrix<T> & m) 00102 { 00103 m.Print (s); 00104 return s; 00105 } 00106 00107 00108 00110 template <class T> 00111 class SymmetricMatrix : public FlatSymmetricMatrix<T> 00112 { 00113 public: 00114 typedef typename mat_traits<T>::TV_COL TV; 00115 00117 SymmetricMatrix (int an) 00118 : FlatSymmetricMatrix<T> (an, new T[an*(an+1)/2]) 00119 { ; } 00120 00122 ~SymmetricMatrix () 00123 { delete [] this->data; } 00124 00126 SymmetricMatrix & operator= (const T & val) 00127 { 00128 int nel = this->n * (this->n+1) / 2; 00129 for (int i = 0; i < nel; i++) 00130 this->data[i] = val; 00131 return *this; 00132 } 00133 }; 00134 00135 00136 00137 00138 00139 00140 00141 } 00142 00143 00144 #endif