/**********************************************************************
* Solve AX = B using a partial pivoting algorithm and reduced storage *
* ------------------------------------------------------------------- *
* SAMPLE RUN:                                                         *
*                                                                     *
* System to solve:                                                    *
*   2.0000  -1.0000   1.0000   7.0000 -12.5400    5.0000              *
*   1.0000   5.0000  -2.0000  -8.0000 100.0000    1.0000              *
*   3.0000  -2.0000   3.0000  45.0000  27.3333    3.0000              *
*  11.0000   0.5500  -2.0000  -4.0000   1.0000    4.0000              *
*  33.0000   2.0000  -3.0000   5.0000   7.3333  -10.0000              *
*                                                                     *
* Solution is:                                                        *
*  2.111496                                                           *
*  -25.829027                                                         *
*  8.174232                                                           *
*  -2.521467                                                          *
*  1.242104                                                           *
*                                                                     *
* ------------------------------------------------------------------- *
* Ref.: "Wassyng, A. - Solving Ax = b: A method with reduced storage  *
*        requirements, SIAM J. Numerical Analysis, vol.19 (1982),     *
*        pp. 197-204".                                                *
*                                                                     *
*                                  C++ Release By J-P Moreau, Paris.  *
*                                          (www.jpmoreau.fr)          *
**********************************************************************/
#include <stdio.h>
#include <math.h>

#define  SIZE  25

const double zero = 0.0;


double  a[SIZE], b[SIZE], soln[SIZE];
double  error, temp;
int     i, j, ierr, row;


void rowk(int n, int k, double *r);

void Swap(double *a, double *b) {
double temp;
  temp=*a;
  *a=*b;
  *b=temp;
}

void dple(int n, double *b, double *c, int *ierr) {
 
/*******************************************************************
          SOLUTION OF LINEAR EQUATIONS WITH REDUCED STORAGE
********************************************************************

  Uses the Henderson-Wassyng partial pivot algorithm.
  Wassyng, A. 'Solving Ax = b: A method with reduced storage requirements',
  SIAM J. Numerical Analysis, vol.19 (1982), pp. 197-204.

  The user must provide a routine ROWK to return the requested row of the
  matrix A.     */
  double bk, cj, dkj;
  double wk[SIZE];
  int    i, iflag, ij, ijold, ik, j, k, kjold, km1, kp1,
         last, lastm1, lcol, lcolp1, m, maxwk, mjold, nm1, np1;
  int    iwk[SIZE];

// SET THE NECESSARY CONSTANTS

  ierr = 0;
  maxwk = (n * n / 4) + n + 3;
  np1 = n + 1;
  k = 1;
  iflag = -1;

// GET THE FIRST COLUMN OF THE TRANSPOSED SYSTEM }

  rowk(n, 1, c);
  bk = b[1];

  if (n <= 1) {
    if (c[1] == zero) goto e130;
    c[1] = bk / c[1];
    return;
  }

// FIND THE PIVOT FOR COLUMN 1

  m = 1;
  for  (i=2; i<=n; i++)
    if (fabs(c[m]) < fabs(c[i]))  m = i;

  iwk[1] = m;
  Swap(&c[m],&c[1]);

  if (c[1] != zero) {
  
// FIND THE FIRST ELEMENTARY MATRIX AND STORE IT IN D
  
    for (i=2; i<=n; i++)
      wk[i-1] = -c[i] / c[1];
    wk[n] = bk / c[1];
  
// K LOOP - EACH K FOR A NEW COLUMN OF THE TRANSPOSED SYSTEM
  
    for (k=2; k<=n; k++) {

	  kp1 = k + 1;
      km1 = k - 1;
    
// GET COLUMN K
    
      rowk(n, k, c);
      for (j=1; j<=km1; j++) {
        m = iwk[j];
        Swap(&c[j],&c[m]);
      }
      bk = b[k];
    
      iflag = -iflag;
      lcol = np1 - k;
      lcolp1 = lcol + 1;
      lastm1 = 1;
      last = maxwk - n + k;
      if (k != 2) {
        lastm1 = maxwk - n + km1;
        if (iflag < 0)  last = last - n + k - 2;
        if (iflag > 0)  lastm1 = lastm1 - n + k - 3;
      }
    
// J LOOP - EFFECT OF COLUMNS 1 TO K-1 OF L-INVERSE
    
      for (j=1; j<=km1; j++) {
        cj = c[j];
        ij = (j-1) * lcolp1;
        if (j == km1)  ij = lastm1 - 1;
      
// I LOOP - EFFECT OF L-INVERSE ON ROWS K TO N+1
      
        for (i=k; i<=n; i++) {
          ij = ij + 1;
          c[i] = c[i] + wk[ij] * cj;
        }
        bk = bk - wk[ij+1] * cj;
      }
    
// K=N CASE
    
      m = k;
      if (k >= n) {
        if (c[k] == zero) goto e130;
        wk[last] = bk / c[k];
      }
      else {

// FIND THE PIVOT
      
        for (i=kp1; i<=n; i++)
          if (fabs(c[m]) < fabs(c[i]))  m = i;
      
        iwk[k] = m;
        Swap(&c[m],&c[k]);
        if (c[k] == zero) goto e130;
      
// FIND THE K-TH ELEMENTARY MATRIX
      
        ik = last;
        for (i=kp1; i<=n; i++) {
          wk[ik] = -c[i] / c[k];
          ik++;
        }
        wk[ik] = bk / c[k];
      }
    
// FORM THE PRODUCT OF THE ELEMENTARY MATRICES
    
      for (j=1; j<=km1; j++) {
        kjold = j * lcolp1 + k - np1;
        mjold = kjold + m - k;
        ij = (j-1) * lcol;
        ijold = ij + j;
        if (j == km1) {
          kjold = lastm1;
          mjold = lastm1 + m - k;
          ijold = lastm1;
        }
      
        ik = last - 1;
        dkj = wk[mjold];
        wk[mjold] = wk[kjold];
        for (i=kp1; i<=np1; i++) {
          ij = ij + 1;
          ijold = ijold + 1;
          ik = ik + 1;
          wk[ij] = wk[ijold] + wk[ik] * dkj;
		}
      }
    } // k loop
  
    last = maxwk;
    if (iflag < 0) last = maxwk - 2;
    wk[n] = wk[last];
  
// INSERT THE SOLUTION IN C
  
    for (i=1; i<=n; i++)  c[i] = wk[i];
  
    nm1 = n - 1;
    for (i=1; i<=nm1; i++) {
      k = n - i;
      m = iwk[k];
      Swap(&c[k],&c[m]);
    }
    return;
  } // if C[1]<>zero

// THE SYSTEM IS SINGULAR

e130: *ierr = k;

} //dple()


void rowk(int n, int k, double *r) {
  if (k==1) {
    r[1]=2.0; r[2]=-1.0; r[3]=1.0; r[4]=7.0; r[5]=-12.54;
  }
  else if (k==2) {
    r[1]=1.0; r[2]=5.0; r[3]=-2.0; r[4]=-8.0; r[5]=100.0;
  }
  else if (k==3) {
    r[1]=3.0; r[2]=-2.0; r[3]=3.0; r[4]=45.0; r[5]=27.3333;
}
  else if (k==4) {
    r[1]=11.0; r[2]=0.55; r[3]=-2.0; r[4]=-4.0; r[5]=1.0;
  }
  else if (k==5) {
    r[1]=33.0; r[2]=2.0; r[3]=-3.0; r[4]=5.0; r[5]=7.3333;
  }
} //rowk()


void main() {

// Define the right-hand side (n=5)
  b[1]=5.0; b[2]=1.0; b[3]=3.0; b[4]=4.0; b[5]=-10.0;

  printf("\n System to solve:\n");
  for (i=1; i<6; i++) {
    rowk(5,i,soln);
    for (j=1; j<6; j++) printf(" %8.4f",soln[j]);
    printf("  %8.4f\n", b[i]);
  }

  dple(5, b, soln, &ierr);

  if (ierr == 0) {
    printf("\n Solution is:\n");
    for (i=1; i<6; i++) printf("  %f\n", soln[i]);
  }
  else
    printf("\n Error = %d\n", ierr);
  printf("\n");
}

// end of file dple.cpp