home *** CD-ROM | disk | FTP | other *** search
/ QBasic & Borland Pascal & C / Delphi5.iso / C / Samples / C-SSP.ARJ / CROUT.C < prev    next >
Encoding:
Text File  |  1984-08-23  |  4.3 KB  |  193 lines
  1.    crout(mc,nr,ncc,a,zmch,dt,ierr,lp)
  2.  
  3.       /* crout (a) operates on a coefficient matrix to solve a system of */
  4.       /* simultaneous equations or to compute an inverse and (2)         */
  5.       /* and (b) computes a determinant.                                 */
  6.  
  7.       /* crout reduces the original matrix and right-hand sides until    */
  8.       /* upon completion the reduced matrix replaces the original        */
  9.       /* matrix and the solutions replace the right-hand sides.          */
  10.  
  11.       int *ierr,mc,nr,ncc,lp[];
  12.       float a[],*dt,zmch;
  13.  
  14.    {
  15.  
  16.       int i,j,k,kf,ks,kx,kkx,ky,m,mtx,mca,mcs,mdn,mtr,mtra,mtt,nrs,ntc,nx,
  17.           ii,iisb,iiad,icf,ics,iec,iia,itemp,iy,jf,lpis,ircnt,irev,ixs,
  18.           kcs,nrtad,krs,ix;
  19.  
  20.       float temp,fj,fmc,x;
  21.       extern double fabs();
  22.  
  23.  
  24.       if(nr > mc) goto er2;
  25.       fmc = mc;
  26.       mtx = mc * nr - 1;
  27.       mca = mc + 1;
  28.       mcs = mc - 1;
  29.       mdn = mc - nr;
  30.       mtr = (nr-1) * mc + nr-1;
  31.       mtra = nr;
  32.       *dt = 1.;
  33.       *ierr = 0;
  34.       nrs = nr - 1;
  35.  
  36.       for(i = 0; i <= nrs; i++)
  37.        lp[i] = i + 1;
  38.  
  39.       if(ncc < 0) goto er2;
  40.       if(ncc == 0)
  41.       {
  42.        ntc = nr + nr - 1;
  43.        mtt = mc * (nr -1) + ntc;
  44.        j = mtra;
  45.        for(k = nr; k <= ntc; k++)
  46.        {
  47.         for(kx = 0; kx <= nr-1; kx++)
  48.         {
  49.          kkx = k + kx * mc;
  50.          a[kkx] = 0.;
  51.         }
  52.         a[j] = 1.;
  53.         j = j + mca;
  54.        }
  55.       }
  56.       else
  57.       {
  58.        ntc = nr + ncc;
  59.        mtt = mc * (nr -1) + ntc;
  60.       }
  61.       if(ntc <= nr) goto er2;
  62.  
  63.       for(i = 0; i <= nrs; i++)
  64.       {
  65.        ii = mc * i + i;
  66.        iisb = ii - mc;
  67.        iiad = ii + 1;
  68.        icf = i + ((nr-1) * mc);
  69.        ics = i;
  70.        iia = ii + mc;
  71.        iec = i * mc + ntc ;
  72.        temp = 0.;
  73.  
  74.        for(j = ii; j <= iec; j++)
  75.        {
  76.         if(i != 0)
  77.         {
  78.          kf = j - mc;
  79.          ks = j - (i* mc) ;
  80.          kx = i * mc;
  81.          for(k = ks; k <= kf; k += mc)
  82.          {
  83.           a[j] = a[j] - a[kx] * a[k];
  84.           kx = kx + 1;
  85.          }
  86.         }
  87.         if(j <= (iec-ncc))
  88.         {
  89.          if(fabs(a[j]) > temp)
  90.          {
  91.           temp = fabs(a[j]);
  92.           fj = j;
  93.           x = fj / fmc + .5;
  94.           nx = x;
  95.           nx = j - nx * mc;
  96.          }
  97.         }
  98.        }
  99.        if(i != nr)
  100.        {
  101.         if(nx != i)
  102.         {
  103.          itemp = lp[nx];
  104.          lp[nx] = lp[i];
  105.          lp[i] = itemp;
  106.          lpis = nx;
  107.  
  108.          for(k = ics; k <= icf; k+=mc)
  109.          {
  110.           kx = lpis + k * mc;
  111.           temp = a[k];
  112.           a[k] = a[lpis];
  113.           a[lpis] = temp;
  114.           lpis = lpis + mc;
  115.           *dt = - *dt;
  116.          }
  117.         }
  118.         *dt = *dt * a[ii];
  119.         if(zmch - fabs(a[ii]) > 0) goto er1;
  120.  
  121.         for(j = iiad; j <= iec; j++)
  122.         {
  123.          a[j] = a[j] / a[ii];
  124.         }
  125.  
  126.         if(i == 0) goto ibl;
  127.         {
  128.          if(i == (nr - 1)) goto ibr;
  129.  
  130.          for(m = iia; m <= icf; m += mc)
  131.          {
  132.           if(m < mc) kx = 0;
  133.            else kx = m / mc;
  134.           kx = kx * mc;
  135.           for(ky = ics; ky <= iisb; ky += mc)
  136.           {
  137.            a[m] = a[m] - a[kx] * a[ky];
  138.            kx = kx + 1;
  139.           }
  140.          }
  141.         }
  142.        }
  143. ibl:  ;
  144.       }
  145. ibr:  nrtad = mtr + 1;
  146.       for(i = 1; i <= nrs; i++)
  147.       {
  148.        irev = nrtad - (i * mc);
  149.        iec = irev / mc;
  150.        iec = iec * mc  + ntc ;
  151.        krs = irev - i;
  152.        for(ircnt = irev; ircnt <= iec; ircnt ++)
  153.        {
  154.         kcs = ircnt + mc;
  155.         if(kcs > mtx) kcs = kcs/mc +1;
  156.         for(k = krs; k <= (iec-mdn); k++)
  157.         {
  158.          a[ircnt] = a[ircnt] - a[kcs] * a[k];
  159.          kcs = kcs + mc;
  160.         if(kcs > mtx) kcs = kcs/mc +1;
  161.         }
  162.        }
  163.       }
  164.  
  165.       for(i = 1; i <= nrs; i++)
  166.       {
  167.        while (lp[i-1] != i)
  168.        {
  169.         nx = lp[i-1];
  170.         lp [i-1] = lp[nx-1];
  171.         lp[nx-1] = nx;
  172.         ixs = nr + (i-1) * mc;
  173.         iy = nr + (nx-1) * mc;
  174.         iec = ixs / mc;
  175.         iec = iec * mc + ntc ;
  176.         for(ix = ixs; ix <= iec; ix++)
  177.         {
  178.          temp = a[ix];
  179.          a[ix] = a[iy];
  180.          a[iy] = temp;
  181.          iy = iy + 1;
  182.         }
  183.        }
  184.       }
  185.       return ;
  186.  
  187. er2:  *ierr = 2;
  188.       *dt = 999999999.;
  189.       return;
  190.  
  191. er1:  *ierr = 1;
  192.    }
  193.