nodematrix.h

Go to the documentation of this file.
00001 #ifndef __nodematrix_h
00002 #define __nodematrix_h
00003 
00004 #include "numfixarray.h"
00005 
00006 
00007 /**************************************************/
00008 
00009 namespace Gascoigne
00010 {
00011 template<int N, class T>
00012 class NodeMatrix : public numfixarray<N*N,T>
00013 {
00014   
00015 public:
00016   
00017   NodeMatrix<N,T>()           : numfixarray<N*N,T>() {}
00018   NodeMatrix<N,T>(const T& A) : numfixarray<N*N,T>(A) {}
00019 
00020   T&       operator()(int i,int j)       { return (*this)[j+i*N];}
00021   const T& operator()(int i,int j) const { return (*this)[j+i*N];}
00022 
00023   T&       value(int i,int j)       { return (*this)[j+i*N];}
00024   const T& value(int i,int j) const { return (*this)[j+i*N];}
00025   
00026   void reserve(int) const {};
00027   void resize(int) const {};
00028 
00029   int m() const { return N;}
00030 
00031   friend std::ostream& operator<<(std::ostream &s, const NodeMatrix<N,T>& A)
00032     {
00033       for(int c=0;c<N*N;c++)  s << A[c] << " ";
00034       s << "\t";
00035       return s;
00036     }
00037 
00038   void identity()
00039     {
00040     numfixarray<N*N,T>::zero();
00041       for(int d=0;d<N;d++)
00042         {
00043           value(d,d) = 1.;
00044         }
00045     }
00046   void zero_component(int c)
00047     {
00048       for(int d=0;d<N;d++)
00049         {
00050           value(c,d) = 0.;
00051         }
00052     }
00053 
00054 /**************************************************/
00055 
00056   void addmult(double k, const NodeMatrix<N,T>& A, const NodeMatrix<N,T>& B)  
00057     // this += k*A*B
00058     {
00059       for(int i=0;i<N;i++)
00060         {
00061           for(int j=0;j<N;j++)
00062             {
00063               for(int k=0;k<N;k++)
00064                 {
00065                   value(i,j) += k*A.value(i,k) * B.value(k,j); 
00066                 }
00067             }
00068         }
00069     }
00070 
00071 /**************************************************/
00072 
00073   void inverse() { inverse(*this);}
00074 
00075 /**************************************************/
00076 
00077   void inverse(const NodeMatrix<N,T>& A)
00078     {
00079       /*for(int i=0;i<N;i++)
00080         {
00081           value(i,i) = 1./A.value(i,i);
00082         }
00083       return;*/
00084 
00085       *this = A;
00086 
00087       
00088       // LU decomposition
00089 
00090       for(int i=1;i<N;i++)
00091         {
00092           for(int k=0;k<i;k++)
00093             {
00094               value(i,k) /= value(k,k);
00095               for(int j=k+1;j<N;j++)
00096                 {
00097                   value(i,j) -= value(i,k)*value(k,j);
00098                 }
00099             }
00100         }
00101 
00102       // Inverse von L
00103 
00104       for(int ncol=0;ncol<N-1;ncol++)
00105         {
00106           for(int i=ncol+1;i<N;i++)
00107             {
00108               value(i,ncol) = -value(i,ncol);
00109               for(int k=ncol+1;k<i;k++)
00110                 {
00111                   value(i,ncol) -= value(i,k)*value(k,ncol);
00112                 }
00113             }
00114         }
00115       
00116 
00117       // Inverse von U
00118 
00119 
00120       for(int nlin=0;nlin<N;nlin++)
00121         {
00122           for(int j=nlin+1;j<N;j++)
00123             {
00124               value(nlin,j) /= -value(nlin,nlin);
00125               for(int k=nlin+1;k<j;k++)
00126                 {
00127                   value(nlin,j) -= value(nlin,k)*value(k,j);
00128                 }
00129               value(nlin,j) /= value(j,j);
00130             }
00131           value(nlin,nlin) = 1./value(nlin,nlin);
00132         }
00133       
00134 
00135       // Inverse von A
00136 
00137       for(int ncol=0;ncol<N;ncol++)
00138         {
00139           for(int i=0;i<ncol+1;i++)
00140             {
00141               for(int k=ncol+1;k<N;k++)
00142                 {
00143                   value(i,ncol) += value(i,k)*value(k,ncol);
00144                 }
00145             }
00146           for(int i=ncol+1;i<N;i++)
00147             {
00148               value(i,ncol) *= value(i,i);
00149               for(int k=i+1;k<N;k++)
00150                 {
00151                   value(i,ncol) += value(i,k)*value(k,ncol);
00152                 }
00153             }
00154         }
00155       
00156     }
00157 
00158   void gauss_jordan()
00159     {
00160       std::vector<int> p(N);
00161 
00162       int i,j,k,r;
00163       double max, hr;
00164       
00165       for (i=0;i<N;i++) p[i] = i;
00166       
00167       for (j=0;j<N;j++)
00168         {
00169           max = fabs(value(j,j));
00170           r = j;
00171           for (i=j+1;i<N;i++)
00172             {
00173               if (fabs(value(i,j)) > max)
00174                 {
00175                   max = fabs(value(i,j));
00176                   r = i;
00177                 }
00178             }
00179           if (r>j)
00180             {
00181               for (k=0;k<N;k++)
00182                 {
00183                   hr = value(j,k) ; value(j,k) = value(r,k) ; value(r,k) = hr;
00184                 }
00185               i = p[j] ; p[j] = p[r] ; p[r] = i;
00186             }
00187           
00188           hr = 1./value(j,j);
00189           value(j,j) = hr;
00190           for (k=0;k<N;k++)
00191             {
00192               if (k==j) continue;
00193               for (i=0;i<N;i++)
00194                 {
00195                   if (i==j) continue;
00196                   value(i,k) -= value(i,j)*value(j,k)*hr;
00197                 }
00198             }
00199           for (i=0;i<N;i++)
00200             {
00201               value(i,j) *= hr;
00202               value(j,i) *= -hr;
00203             }
00204           value(j,j) = hr;
00205         }
00206       std::vector<double> hv(N);
00207       for (i=0;i<N;i++)
00208         {
00209           for (k=0;k<N;k++) hv[p[k]] = value(i,k);
00210           for (k=0;k<N;k++) value(i,k) = hv[k];
00211         }
00212     }
00213   
00214 };
00215 }
00216 
00217 #endif
00218 

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