1 #ifdef HAVE_STDLIB_H
2 #include <stdlib.h>
3 #endif
4 #include "yagi.h"
5 
mprove(double ** a,double ** alud,int n,int * indx,double * b,double * x)6 void mprove(double **a,double **alud,int n,int *indx,double *b,double *x)
7 {
8 	int j,i;
9 	double sdp;
10 	double *r;
11 
12 	r=dvector(1,n);
13 	for (i=1;i<=n;i++) {
14 		sdp = -b[i];
15 		for (j=1;j<=n;j++) sdp += a[i][j]*x[j];
16 		r[i]=sdp;
17 	}
18 	lubksb(alud,n,indx,r);
19 	for (i=1;i<=n;i++) x[i] -= r[i];
20 	free_dvector(r,1,n);
21 }
22