home *** CD-ROM | disk | FTP | other *** search
/ QBasic & Borland Pascal & C / Delphi5.iso / C / Samples / C-SSP.ARJ / MATINV.C < prev    next >
Encoding:
Text File  |  1984-08-27  |  3.8 KB  |  179 lines
  1.    matinv(isol,idsol,nr,nc,a,mca,kwa,det)
  2.  
  3.       /* this subroutine finds the inverse and/or solves
  4.          simultaneous equations, or neither, and
  5.          calculates a determinant of a real matrix.*/
  6.  
  7.       int *isol,*idsol,nr,nc,mca,kwa[];
  8.       float a[],*det;
  9.  
  10.    {
  11.       int i,ic,iec,ii,ip,ir,iric,ibmp,
  12.           j,jbmp,jj,jr,k,kbmp,kk,kser,l,
  13.           mad,mdiv,mser,mz,nes,net;
  14.  
  15.       float piv,psto;
  16.       extern double fabs();
  17.  
  18.       ir = nr;
  19.       *isol = 1;
  20.       *idsol = 1;
  21.       if(nr <= 0) goto er1;
  22.       if((ir - mca) > 0) goto er1;
  23.       ic = abs(nc);
  24.       if ((ic - ir) <  0) ic = ir;
  25.       ibmp = mca;
  26.       jbmp = 1;
  27.       kbmp = jbmp + ibmp;
  28.       nes = (ir-1) * ibmp + ir - 1;
  29.       net = nes + ic - ir;
  30.       if(nc == 0) goto er1;
  31.       if(nc < 0)
  32.       {
  33.        mdiv = 1;
  34.        iric = ir - ic;
  35.       }
  36.       else mdiv = 0;
  37.       mad = mdiv;
  38.       mser = 0;
  39.       kser = (ir-1) * mca;
  40.       mz = 0;
  41.       *det = 1.;
  42. reit: piv = 0.;
  43.       for(i = mser; i <= kser; i+=mca)
  44.       {
  45.        if((fabs(a[i])-piv) > 0.)
  46.        {
  47.        piv = fabs(a[i]);
  48.        ip = i;
  49.        }
  50.       }
  51.       if(piv == 0.) goto er2;
  52.       if(nc > 0)
  53.       {
  54.       i = ip / mca;
  55.       j = mser - (mser / mca);
  56.       jj = mser/mca;
  57.       ii = jj + (ip -mser)/mca;
  58.       kwa[jj] = ii;
  59.       }
  60.       else
  61.       {
  62.        i = ip;
  63.        j = mser;
  64.       }
  65.  
  66.       if((ip - (mser/mca)) < 0) goto er1;
  67.        else
  68.        {
  69.         if((ip - mser) > 0)
  70.         {
  71.          iec = (i / mca) * mca + ic - 1;
  72.          for(k = ip; k <= iec; k++)
  73.          {
  74.           psto = a[i];
  75.           a[i] = a[j];
  76.           a[j] = psto;
  77.           i = i + 1;
  78.           j = j + 1;
  79.          }
  80.          *det = - *det;
  81.         }
  82.         psto = a[mser];
  83.         *det = *det * psto;
  84.        }
  85.       if (*det == 0.)
  86.       {
  87.        *idsol = 1;
  88.        *isol = 2;
  89.        return;
  90.       }
  91.       psto = 1. / psto;
  92.       a[mser] = 1.;
  93.       i = mdiv;
  94.       iec = (i / mca) * mca + ic -1;
  95.       for(k = i; k <= iec; k++)
  96.       {
  97.        a[i] = a[i] * psto;
  98.        i = i + 1;
  99.       }
  100.       for(k = mz; k <= kser; k+=mca)
  101.       {
  102.        if(mz != mser)
  103.        {
  104.         i = mad;
  105.         j = mdiv;
  106.         psto = a[mz];
  107.         if(psto == 0.) goto kbr;
  108.         a[mz] = 0.;
  109.         iec = (j / mca) * mca + ic - 1;
  110.         for(kk = j; kk <= iec; kk++)
  111.         {
  112.          a[i] = a[i] - a[j] * psto;
  113.          j = j + 1;
  114.          i = i + 1;
  115.         }
  116.        }
  117. kbr:   mad = mad + mca;
  118.        if(mad > net)
  119.         mad = mad - (mad / mca) * mca + 1;
  120.        mz = mz + mca;
  121.        if( mz > net)
  122.         mz = mz - (mz / mca) * mca + 1;
  123.       }
  124.       kser = kser + 1;
  125.       if ((kser-nes) >  0) goto invr;
  126.       mser = mser + kbmp;
  127.       if(nc >= 0)
  128.       {
  129.        mdiv = mdiv + mca;
  130.        if( mdiv > net)
  131.         mdiv = mdiv - (mdiv / mca) * mca + 1;
  132.        mz = mser - (mser / mca) * mca;
  133.        mad = 0;
  134.        goto reit;
  135.       }
  136.       mdiv = mdiv + kbmp;
  137.       if(iric == 0)
  138.       {
  139.        mz = mser + mca;
  140.        if( mz > net)
  141.         mz = mz - (mz / mca) * mca + 1;
  142.       }
  143.         else mz = mser - (mser / mca) * mca;
  144.       mad = mz + 1;
  145.       goto reit;
  146. invr: if(nc <  0) return;
  147.       jr = ir;
  148.       for(kk = jr; kk >= 1; kk--)
  149.       {
  150.        if((kwa[jr - 1] - (jr-1)) < 0) goto er1;
  151.        if((kwa[jr - 1] - (jr-1)) > 0)
  152.        {
  153.         k = jr - 1;
  154.         j = k - 1 + (ir-1) * mca;
  155.         l = kwa[jr-1] - 1 + (ir-1) * mca;
  156.         for(kk = j; kk > k; kk+=mca)
  157.         {
  158.         psto = a[l-1];
  159.         a[l-1] = a[j-1];
  160.         a[j-1] = psto;
  161.         j = j - mca;
  162.         l = l - mca;
  163.         }
  164.         goto inv2;
  165.        }
  166. inv2:  ;
  167.       }
  168.       return;
  169. er1:  *isol = 3;
  170.       return;
  171. er2:  *det = 0.;
  172.       *isol = 2;
  173.       *idsol = 1;
  174.       return;
  175. er3:  *isol = 2;
  176.       *idsol = 2;
  177.       return;
  178.    }
  179.