home *** CD-ROM | disk | FTP | other *** search
/ Aminet 18 / aminetcdnumber181997.iso / Aminet / dev / c / math_classes.lha / math_classes / matrix / libmatrix / lumat_ludecomp.cc < prev    next >
Encoding:
C/C++ Source or Header  |  1996-09-18  |  1.0 KB  |  66 lines
  1. #include "Matrix.h"
  2.  
  3. LUmatrix LUdecompose(const Matrix& m1)
  4.     {
  5.     LUmatrix m;
  6.     unsigned int maxr;
  7.     mtrxtype maxv;
  8.     unsigned int i,j,k,i2;
  9.     unsigned int t2;
  10.     mtrxtype t1,t3;
  11.     if(m1.cols()!=m1.rows()) 
  12.         {
  13.         error("dimensions do not match");
  14.         }
  15.     m.lu=m1;
  16.     m.p=Pmatrix(m1.rows());
  17.     for(k=1;k<m.lu.rows();k++)
  18.         {
  19.         maxr=k;
  20.         maxv=m.lu(k,k);
  21.         for(i2=k+1;i2<=m.lu.rows();i2++)
  22.             {
  23.             t3=m.lu(i2,k);
  24.             if (abs(t3)>abs(maxv))
  25.                 {
  26.                 maxr=i2;
  27.                 maxv=t3;
  28.                 }
  29.             }
  30.         if (maxr!=k)
  31.             {
  32.             m.lu.rowexchange(k,maxr);
  33.             t2=m.p[k];
  34.             m.p[k]=m.p[maxr];
  35.             m.p[maxr]=t2;
  36.             }
  37.         for(j=k+1;j<=m.lu.rows();j++)
  38.             {
  39.             t1=m.lu(j,k)/m.lu(k,k);
  40.             for(i=k+1;i<=m.lu.cols();i++)
  41.                 {
  42.                 t3=m.lu(j,i)-t1*m.lu(k,i);
  43.                 m.lu(j,i)=t3;
  44.                 }
  45.             }
  46.         }
  47.     for(i=1;i<=m1.rows();i++)
  48.         {
  49.         m.lu.setrow(i,m1.getrow(m.p[i]));
  50.         }
  51.     for(k=1;k<m.lu.rows();k++)
  52.         {
  53.         for(j=k+1;j<=m.lu.rows();j++)
  54.             {
  55.             t1=m.lu(j,k)/m.lu(k,k);
  56.             for(i=k+1;i<=m.lu.cols();i++)
  57.                 {
  58.                 t3=m.lu(j,i)-t1*m.lu(k,i);
  59.                 m.lu(j,i)=t3;
  60.                 }
  61.             m.lu(j,k)=t1;
  62.             }
  63.         }
  64.     return m;
  65.     }
  66.