home *** CD-ROM | disk | FTP | other *** search
/ Aminet 18 / aminetcdnumber181997.iso / Aminet / dev / c / math_classes.lha / math_classes / matrix / libmatrix / lusolve.cc < prev    next >
Encoding:
C/C++ Source or Header  |  1996-10-05  |  774 b   |  47 lines
  1. #include "Matrix.h"
  2.  
  3. Matrix LUsolve(const LUmatrix& m1,const Matrix& b)
  4.     {
  5.     Matrix x(b.rows(),b.cols());
  6.     unsigned int i,j,k;
  7.     mtrxtype t1;
  8.     mtrxtype t2;
  9.     Matrix v(b.rows(),b.cols());
  10.     Matrix y(b.rows(),b.cols());
  11.     Matrix m(m1.lu.rows(),m1.lu.cols());
  12.     if(m.rows()!=b.rows() || m.cols()!=m.rows()) 
  13.         {
  14.         error("dimensions do not match");
  15.         }
  16.     for(k=1;k<=b.cols();k++)    
  17.         {
  18.         for(i=1;i<=b.rows();i++)
  19.             {
  20.             v(i,k)=b(m1.p[i],k);
  21.             }
  22.         m=getl(m1);
  23.         for(j=1;j<=m.rows();j++)
  24.             {
  25.             t1=0.0;
  26.             for(i=1;i<j;i++)
  27.                 {
  28.                 t1+=m(j,i)*y(i,k);
  29.                 }
  30.             t2=(v(j,k)-t1)/m(j,j);
  31.             y(j,k)=t2;
  32.             }
  33.         m=getu(m1);
  34.         for(j=m.rows();j>=1;j--)
  35.             {
  36.             t1=0.0;
  37.             for(i=m.cols();i>j;i--)
  38.                 {
  39.                 t1+=m(j,i)*x(i,k);
  40.                 }
  41.             t2=(y(j,k)-t1)/m(j,j);
  42.             x(j,k)=t2;
  43.             }
  44.         }
  45.     return x;
  46.     }
  47.