fmatrixblock.h

Go to the documentation of this file.
00001 #ifndef __fmatrixblock_h
00002 #define __fmatrixblock_h
00003 
00004 #include  "entrymatrix.h"
00005 #include  "nodematrix.h"
00006 
00007 /*****************************************************/
00008 
00009 namespace Gascoigne
00010 {
00011 template<int N>
00012 class FMatrixBlock : public NodeMatrix<N,float>
00013 {
00014   typedef typename NodeMatrix<N,float>::iterator        iterator;
00015   typedef typename NodeMatrix<N,float>::const_iterator  const_iterator;
00016 
00017   typedef nvector<double>::iterator        viterator;
00018   typedef nvector<double>::const_iterator  const_viterator;
00019 
00020 public:
00021 
00022   int ncomp() const { return N;}
00023 
00024   void   operator *= (const FMatrixBlock<N>&);
00025   void   operator *= (double s);
00026 
00027   void   transpose();
00028   void   transpose(FMatrixBlock<N>& A);
00029   void   copy_transpose(const FMatrixBlock<N>& A);
00030 
00031   void   zero_row(int);
00032   void   uno_diag(int);
00033   float& diag(int i);
00034 
00035   void   DirichletRow (const std::vector<int>& cv);
00036   void   DirichletCol (const std::vector<int>& cv);
00037   void   DirichletDiag(const std::vector<int>& cv);
00038  
00039 
00040 
00041   void   entry     (const nmatrix<double>&);
00042   void   entry     (int i, int j, const EntryMatrix&, double s=1.);
00043   void   dual_entry(int i, int j, const EntryMatrix&, double s=1.);
00044   void   inverse ();
00045   void   vmult   (viterator) const;
00046   void   mult    (FMatrixBlock<N>&, const FMatrixBlock<N>&) const; 
00047 
00048   void   submult(const FMatrixBlock<N>& B, const FMatrixBlock<N>& C)
00049   {
00050     // this -= B*C
00051     nvector<float>::iterator p(numfixarray<N*N,float>::begin());
00052     for (char i=0; i<N; i++)
00053       {
00054         for (char j=0; j<N; j++)
00055           {
00056             nvector<float>::const_iterator pC(C.begin()+j);
00057             nvector<float>::const_iterator pB(B.begin()+i*N);
00058             nvector<float>::const_iterator qB(pB+N);
00059             //for (int k=0; k<N; k++)
00060             for (; pB!=qB; pB++)
00061               {
00062                 //value(i,j) -= B(i,k) * C(k,j);
00063                 *p -= *pB * *pC;
00064                 pC += N;
00065               }
00066             p++;
00067         }      
00068       }
00069   }
00070 
00071   void add(double s, const FMatrixBlock<N>& A)
00072     {
00073       for (int i=0; i<N; i++)
00074         {
00075           for (int j=0; j<N; j++)
00076             {
00077               NodeMatrix<N,float>::value(i,j) += s*A(i,j);
00078             }
00079         }
00080     }
00081   void adddiag(const nvector<double>& s, double l)
00082     {
00083       for (int i=0; i<N; i++)
00084         {
00085           NodeMatrix<N,float>::value(i,i) += s[i]*l;
00086         }
00087     }
00088   void add(double s, const TimePattern& TP);
00089 
00090   void cadd(double s, viterator p, const_viterator q0) const
00091     {
00092       const_iterator pm = numfixarray<N*N,float>::begin();
00093       const_viterator pend = p+N;
00094       for ( ; p!=pend; p++)
00095         {
00096           double sum = 0.;
00097           const_viterator qend(q0+N);
00098           for (const_viterator q=q0; q!=qend; q++)
00099             {
00100               sum +=  *pm++ * *q;
00101             }
00102           *p += s*sum;
00103         }
00104     }
00105   void caddtrans(double s, viterator p, const_viterator q0) const
00106     {
00107       const_iterator pm = numfixarray<N*N,float>::begin();
00108 
00109       for (int k=0; k<N; k++)
00110         {
00111           for (int h=0; h<N; h++)
00112             {
00113               //sum += M(k,h) * (*(q0+h));
00114               *p++ += *pm++ * *q0;
00115             }
00116           p -= N;
00117           q0++;
00118         }
00119     }
00120   void subtract(viterator p0, const_viterator q0) const
00121     {
00122       const_iterator pm = numfixarray<N*N,float>::begin();
00123 
00124       for (viterator p(p0); p!=p0+N; p++)
00125         {
00126           for (const_viterator q(q0); q!=q0+N; q++)
00127             {
00128               *p -= *pm++ * *q;
00129             }
00130         }
00131     };
00132   std::ostream& print(std::ostream& s) const;
00133 
00134   // Zugriff auf Inhalt ueber ganzen Vektor, damits auch ohne
00135   // Struktur geht.
00136   void vector_get(nvector<float>& v) const
00137     {
00138       v.resize(numfixarray<N*N,float>::size());
00139       for (int i=0;i<numfixarray<N*N,float>::size();++i)
00140         v[i]=NodeMatrix<N,float>::operator[](i);
00141     }
00142   void vector_set(nvector<float>& v)
00143     {
00144       assert(v.size()==this->size());
00145       for (int i=0;i<numfixarray<N*N,float>::size();++i)
00146         NodeMatrix<N,float>::operator[](i)=v[i];
00147     }
00148   void vector_add(double d, nvector<float>& v)
00149     {
00150       assert(v.size()==N*N);
00151       for (int i=0;i<NodeMatrix<N,float>::size();++i)
00152         NodeMatrix<N,float>::operator[](i)+=d*v[i];
00153     }
00154 };
00155 }
00156 
00157 #endif
00158 
00159 

Generated on Thu Sep 14 10:34:37 2006 for Gascoigne by  doxygen 1.4.7