home *** CD-ROM | disk | FTP | other *** search
/ Aminet 18 / aminetcdnumber181997.iso / Aminet / dev / c / math_classes.lha / math_classes / matrix / Matrix.cc < prev    next >
Encoding:
C/C++ Source or Header  |  1996-09-17  |  13.9 KB  |  976 lines
  1.  
  2. #include <iostream.h>
  3. #include <iomanip.h>
  4. #include <stdlib.h>
  5. #include "matrix.h"
  6.  
  7.  
  8. double abs(double a)
  9.     {
  10.     return (a>0 ? a : -a);
  11.     }
  12.     
  13.  
  14. Matrix::Matrix(const Vector& v,int m)
  15.     {
  16.     unsigned int i;
  17.     if(v.size()==0)error("bad size");
  18.     if(m==0)  //row vector
  19.         {
  20.         r=1;
  21.         c=v.size();
  22.         vals=new mtrxtype[v.size()] ;
  23.         if(vals==0)error("out of mem");
  24.         for(i=1;i<=cols();i++)
  25.             {
  26.             (*this)(1,i)=v(i);
  27.             }
  28.         }
  29.     else if(m==1) //column vector
  30.         {
  31.         r=v.size();
  32.         c=1;
  33.         vals=new mtrxtype[v.size()] ;
  34.         if(vals==0)error("out of mem");
  35.         for(i=1;i<=rows();i++)
  36.             {
  37.             (*this)(i,1)=v(i);
  38.             }
  39.         }
  40.     else
  41.         {
  42.         error("bad mode");
  43.         }
  44.     }
  45.  
  46.  
  47. Vector::Vector(const Matrix& m1)
  48.     {
  49.     unsigned int i;
  50.     if(m1.rows()!=1 && m1.cols()!=1)
  51.         {
  52.         error("dimensions do not match");
  53.         }
  54.     if(m1.rows()==1)
  55.         {
  56.         len=m1.cols();
  57.         vals=new vectype[len];
  58.         for(i=1;i<=size();i++)
  59.             {
  60.             (*this)(i)=m1(1,i);
  61.             }
  62.         }
  63.     else if(m1.cols()==1)
  64.         {
  65.         len=m1.rows();
  66.         vals=new vectype[len];
  67.         for(i=1;i<=size();i++)
  68.             {
  69.             (*this)(i)=m1(i,1);
  70.             }
  71.         }
  72.     }
  73.  
  74.  
  75. Matrix::Matrix()
  76.     {
  77.     r=0;
  78.     c=0;
  79.     vals=0;
  80.     }
  81.     
  82.     
  83. Matrix::Matrix(unsigned int ro,unsigned int co)
  84.     {
  85.     if(ro==0 || co==0) error("bad size");
  86.     r=ro;
  87.     c=co;
  88.     vals=new mtrxtype[ro*co];
  89.     if(vals==0) error("out of mem");
  90.     for(unsigned int i=0;i<ro*co;i++)
  91.         {
  92.         vals[i]=0.0;
  93.         }
  94.     }
  95.  
  96.  
  97. Matrix::Matrix(const Matrix& m)
  98.     {
  99.     r=m.rows();
  100.     c=m.cols();
  101.     vals=new mtrxtype[rows()*cols()];
  102.     if(vals==0) error("out of mem");
  103.     for(unsigned int i=1;i<=rows();i++)
  104.         {
  105.         for(unsigned int j=1;j<=cols();j++)
  106.             {
  107.             (*this)(i,j)=m(i,j);
  108.             }
  109.         }
  110.     }
  111.     
  112. Matrix::Matrix(const Pmatrix& p)
  113.     {
  114.     unsigned int i;
  115.     r=p.rows();
  116.     c=p.rows();
  117.     vals=new mtrxtype [r*c];
  118.     if(vals==0) error("out of mem");
  119.     for(i=0;i<rows()*cols();i++)
  120.         {
  121.         vals[i]=0.0;
  122.         }
  123.     for(i=1;i<=rows();i++)
  124.         {
  125.         (*this)(i,p[i])=1.0;
  126.         }
  127.     }
  128.  
  129.  
  130. Matrix::Matrix(unsigned int ro,unsigned int co,mtrxtype *m)
  131.     {
  132.     unsigned int p=0;
  133.     unsigned int i,j;
  134.     if(ro==0 || co ==0) error("bad size");
  135.     r=ro;
  136.     c=co;
  137.     vals=new mtrxtype[ro*co];
  138.     for(i=1;i<=rows();i++)
  139.         {
  140.         for(j=1;j<=cols();j++)
  141.             {
  142.             (*this)(i,j)=m[p++];
  143.             }
  144.         }
  145.     }
  146.  
  147.  
  148. Matrix::~Matrix()
  149.     {
  150.     if(vals!=0)
  151.     delete[] vals;
  152.     }
  153.  
  154.  
  155. int operator==(const Matrix& m1,const Matrix& m2)
  156.     {
  157.     if(m1.rows()!=m2.rows() || m1.cols()!=m2.cols()) return(0);
  158.     for(unsigned int i=1;i<=m1.rows();i++)
  159.         {
  160.         for(unsigned int j=1;j<=m1.cols();j++)
  161.             {
  162.             if(m1(i,j)!=m2(i,j))
  163.                 {
  164.                 return(0);
  165.                 }
  166.             }
  167.         }
  168.     return(1);
  169.     }
  170.  
  171.  
  172. int operator!=(const Matrix& m1,const Matrix& m2)
  173.     {
  174.     if(m1.rows()!=m2.rows() || m1.cols()!=m2.cols()) return(1);
  175.     for(unsigned int i=1;i<=m1.rows();i++)
  176.         {
  177.         for(unsigned int j=1;j<=m1.cols();j++)
  178.             {
  179.             if(m1(i,j)!=m2(i,j))
  180.                 {
  181.                 return(1);
  182.                 }
  183.             }
  184.         }
  185.     return(0);
  186.     }
  187.  
  188.  
  189. Matrix operator+ (const Matrix& m1,const Matrix& m2)
  190.     {
  191.     Matrix m(m1.rows(),m1.cols());
  192.     if(m1.rows()!=m2.rows() || m1.cols()!=m2.cols())
  193.         {
  194.         error("dimensions do not match");
  195.         }
  196.     for(unsigned int i=1;i<=m.rows();i++)
  197.         {
  198.         for(unsigned int j=1;j<=m.cols();j++)
  199.             {
  200.             m(i,j)=m1(i,j)+m2(i,j);
  201.             }
  202.         }
  203.     return m;
  204.     }    
  205.  
  206.  
  207. Matrix operator- (const Matrix& m1,const Matrix& m2)
  208.     {
  209.     Matrix m(m1.rows(),m1.cols());
  210.     if(m1.rows()!=m2.rows() || m1.cols()!=m2.cols())
  211.         {
  212.         error("dimensions do not match");
  213.         }
  214.     for(unsigned int i=1;i<=m.rows();i++)
  215.         {
  216.         for(unsigned int j=1;j<=m.cols();j++)
  217.             {
  218.             m(i,j)=m1(i,j)-m2(i,j);
  219.             }
  220.         }
  221.     return m;
  222.     }    
  223.  
  224.  
  225. void print(const Matrix& m)
  226.     {
  227.     for(unsigned int i=1;i<=m.rows();i++)
  228.         {
  229.         for(unsigned int j=1;j<=m.cols();j++)
  230.             {
  231.             cout.setf(ios::left);
  232.             cout.setf(ios::adjustfield);
  233.             cout.setf(ios::fixed);
  234.             cout <<setw(8) <<setprecision(4) <<m(i,j);
  235.             }
  236.         cout << "\n";
  237.         }
  238.     }
  239.  
  240.  
  241. LUmatrix& LUmatrix::operator= (const LUmatrix& m)
  242.     {
  243.     lu=m.lu;
  244.     p=m.p;
  245.     return *this;
  246.     }
  247.  
  248.  
  249. Matrix& Matrix::operator= (const Matrix& m)
  250.     {
  251.     if(rows()==0)
  252.         {
  253.         r=m.rows();
  254.         c=m.cols();
  255.         vals=new mtrxtype[rows()*cols()];
  256.         if(vals==0) error("out of mem");
  257.         for(unsigned int i=1;i<=rows();i++)
  258.             {
  259.             for(unsigned int j=1;j<=cols();j++)
  260.                 {
  261.                 (*this)(i,j)=m(i,j);
  262.                 }
  263.             }
  264.         return *this;
  265.         }
  266.     else
  267.         {
  268.         if(rows()!=m.rows() || cols()!=m.cols())error("dimensions do not match");
  269.         for(unsigned int i=1;i<=rows();i++)
  270.             {
  271.             for(unsigned int j=1;j<=cols();j++)
  272.                 {
  273.                 (*this)(i,j)=m(i,j);
  274.                 }
  275.             }
  276.         return *this;
  277.         }
  278.     }    
  279.  
  280.  
  281. Matrix& Matrix::operator+=(const Matrix& m)
  282.     {
  283.     if(rows()!=m.rows() || cols()!=m.cols())error("dimensions do not match");
  284.     for(unsigned int i=0;i<=rows();i++)
  285.         {
  286.         for(unsigned int j=0;j<=cols();j++)
  287.             {
  288.             (*this)(i,j)+=m(i,j);
  289.             }
  290.         }
  291.     return *this;
  292.     }
  293.  
  294. Matrix& Matrix::operator-=(const Matrix& m)
  295.     {
  296.     if(rows()!=m.rows() || cols()!=m.cols())error("dimensions do not match");
  297.     for(unsigned int i=0;i<=rows();i++)
  298.         {
  299.         for(unsigned int j=0;j<=cols();j++)
  300.             {
  301.             (*this)(i,j)-=m(i,j);
  302.             }
  303.         }
  304.     return *this;
  305.     }
  306.  
  307.  
  308.  
  309.  
  310. Matrix& Matrix::operator*=(mtrxtype a)
  311.     {
  312.     for(unsigned int i=1;i<=rows();i++)
  313.         {
  314.         for(unsigned int j=1;j<=rows();j++)
  315.             {
  316.             (*this)(i,j)*=a;
  317.             }
  318.         }
  319.     return *this;
  320.     }
  321.  
  322.  
  323. Vector operator* (const Matrix& m1,const Vector& v1)
  324.     {
  325.     Vector a(m1.rows());
  326.     unsigned int i,k;
  327.     mtrxtype t;
  328.     if(m1.cols()!=v1.size())
  329.         {
  330.         error("dimensions do not match");
  331.         }
  332.     for(i=1;i<=m1.rows();i++)
  333.         {
  334.         t=0.0;
  335.         for(k=1;k<=m1.cols();k++)
  336.             {
  337.             t+=m1(i,k)*v1(k);
  338.             }
  339.         a(i)=t;
  340.         }
  341.     return a;
  342.     }
  343.  
  344.  
  345. Vector operator* (const Vector& v1,const Matrix& m1)
  346.     {
  347.     Vector a(m1.cols());
  348.     unsigned int i,k;
  349.     mtrxtype t;
  350.     if(v1.size()!=m1.rows())
  351.         {
  352.         error("dimensions do not match");
  353.         }
  354.     for(i=1;i<=m1.cols();i++)
  355.         {
  356.         t=0.0;
  357.         for(k=1;k<=m1.rows();k++)
  358.             {
  359.             t+=m1(k,i)*v1(k);
  360.             }
  361.         a(i)=t;
  362.         }
  363.     return a;
  364.     }
  365.  
  366.  
  367. Matrix operator* (const Matrix& m1,const Matrix& m2)
  368.     {
  369.     Matrix a(m1.rows(),m2.cols());
  370.     unsigned int i,j,k;
  371.     mtrxtype t;
  372.     if(m1.cols() != m2.rows()) 
  373.         {
  374.         error("dimensions do not match");
  375.         }
  376.     for(i=1;i<=m1.rows();i++)
  377.         {
  378.         for(j=1;j<=m2.cols();j++)
  379.             {
  380.             t=0.0;
  381.             for(k=1;k<=m2.rows();k++)
  382.                 {
  383.                 t+=m1(i,k)*m2(k,j);
  384.                 }
  385.             a(i,j)=t;
  386.             }
  387.         }
  388.     return a;
  389.     }
  390.  
  391.  
  392. Matrix operator* (const Matrix& m,mtrxtype d)
  393.     {
  394.     Matrix a(m);
  395.     unsigned int i,j;
  396.     for(i=1;i<=m.rows();i++)
  397.         {
  398.         for(j=1;j<=m.cols();j++)
  399.             {
  400.             a(i,j)*=d;
  401.             }
  402.         }    
  403.     return a;
  404.     }
  405.     
  406.  
  407. Matrix operator* (mtrxtype d,const Matrix& m)
  408.     {
  409.     Matrix a(m);
  410.     unsigned int i,j;
  411.     for(i=1;i<=m.rows();i++)
  412.         {
  413.         for(j=1;j<=m.cols();j++)
  414.             {
  415.             a(i,j)*=d;
  416.             }
  417.         }    
  418.     return a;
  419.     }
  420.  
  421.  
  422. void Matrix::rowexchange(unsigned int r1,unsigned int r2)
  423.     {
  424.     unsigned int i;
  425.     mtrxtype a;
  426.     if(r1 > rows() || r1==0 || r2 > rows() || r2==0) error("bad subscript\n");
  427.     for(i=1;i<=cols();i++)
  428.         {
  429.         a=(*this)(r1,i);
  430.         (*this)(r1,i)=(*this)(r2,i);
  431.         (*this)(r1,i)=a;
  432.         }
  433.     }
  434.  
  435.  
  436. void Matrix::setcol(unsigned int co,Vector v)
  437.     {
  438.     for(unsigned int i=1;i<=rows();i++)
  439.         {
  440.         (*this)(i,co)=v(i);
  441.         }
  442.     }
  443.  
  444.  
  445. Vector Matrix::getcol (unsigned int co) const
  446.     {
  447.     Vector v(rows());
  448.     for(unsigned int i=1;i<=rows();i++)
  449.         {
  450.         v(i)=(*this)(i,co);
  451.         }
  452.     return v;
  453.     }
  454.     
  455. void Matrix::setrow(unsigned int ro,Vector v)
  456.     {
  457.     for(unsigned int i=1;i<=cols();i++)
  458.         {
  459.         (*this)(ro,i)=v(i);
  460.         }
  461.     }
  462.  
  463.  
  464. Vector Matrix::getrow (unsigned int ro) const
  465.     {
  466.     Vector v(rows());
  467.     for(unsigned int i=1;i<=cols();i++)
  468.         {
  469.         v(i)=(*this)(ro,i);
  470.         }
  471.     return v;
  472.     }
  473.  
  474.  
  475. Matrix transpose(const Matrix& m)
  476.     {
  477.     Matrix a(m.cols(),m.rows());
  478.     unsigned int i,j;
  479.     if(m.cols() ==0 || m.rows()==0)error("bad size");
  480.     for (i=1;i<=m.rows();i++)
  481.         {
  482.         for(j=1;j<=m.cols();j++)
  483.             {
  484.             a(j,i)=m(i,j);
  485.             }
  486.         }
  487.     return a;
  488.     }
  489.     
  490.     
  491. Vector solve(const Matrix& m1,const Vector& v1)
  492.     {
  493.     Vector x(v1.size());
  494.     unsigned int maxr;
  495.     mtrxtype maxv;
  496.     unsigned int i,j,k,i2;
  497.     mtrxtype t1,t3;
  498.     mtrxtype t2;
  499.     Matrix m(m1);
  500.     Vector v(v1);
  501.     if(m1.rows()!=v.size() || m1.cols()!=m1.rows()) 
  502.         {
  503.         error("dimensions do not match");
  504.         }
  505.     for(k=1;k<m.rows();k++)
  506.         {
  507.         maxr=k;
  508.         maxv=m(k,k);
  509.         for(i2=k;i2<=m.rows();i2++)
  510.             {
  511.             t3=m(i2,k);
  512.             if (abs(t3)>abs(maxv))
  513.                 {
  514.                 maxr=i2;
  515.                 maxv=t3;
  516.                 }
  517.             }
  518.         if (maxr!=k)
  519.             {
  520.             m.rowexchange(k,maxr);
  521.             t2=v(k);
  522.             v(k)=v(maxr);
  523.             v(maxr)=t2;
  524.             }
  525.         for(j=k+1;j<=m.rows();j++)
  526.             {
  527.             t1=m(j,k)/m(k,k);
  528.             for(i=k;i<=m.cols();i++)
  529.                 {
  530.                 t3=m(j,i)-t1*m(k,i);
  531.                 m(j,i)=t3;
  532.                 }
  533.             t2=v(j)-t1*v(k);
  534.             v(j)=t2;
  535.             }
  536.         }
  537.     for(j=m.rows();j>=1;j--)
  538.         {
  539.         t1=0.0;
  540.         for(i=m.cols();i>j;i--)
  541.             {
  542.             t1+=m(j,i)*x(i);
  543.             }
  544.         t2=(v(j)-t1)/m(j,j);
  545.         x(j)=t2;
  546.         }
  547.     return v;
  548.     }
  549.  
  550.  
  551. mtrxtype determinant(const Matrix& m1)
  552.     {
  553.     unsigned int k,i2,i,j;
  554.     unsigned int c;
  555.     unsigned int maxr;
  556.     mtrxtype maxv;
  557.     Matrix m(m1);
  558.     mtrxtype t1,t2,t3;
  559.     if(m1.rows()!=m1.cols())error("dimensions do not match");
  560.     c=1;
  561.     for(k=1;k<m.rows();k++)
  562.         {
  563.         maxr=k;
  564.         maxv=m(k,k);
  565.         for(i2=k;i2<=m.rows();i2++)
  566.             {
  567.             t3=m(k,i2);
  568.             if (abs(t3)>abs(maxv))
  569.                 {
  570.                 maxr=i2;
  571.                 maxv=t3;
  572.                 }
  573.             }
  574.         if (maxr!=k)
  575.             {
  576.             m.rowexchange(k,maxr);
  577.             c*=-1;
  578.             }
  579.         for(j=k+1;j<=m.rows();j++)
  580.             {
  581.             t1=m(j,k)/m(k,k);
  582.             for(i=k;i<=m.cols();i++)
  583.                 {
  584.                 t2=m(j,i)-t1*m(k,i);
  585.                 m(j,i)=t2;
  586.                 }
  587.             }
  588.         }
  589.     t1=1.0;
  590.     for(i=1;i<=m.cols();i++)
  591.         {
  592.         t1*=m(i,i);
  593.         }
  594.     return(t1*c);
  595.     }
  596.     
  597.     
  598. Matrix inverse(const Matrix& a)
  599.     {
  600.     Matrix m(a.rows(),a.cols());
  601.     unsigned int i,j;
  602.     Vector v(a.cols());
  603.     Vector x(a.cols());
  604.     LUmatrix b;
  605.     if(a.rows()!=a.cols())
  606.         {
  607.         error("dimensions do not match");
  608.         }
  609.     b=LUdecompose(a);
  610.     for(i=1;i<=m.rows();i++)
  611.         {
  612.         for(j=1;j<=v.size();j++)
  613.             {
  614.             v(j)=0.0;
  615.             }
  616.         v(i)=1.0;
  617.         x=LUsolve(b,v);
  618.         m.setcol(i,x);
  619.         }
  620.     return m;
  621.     }
  622.     
  623.  
  624. void Matrix::setunity(void)
  625.     {
  626.     if(rows()!=cols()) error("dimensions do not match");
  627.     for(unsigned int i=1;i<=rows();i++)
  628.         {
  629.         for(unsigned int j=1;j<=cols();j++)
  630.             {
  631.             if(i==j)
  632.                 {
  633.                 (*this)(i,j)=1.0;
  634.                 }
  635.             else
  636.                 {
  637.                 (*this)(i,j)=0.0;
  638.                 }
  639.             }
  640.         }
  641.     }
  642.  
  643.  
  644. Matrix unity(unsigned int n)
  645.     {
  646.     Matrix a(n,n);
  647.     if(n==0)error("bad size");
  648.     for(unsigned int i=1;i<=n;i++)
  649.         {
  650.         a(i,i)=1.0;
  651.         }
  652.     return a;
  653.     }
  654.  
  655.  
  656. void swap(Matrix& m1,Matrix& m2)
  657.     {
  658.     unsigned int a,b;
  659.     mtrxtype * c;
  660.     a=m1.r;
  661.     b=m1.c;
  662.     c=m1.vals;
  663.     m1.r=m2.r;
  664.     m1.c=m2.c;
  665.     m1.vals=m2.vals;
  666.     m2.r=a;
  667.     m2.c=b;
  668.     m2.vals=c;
  669.     }
  670.  
  671.  
  672. LUmatrix LUdecompose(const Matrix& m1)
  673.     {
  674.     LUmatrix m;
  675.     unsigned int maxr;
  676.     mtrxtype maxv;
  677.     unsigned int i,j,k,i2;
  678.     unsigned int t2;
  679.     mtrxtype t1,t3;
  680.     if(m1.cols()!=m1.rows()) 
  681.         {
  682.         error("dimensions do not match");
  683.         }
  684.     m.lu=m1;
  685.     m.p=Pmatrix(m1.rows());
  686.     for(k=1;k<m.lu.rows();k++)
  687.         {
  688.         maxr=k;
  689.         maxv=m.lu(k,k);
  690.         for(i2=k+1;i2<=m.lu.rows();i2++)
  691.             {
  692.             t3=m.lu(i2,k);
  693.             if (abs(t3)>abs(maxv))
  694.                 {
  695.                 maxr=i2;
  696.                 maxv=t3;
  697.                 }
  698.             }
  699.         if (maxr!=k)
  700.             {
  701.             m.lu.rowexchange(k,maxr);
  702.             t2=m.p[k];
  703.             m.p[k]=m.p[maxr];
  704.             m.p[maxr]=t2;
  705.             }
  706.         for(j=k+1;j<=m.lu.rows();j++)
  707.             {
  708.             t1=m.lu(j,k)/m.lu(k,k);
  709.             for(i=k+1;i<=m.lu.cols();i++)
  710.                 {
  711.                 t3=m.lu(j,i)-t1*m.lu(k,i);
  712.                 m.lu(j,i)=t3;
  713.                 }
  714.             }
  715.         }
  716.     for(i=1;i<=m1.rows();i++)
  717.         {
  718.         m.lu.setrow(i,m1.getrow(m.p[i]));
  719.         }
  720.     for(k=1;k<m.lu.rows();k++)
  721.         {
  722.         for(j=k+1;j<=m.lu.rows();j++)
  723.             {
  724.             t1=m.lu(j,k)/m.lu(k,k);
  725.             for(i=k+1;i<=m.lu.cols();i++)
  726.                 {
  727.                 t3=m.lu(j,i)-t1*m.lu(k,i);
  728.                 m.lu(j,i)=t3;
  729.                 }
  730.             m.lu(j,k)=t1;
  731.             }
  732.         }
  733.     return m;
  734.     }
  735.     
  736.     
  737. Matrix getl(const LUmatrix& m)
  738.     {
  739.     Matrix a(m.lu);
  740.     unsigned int i,j;
  741.     for(i=1;i<=a.rows();i++)
  742.         {
  743.         for(j=i+1;j<=a.cols();j++)
  744.             {
  745.             a(i,j)=0.0;
  746.             }
  747.         }
  748.     for(i=1;i<=a.rows();i++)
  749.         {
  750.         a(i,i)=1.0;
  751.         }
  752.     return a;
  753.     }
  754.     
  755.     
  756. Matrix getu(const LUmatrix& m)
  757.     {
  758.     Matrix a(m.lu);
  759.     unsigned int i,j;
  760.     for(i=2;i<=a.rows();i++)
  761.         {
  762.         for(j=1;j<i;j++)
  763.             {
  764.             a(i,j)=0.0;
  765.             }
  766.         }
  767.     return a;
  768.     }
  769.  
  770.  
  771. Matrix getp(const LUmatrix& m)
  772.     {
  773.     Matrix a(m.p);
  774.     return a;
  775.     }
  776.  
  777.  
  778. Vector solve(const LUmatrix& m1,const Vector& b)
  779.     {
  780.     return LUsolve(m1,b);
  781.     }
  782.  
  783.  
  784. Vector LUsolve(const LUmatrix& m1,const Vector& b) return x(b.size());
  785.     {
  786.     Vector x(b.size());
  787.     unsigned int i,j;
  788.     mtrxtype t1;
  789.     mtrxtype t2;
  790.     Vector v(b.size());
  791.     Vector y(b.size());
  792.     Matrix m(m1.lu.rows(),m1.lu.cols());
  793.     if(m.rows()!=b.size() || m.cols()!=m.rows()) 
  794.         {
  795.         error("dimensions do not match");
  796.         }
  797.     for(i=1;i<=b.size();i++)
  798.         {
  799.         v(i)=b(m1.p[i]);
  800.         }
  801.     m=getl(m1);
  802.     for(j=1;j<=m.rows();j++)
  803.         {
  804.         t1=0.0;
  805.         for(i=1;i<j;i++)
  806.             {
  807.             t1+=m(j,i)*y(i);
  808.             }
  809.         t2=(v(j)-t1)/m(j,j);
  810.         y(j)=t2;
  811.         }
  812.     m=getu(m1);
  813.     for(j=m.rows();j>=1;j--)
  814.         {
  815.         t1=0.0;
  816.         for(i=m.cols();i>j;i--)
  817.             {
  818.             t1+=m(j,i)*x(i);
  819.             }
  820.         t2=(y(j)-t1)/m(j,j);
  821.         x(j)=t2;
  822.         }
  823.     return x;
  824.     }
  825.  
  826.  
  827. Pmatrix::Pmatrix()
  828.     {
  829.     rws=0;
  830.     r=NULL;
  831.     }
  832.  
  833.  
  834. Pmatrix::Pmatrix(unsigned int s)
  835.     {
  836.     unsigned int i;
  837.     if(s<1)error("bad size");
  838.     rws=s;
  839.     r=new unsigned int[s];
  840.     if(r==NULL)error("out of store");
  841.     for (i=0;i<s;i++)
  842.         {
  843.         r[i]=i+1;
  844.         }
  845.     }
  846.     
  847.     
  848. Pmatrix::Pmatrix(const Pmatrix& m)
  849.     {
  850.     unsigned int i;
  851.     if(m.rows()==0)
  852.         {
  853.         rws=0;
  854.         r=NULL;
  855.         }
  856.     else
  857.         {
  858.         rws=m.rws;
  859.         r=new unsigned int[rws];
  860.         if(r==NULL)error("out of store");
  861.         for(i=0;i<rws;i++)
  862.             {
  863.             r[i]=m.r[i];
  864.             }
  865.         }
  866.     }
  867.     
  868.     
  869. Pmatrix::~Pmatrix()
  870.     {
  871.     if (rws!=0)
  872.         {
  873.         delete [] r;
  874.         }
  875.     }
  876.  
  877.  
  878. Pmatrix& Pmatrix::operator= (const Pmatrix& m)
  879.     {
  880.     unsigned int i;
  881.     if(rws!=0) 
  882.         {
  883.         delete[] r;
  884.         }
  885.     r=new unsigned int[rws=m.rws];
  886.     for(i=0;i<rws;i++)
  887.         {
  888.         r[i]=m.r[i];
  889.         }
  890.     return *this;
  891.     }
  892.     
  893.     
  894. Pmatrix inverse(const Pmatrix& m)
  895.     {
  896.     Pmatrix a(m.rws);
  897.     unsigned int i;
  898.     for(i=0;i<m.rws;i++)
  899.         {
  900.         a.r[m.r[i]]=i;
  901.         }
  902.     return a;
  903.     }
  904.  
  905.  
  906. Pmatrix transpose(const Pmatrix& m)
  907.     {
  908.     Pmatrix a(m.rws);
  909.     unsigned int i;
  910.     for(i=0;i<m.rws;i++)
  911.         {
  912.         a.r[m.r[i]]=i;
  913.         }
  914.     return a;
  915.     }
  916.  
  917.  
  918. Matrix operator* (const Pmatrix& p,const Matrix& m)
  919.     {
  920.     Matrix a(m.rows(),m.cols());
  921.     unsigned int i;
  922.     if(m.rows()!=p.rows() || m.rows() !=m.cols()) error("dimensions do not match");
  923.     for(i=1;i<=p.rows();i++)
  924.         {
  925.         a.setrow(i,m.getrow(p[i]));
  926.         }
  927.     return a;
  928.     }
  929.     
  930.     
  931. Vector operator* (const Pmatrix& p,const Vector& v)
  932.     {
  933.     Vector a(v.size());
  934.     unsigned int i;
  935.     if(v.size()!=p.rws) error("dimensions do not match");
  936.     for(i=1;i<=p.rws;i++)
  937.         {
  938.         a(i)=v(p[i]);
  939.         }
  940.     return v;
  941.     }
  942.  
  943.  
  944. Pmatrix operator* (const Pmatrix& p1,const Pmatrix& p2)
  945.     {
  946.     Pmatrix a(p1.rws);
  947.     unsigned int i;
  948.     if(p1.rws!=p2.rws)error("dimensions do not match");
  949.     for(i=0;i<p1.rws;i++)
  950.         {
  951.         a.r[i]=p2.r[p1.r[i]];
  952.         }
  953.     return a;
  954.     }
  955.     
  956.  
  957. unsigned int& Pmatrix::operator[](unsigned int i) const
  958.     {
  959. #ifdef CHECKRANGE
  960.     if(i>rws || i<1)error("bad subscript");
  961. #endif
  962.     return r[i-1];
  963.     }
  964.     
  965. mtrxtype Pmatrix::operator()(unsigned int i,unsigned int j) const
  966.     {
  967. #ifdef CHECKRANGE
  968.     if(i>rows() || i==0 || j>rows() || j==0)error("bad subscript");
  969. #endif    
  970.     if(j==(*this)[i])
  971.         return 1.0;
  972.     else
  973.         return 0.0;
  974.     }
  975.         
  976.