home *** CD-ROM | disk | FTP | other *** search
/ QBasic & Borland Pascal & C / Delphi5.iso / C / Samples / C-SSP.ARJ / RANK.C < prev    next >
Encoding:
Text File  |  1984-08-24  |  5.6 KB  |  250 lines
  1.    rank(mr,nr,nc,ncc,a,zmch,cons,nrank,ierr,lp,mp)
  2.  
  3.       /* purpose...solution of m equations in n unknowns.*/
  4.  
  5.       int *ierr,lp[],mp[],mr,nr,nc,ncc,*nrank;
  6.       float a[],zmch,cons;
  7.  
  8.    {
  9.       int m,mdn,mtr,mtt,mtra,mtxi,ncs,nrs,ntc,nx,ny,
  10.           i,ii,iia,ics,icf,iiad,iisb,irev,is,ix,itemp,
  11.           ia,iec,j,lpis,ntrad,nfin,nranka,nranks,
  12.           k,kf,ks,krs,kcs,kx,ky,ircnt,nrtad,ntr;
  13.       float temp,fj,fmr,x;
  14.       extern double fabs();
  15.  
  16.       if(nr > mr || ncc > mr)
  17.       {
  18.       *ierr = 2;
  19.       return;
  20.       }
  21.       fmr = mr;
  22.       mdn = mr - nc;
  23.       mtr = (nc-1) * mr + nc - 1;
  24.       if(nr < nc) ny = nr;
  25.        else ny = nc;
  26.       *nrank = 0;
  27.       *ierr = 0;
  28.       nrs = nr - 1;
  29.       ncs = nc - 1;
  30.       for(i = 0; i <= nrs; i++)
  31.        mp[i] = i + 1;
  32.       for(i = 0; i <= ncs; i++)
  33.        lp[i] = i + 1;
  34.       ntc = nc + ncc;
  35.       mtt = mr * nrs + ntc - 1;
  36.       ix = nr;
  37.       for(i = 0; i <= ny -1; i++)
  38.       {
  39.        ii = mr * i + i;
  40.        iisb = ii - mr;
  41.        iiad = ii + 1;
  42.        ics = i;
  43.        icf = i + (nr - 1) * mr;
  44.        iec = i * mr + ntc - 1;
  45.        iia = ii + mr;
  46.  reit: temp = 0.;
  47.  
  48.        for(j = ii; j <= iec; j++)
  49.        {
  50.         if(i != 0)
  51.         {
  52.          kf = j - mr;
  53.          ks = j - (i * mr);
  54.          kx = i * mr;
  55.          for(k = ks; k <= kf; k+= mr)
  56.          {
  57.           a[j] = a[j] - a[kx] * a[k];
  58.           kx = kx + 1;
  59.          }
  60.         }
  61.         if(j <= iec-ncc)
  62.         {
  63.          if(fabs(a[j]) > temp)
  64.          {
  65.          temp = fabs(a[j]);
  66.  
  67.          fj = j;
  68.          x = fj / fmr + .5;
  69.          nx = x;
  70.          nx = j - nx * mr;
  71.          }
  72.         }
  73.        }
  74.        if(i != nc)
  75.        {
  76.         if(nx != i)
  77.         {
  78.          itemp = lp[nx];
  79.          lp[nx] = lp[i];
  80.          lp[i] = itemp;
  81.          lpis = nx;
  82.          for(k = ics; k <= icf; k += mr)
  83.          {
  84.           temp = a[k];
  85.           a[k] = a[lpis];
  86.           a[lpis] = temp;
  87.           lpis = lpis + mr;
  88.          }
  89.         }
  90.        }
  91.        if((zmch - fabs(a[ii])) > 0)
  92.        {
  93.         mtxi = nr + i * mr;
  94.         if(ncc != 0)
  95.         {
  96.          for(j = mtxi; j <= iec; j++)
  97.          {
  98.           if(abs(a[j]) > cons)
  99.           {
  100.            mp[i] = 0 - abs(mp[i]);
  101.            a[j] = mp[i];
  102.            *ierr = 1;
  103.           }
  104.            else a[j] = abs(mp[i]);
  105.          }
  106.          if(mp[i] > 0) mp[i] = 0;
  107.         }
  108.         if(i == (ix-1)) goto r78;
  109.         itemp = mp[i];
  110.         mp[i] = mp[ix-1];
  111.         mp[ix-1] = itemp;
  112.         k = i;
  113.         for(j = ix -1; j <= iec; j++)
  114.         {
  115.          temp = a[k];
  116.          a[k] = a[j];
  117.          a[j] = temp;
  118.          k = k + 1;
  119.         }
  120.         ix = ix - 1;
  121.         goto reit;
  122.        }
  123.        *nrank = *nrank + 1;
  124.        for(j = iiad; j <= iec; j++)
  125.        {
  126.         a[j] = a[j] / a[ii];
  127.        }
  128.        if(i == 0) goto ibr;
  129.        {
  130.         if(i == (ix-1)) goto r78;
  131.         for(m = iia; m <= icf; m += mr)
  132.         {
  133.          if(m < mr) kx = 0;
  134.           else kx = m / mr;
  135.          kx = kx * mr;
  136.          for(ky = ics; ky <= iisb; ky += mr)
  137.          {
  138.           a[m] = a[m] - a[kx] * a[ky];
  139.           kx = kx + 1;
  140.          }
  141.         }
  142.        }
  143. ibr:  ;
  144.       }
  145.       if(ncc == 0) return;
  146.       nranka = *nrank + 1;
  147.       for(i = nranka; i <= ix; i++)
  148.       {
  149.        ia = nc + (i - 1) * mr;
  150.        iec = ia / mr +1;
  151.        iec = iec * (mr - 1) + ntc - 1;
  152.        for( j = ia; j <= iec; j++)
  153.        {
  154.         if( j < mr) ky = 0;
  155.          else ky = j / mr;
  156.         ky = ky * mr;
  157.  
  158.         for(kx = i; kx <= (iec-mdn); kx++)
  159.         {
  160.          a[j] = a[j] - a[kx] * a[ky];
  161.          ky = ky + mr;
  162.         }
  163.         if(fabs(a[j]) > cons)
  164.         {
  165.          mp[i-1] = 0 - abs(mp[i-1]);
  166.          a[j] = mp[i-1];
  167.          *ierr = 1;
  168.         }
  169.          else a[j] = abs(mp[i-1]);
  170.        }
  171.        if(mp[i-1] >  0) mp[i-1] = 0;
  172.       }
  173. r78:  if(ncc == 0) return;
  174.       nrtad = nc + (*nrank - 1) * mr;
  175.       nranks = *nrank - 1;
  176.       nfin = nc - *nrank;
  177.       ntr = mtr - nfin;
  178.       for(i = 1; i <= nranks; i++)
  179.       {
  180.        irev = nrtad - i * mr;
  181.        iec = irev / mr ;
  182.        iec = iec * mr + nc + ncc - 1;
  183.        krs = irev - i - nfin;
  184.        for(ircnt = irev; ircnt <= iec; ircnt ++)
  185.        {
  186.         kcs = ircnt + mr;
  187.         for(k = krs; k <= (iec-ncc); k++)
  188.         {
  189.         a[ircnt] = a[ircnt] - a[kcs] * a[k];
  190.         kcs = kcs + mr;
  191.         }
  192.        }
  193.       }
  194.  
  195.       for(i = 1; i <= nc; i++)
  196.       {
  197.        ics =lp[i-1] - 1;
  198.        icf = ics + (ncc - 1) * mr;
  199.        k = nc + ((i-1) * mr) - 1;
  200.        for(j = ics; j <= icf; j += mr)
  201.        {
  202.         k = k + 1;
  203.         if(i <= *nrank)
  204.         {
  205.         a[j] = a[k];
  206.         a[k] = mp[i-1];
  207.         }
  208.          else a[j] = 0;
  209.        }
  210.       }
  211.       for(i = 1; i <= nrs; i++)
  212.       {
  213.        mtxi = nc + (i - 1) * mr;
  214.         ix = fabs(a[mtxi]);
  215.         if(ix = i) goto ixe;
  216.         itemp = mp[i-1];
  217.         mp[i-1] = mp[ix-1];
  218.         mp[ix-1] = itemp;
  219.         j = nc + ix * mr;
  220.         iec = mtxi / mr ;
  221.         iec = iec * mr + nc + ncc - 1;
  222.         for(k = mtxi; k <= iec; k++)
  223.         {
  224.          temp = a[k];
  225.          a[k] = a[j];
  226.          a[j] = temp;
  227.          j = j + 1;
  228.         }
  229.  ixe: ;
  230.       }
  231.       if(*nrank != nc)
  232.       {
  233.        nranka = *nrank +1;
  234.        for(i = nranka; i <= nc; i++)
  235.         lp[i+1] = -lp[i+1];
  236.       }
  237.       ncs = nc - 1;
  238.       for(i = 1; i <= ncs; i++)
  239.       {
  240.        while (abs(lp[i-1]) != i)
  241.        {
  242.         ix = abs(lp[i-1]);
  243.         itemp = lp[i-1];
  244.         lp[i-1] = lp[ix-1];
  245.         lp[ix-1] = itemp;
  246.        }
  247.       }
  248.    }
  249.  
  250.