NGSolve  4.9
basiclinalg/symmetricmatrix.hpp
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