/********************************************************************** * 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 #include #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