Nadesłany przez Codename, 08 kwietnia 2011 03:26
Kod przedstawiony poniżej przedstawia główną część rozwiązania problemu.Pobierz pełne rozwiązanie.
Jeżeli nie odpowiada Ci sposób formatowania kodu przez autora skorzystaj z pretty printer'a i dostosuj go automatycznie do siebie.
gs_2_c.c:
/* Metoda Gaussa-Seidela */
/* Jacek Czarniecki */
/* www.algorytm.org */
#include <stdio.h>
#define MAX 100
const int MAX_ITER = 200; // maksymalna liczba iteracji
const float EPS = 0.000001; // przybliżenie
const float INIT_X = 0.0; // wartości wektora poczatkowego
float L[MAX][MAX],D[MAX],U[MAX][MAX],b[MAX],x[MAX],tmp[MAX]; int n;
float I[MAX][MAX],TL[MAX][MAX];
float min(float x, float y){ return x<y?x:y; }
float f_abs(float x){ return x>0?x:-x; }
int read(){
int i, j;
scanf("%i",&n); // wczytanie rozmiaru macierzy
for(i=0;i<n;i++){
for(j=0;j<n;j++){ // wczytywanie bezpośrednio do L,D,U
if (i < j) {scanf("%f",&U[i][j]);}
else if (i > j) {scanf("%f",&L[i][j]);}
else {scanf("%f",&D[i]); }
}
} //na przekątnej nie powinno być zer
for(i=0;i<n;i++) scanf("%f",&b[i]);
return 0;
}
int copy(float *a, float *b,int n){
while(n--) b[n] = a[n];
return 0;
}
float cmp(float *a, float *b,int n, float eps){
float w,v=0; int flag = 0;
// sprawdzenie warunku zakonczenia: ||x(k)-x(k-1)|| <= epsilon
while(n--) { if(f_abs(a[n]-b[n]) > eps) {flag = 1; break;} }
return flag;
}
int test_zero(float *a, int n){
while(n--) if(a[n]==0) return 0;
return 1; // czy na przekatnej wystepuja zera?
}
float norma1(){
float n1,h; int i,j;
for (n1=i=0; i < n; i++) {
for (h=j=0; j < n; j++)
h += f_abs(TL[j][i]);
if (h > n1) n1 = h;
} // norma "jeden"
return n1;
}
float norma2(){
float n2,h; int i,j;
for (n2=i=0; i < n; i++) {
for (h=j=0; j < n; j++)
h += f_abs(TL[i][j]);
if (h > n2) n2 = h;
} // norma "nieskonczonosc"
return n2;
}
int zbiega(){
int i,j,k; float h;
// wyznaczenie (L+D)^-1 * N
for (i=0;i<n;i++) for (j=0;j<i;TL[i][j]=L[i][j++]); // TL <- L - kopia macierzy L
for (i=0;i<n;TL[i][i]=D[i++]); // TL <- L + D
for (i=0;i<n;i++) for (j=0;j<n;j++) I[i][j]=(i!=j?0:1);
for (i=0;i<n;i++){ // TL <- (L + D)^-1
for (j=0;j<=i;j++) I[i][j]/=D[i];
for (k=i+1;k<n;k++)
for (j=0;j<=k;j++) I[k][j]-=(I[i][j]*TL[k][i]);
}
for (i = 0; i < n; i++) { // TL <- (L + D)^-1 * N
for (j = 0; j < n; j++) {
TL[i][j] = 0;
for (k = 0; k <= i; k++)
TL[i][j] += (I[i][k] * U[k][j]);
}
}
return (min(norma1(), norma2()) <= 1)?1:0;
}
int solve(){
int i,j,k,counter;
read(); // wczytywanie macierzy
if(test_zero(D, n)) { // sprawdzanie przekątnej
if(zbiega()){ // warunek zbieżności
for (i=0; i<n; i++) {
D[i] = 1/D[i]; // D <- D^-1
b[i] *= D[i]; // b <- D^-1 * b
for (j=0; j<i; j++) L[i][j] *= D[i]; // L <- D^-1 * L
for (j=i+1; j<n; j++) U[i][j] *= D[i]; // U <- D^-1 * U
x[i] = INIT_X; // x <- [0,...,0]
}
counter = 0;
do{ // kolejne iteracje
copy(x,tmp,n); // tmp[] <- x[]
for (i=0; i<n; i++) x[i] = b[i];
for (i=0; i<n; i++) { // x = D^-1*b -
for (j=0; j<i; j++) x[i] -= L[i][j]*x[j]; // D^-1*L * x -
for (j=i+1; j<n; j++) x[i] -= U[i][j]*tmp[j]; // D^-1*U * x
}
counter++;
}while( cmp(x,tmp,n,EPS) && counter < MAX_ITER); // warunki zakonczenia iteracji
for(i=0;i<n;i++) printf("%.0f ",x[i]); putchar('\n'); // wypisanie wyniku
} else printf("Nie spelniono warunku zbieznosci\n");
} else printf("Na przekatnej wystepuja zera\n");
return 0;
}
int main(){
int t;
scanf("%i",&t);
while(t--) solve();
return 0;
}

