[WikiDyd] [TitleIndex] [WordIndex

Języki i metodyka programowania - kurs w sem. zimowym 2010/2011

Wykład XI: 14 grudnia 2010

Pisaliśmy program, który ma rozwiązywać układ równań algebraicznych metodą eliminacji Gaussa.

Rozpoczęliśmy od pliku nagłówkowego modułu obrabiającego macierze:

   1 #ifndef _MATRIX_H_
   2 #define _MATRIX_H_
   3 
   4 #include <stdio.h>
   5 
   6 typedef struct {
   7         int rn;
   8     int cn;
   9     double *e;
  10 } matrix_t;
  11 
  12 matrix_t * make_matrix( int rn, int cn );  
  13 /* tworzenie macierzy o zadanych wymiarach */
  14 
  15 matrix_t * read_matrix( FILE *in );        
  16 /* czytanie macierzy */
  17 
  18 void write_matrix( matrix_t *, FILE *out ); 
  19 /* pisanie (kompatybilne z czytaniem) */
  20 
  21 matrix_t * copy_matrix( matrix_t *s );  
  22 /* kopia macierzy */
  23 
  24 matrix_t * transpose_matrix( matrix_t * s ); 
  25 /* transpozycja */
  26 
  27 matrix_t * mull_matrix( matrix_t *, matrix_t * ); 
  28 /* mnożenie macierzy */
  29 
  30 matrix_t * lu_matrix( matrix_t * ); 
  31 /* metoda Gaussa przekształcenie macierzy prostokątnej (liczba w. <= liczba kol.) w postać schodkową */
  32 
  33 #endif
  34 

Implementacja:

   1 #include "matrix.h"
   2 #include <stdlib.h>
   3 #include <stdio.h>
   4 #include <string.h>
   5 
   6 matrix_t *make_matrix( int rn, int cn ) {
   7         matrix_t *new_mat= malloc( sizeof *new_mat );
   8     if( new_mat == NULL )
   9                 return NULL;
  10         if( (new_mat->e= malloc( (size_t)rn*(size_t)cn * sizeof * new_mat->e )) == NULL ) {
  11                 free( new_mat );
  12                 return NULL;
  13         }
  14     new_mat->rn= rn;
  15     new_mat->cn= cn;
  16         return new_mat;
  17 }
  18 
  19 void free_matrix( matrix_t *m ) {
  20         free( m->e );
  21         free( m );
  22 }
  23 
  24 matrix_t * read_matrix( FILE *in ) {
  25         int rn, cn;
  26     int i, j;
  27     matrix_t *new_mat;
  28         if( fscanf( in, "%d %d", &rn, &cn ) != 2 )
  29                 return NULL;
  30 
  31         if( (new_mat= make_matrix( rn, cn )) == NULL )
  32                 return NULL;
  33         for( i= 0; i < rn; i++ )
  34                 for( j= 0; j < cn; j++ )
  35                         if( fscanf( in, "%lf", &new_mat->e[i*cn+j] ) != 1 ) {
  36                                 free_matrix( new_mat );
  37                                 return NULL;
  38                         }
  39 
  40         return new_mat;
  41 }
  42 
  43 void write_matrix( matrix_t *m, FILE *out ) {
  44         int i,j;
  45         fprintf( out, "%d %d\n", m->rn, m->cn );
  46         for( i= 0; i < m->rn; i++ ) {
  47                 for( j= 0; j < m->cn-1; j++ )
  48                         fprintf( out, "%g ", m->e[i*m->cn+j] );
  49                 fprintf( out, "%g\n", m->e[i*m->cn+j] );
  50         }
  51 }
  52 
  53 matrix_t *copy_matrix( matrix_t *s ) {
  54         matrix_t *d= make_matrix( s->rn, s->cn );
  55         if( d != NULL ) {
  56         memcpy( d->e, s->e, s->rn*s->cn * sizeof * s->e );
  57                 /* int i;
  58                 for( i= 0; i < s->rn*s->cn; i++ )
  59                         *(d->e+i)= *(s->e+i);
  60         */
  61                 /* d->rn= s->rn; d->cn= s->cn; done in make_matrix */
  62         }
  63         return d;
  64 }
  65 
  66 matrix_t *transpose_matrix( matrix_t *s ) {
  67         matrix_t *d= make_matrix( s->cn, s->rn );
  68         if( d != NULL ) {
  69                 int i,j;
  70                 for( i= 0; i < s->rn; i++ )
  71                         for( j= 0; j < s->cn; j++ )
  72                           *(d->e+j*d->cn+i)= *(s->e+i*s->cn+j);
  73                 /* d->rn= s->cn; d->cn= s->rn; done in make_matrix */
  74         }
  75         return d;
  76 }
  77 
  78 matrix_t * mull_matrix( matrix_t *a, matrix_t *b ) {
  79     if( a->cn != b->rn )
  80                 return NULL;
  81     else {
  82             matrix_t *c = make_matrix( a->rn, b->cn );
  83         int i,j,k;
  84                 if( c == NULL )
  85                         return NULL;
  86         for( i= 0; i < c->rn; i++ ) {
  87                         for( j= 0; j < c->cn; j++ ) {
  88                                 double s= 0;
  89                                 for( k= 0; k < a->cn; k++ )
  90                                         s += *(a->e+i*a->cn+k) * *(b->e+k*b->cn+j);
  91                                 *(c->e+i*c->cn+j)= s;
  92                         }
  93                 }
  94                 return c;
  95         }
  96 }
  97 
  98 matrix_t *lu_matrix( matrix_t *a ) {
  99         matrix_t *c = copy_matrix( a );
 100         if( c != NULL ) {
 101                 int i,j,k;
 102         int cn= c->cn;
 103         int rn= c->rn;
 104         double *e = c->e;
 105                 for( k= 0; k < cn-1; k++ ) { /* eliminujemy (zerujemy) kolumnę nr k */
 106                         for( i= k+1; i < rn; i++ ) { /* pętla po kolejnych
 107                                             wierszach poniżej diagonalii k,k */
 108                                 double d= *(e+i*cn+k) / *(e+k*cn+k);
 109                                 for( j= k; j < cn; j++ )
 110                    *(e+i*cn+j) -= d * *(e+k*cn+j);
 111             }
 112                 }
 113         }
 114         return c;
 115 }

Ponieważ pisaliśmy tę funkcjonalność stopniowo, a więc także stopniowo ją testowaliśmy.

Najpierw czytanie, pisanie, kopia i transpozycja:

   1 #include "matrix.h"
   2 #include <stdio.h>
   3 
   4 int
   5 main( int argc, char **argv ) {
   6         FILE *in;
   7         if( argc > 1 && (in= fopen( argv[1], "r" )) != NULL ) {
   8                 matrix_t *m= read_matrix( in );
   9                 if( m != NULL ) {
  10                         printf( "Kopia:\n" );
  11                         write_matrix( copy_matrix(m), stdout );
  12                         printf( "\nTranspozycja:\n" );
  13                         write_matrix( transpose_matrix(m), stdout );
  14                 }
  15                 return 0;
  16         } else
  17                 return 1;
  18 }

Potem mnożenie:

   1 #include <stdio.h>
   2 #include "matrix.h"
   3 
   4 int
   5 main( int argc, char ** argv ) {
   6         if( argc == 3 ) {
   7                 matrix_t *a= read_matrix( fopen( argv[1], "r" ) );
   8                 matrix_t *b= read_matrix( fopen( argv[2], "r" ) );
   9                 matrix_t *c= mull_matrix( a, b );
  10         if( c == NULL )
  11                         printf( "nie da rady pomnożyć\n" );
  12                 else
  13                 write_matrix( c, stdout );
  14                 return 0;
  15         }
  16         else
  17                 return 1;
  18 }

I wreszcie eliminacja Gaussa:

   1 #include "matrix.h"
   2 #include <stdio.h>
   3 
   4 int
   5 main( int argc, char **argv ) {
   6         FILE *in;
   7         if( argc > 1 && (in= fopen( argv[1], "r" )) != NULL ) {
   8                 matrix_t *m= read_matrix( in );
   9                 if( m != NULL ) {
  10                         matrix_t *c = lu_matrix( m );
  11                         if( c != NULL ) {
  12                                 printf( "\nPo LU:\n" );
  13                                 write_matrix( c, stdout );
  14                         } else
  15                                 printf( "nie lula!\n" );
  16                 }
  17                 return 0;
  18         } else
  19                 return 1;
  20 }

Do testowania używaliśmy następujących plików z danymi:

2 3
1 2 3
10 20 30.3333333333333

10 11
1 2 3 4 5 6 7 8 9 0 11
2 2 3 4 5 6 7 8 9 0 12
3 2 3 4 5 6 7 8 9 0 13
4 2 3 4 5 6 7 8 9 0 14
5 2 3 4 5 6 7 8 9 0 15
6 2 3 4 5 6 7 8 9 0 16
7 2 3 4 5 6 7 8 9 0 17
8 2 3 4 5 6 7 8 9 0 18
9 2 3 4 5 6 7 8 9 0 19
0 2 3 4 5 6 7 8 9 0 20

3 2
1 10
2 20
3 30.3333

2 3
1 1 2
1 -1 0

Ten ostatni plik to macierzowy zapis układu równań

x + y = 2
x - y = 1

Używaliśmy go do testowania rozkładu Gaussa, ale stwierdziliśmy też, że prawdziwe testy wymagają generacji większych zbiorów danych. Napisaliśmy więc program, który generuje układ równań o dowolnej (niemal!) wielkości i rozwiązaniu x takim, że x_i = i.

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 
   4 int
   5 main( int argc, char **argv ) {
   6         if( argc > 1 ) {
   7                 int n= atoi( argv[1] );
   8                 int i,j;
   9                 double s;
  10                 printf( "%d %d\n", n, n+1 );
  11                 for( i= 0; i < n; i++ ) {
  12                         s= 0;
  13                         for( j= 0; j < n; j++ ) {
  14                                 double aij= (double)rand() / RAND_MAX;
  15                                 s+= aij * j;
  16                                 printf( "%g ", aij );
  17                         }
  18                         printf( "%g\n", s );
  19                 }
  20         }
  21         return 0;
  22 }

Bardzo proszę napisać w domu funkcję, która wykona drugą część metod Gaussa (podstawienie wstecz).


2015-09-23 06:44