[WikiDyd] [TitleIndex] [WordIndex

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

Wykład XII: 21 grudnia 2010

Skończyliśmy program, który ma rozwiązywać układ równań algebraicznych metodą eliminacji Gaussa. Dodaliśmy funkcję wykonującą podstawienie wstecz oraz dopisaliśmy wersję wykorzystującą wybór elementu głównego.

matrix.h

   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 
  14 matrix_t * read_matrix( FILE *in );
  15 
  16 void write_matrix( matrix_t *, FILE *out );
  17 
  18 matrix_t * copy_matrix( matrix_t *s );
  19 
  20 matrix_t * transpose_matrix( matrix_t * s );
  21 
  22 void xchg_rows( matrix_t *m, int i, int j );
  23 
  24 void xchg_cols( matrix_t *m, int i, int j );
  25 
  26 matrix_t * mull_matrix( matrix_t *, matrix_t * );
  27 
  28 matrix_t * ge_matrix( matrix_t * );
  29 
  30 int bs_matrix( matrix_t * );
  31 
  32 matrix_t * pivot_ge_matrix( matrix_t *, int *row_per );
  33 
  34 matrix_t * symm_pivot_ge_matrix( matrix_t *, int *per );
  35 
  36 int *pivot_get_inv_per( matrix_t *, int *row_per );
  37 
  38 #endif
  39 

matrix.c

   1 #include "matrix.h"
   2 #include <stdlib.h>
   3 #include <stdio.h>
   4 #include <string.h>
   5 
   6 matrix_t *
   7 make_matrix (int rn, int cn)
   8 {
   9   matrix_t *new_mat = malloc (sizeof *new_mat);
  10   if (new_mat == NULL)
  11     return NULL;
  12   if ((new_mat->e =
  13        malloc ((size_t) rn * (size_t) cn * sizeof *new_mat->e)) == NULL) {
  14     free (new_mat);
  15     return NULL;
  16   }
  17   new_mat->rn = rn;
  18   new_mat->cn = cn;
  19   return new_mat;
  20 }
  21 
  22 void
  23 free_matrix (matrix_t * m)
  24 {
  25   free (m->e);
  26   free (m);
  27 }
  28 
  29 matrix_t *
  30 read_matrix (FILE * in)
  31 {
  32   int rn, cn;
  33   int i, j;
  34   matrix_t *new_mat;
  35   if (fscanf (in, "%d %d", &rn, &cn) != 2)
  36     return NULL;
  37 
  38   if ((new_mat = make_matrix (rn, cn)) == NULL)
  39     return NULL;
  40   for (i = 0; i < rn; i++)
  41     for (j = 0; j < cn; j++)
  42       if (fscanf (in, "%lf", &new_mat->e[i * cn + j]) != 1) {
  43         free_matrix (new_mat);
  44         return NULL;
  45       }
  46 
  47   return new_mat;
  48 }
  49 
  50 void
  51 write_matrix (matrix_t * m, FILE * out)
  52 {
  53   int i, j;
  54   if (m == NULL) {
  55     fprintf (out, "Matrix is NULL\n");
  56     return;
  57   }
  58 
  59   fprintf (out, "%d %d\n", m->rn, m->cn);
  60   for (i = 0; i < m->rn; i++) {
  61     for (j = 0; j < m->cn - 1; j++)
  62       fprintf (out, "%g ", m->e[i * m->cn + j]);
  63     fprintf (out, "%g\n", m->e[i * m->cn + j]);
  64   }
  65 }
  66 
  67 matrix_t *
  68 copy_matrix (matrix_t * s)
  69 {
  70   matrix_t *d = NULL;
  71   if (s != NULL)
  72     d = make_matrix (s->rn, s->cn);
  73   if (d != NULL) {
  74     memcpy (d->e, s->e, s->rn * s->cn * sizeof *s->e);
  75     /* int i;
  76        for( i= 0; i < s->rn*s->cn; i++ )
  77        *(d->e+i)= *(s->e+i);
  78      */
  79     /* d->rn= s->rn; d->cn= s->cn; done in make_matrix */
  80   }
  81   return d;
  82 }
  83 
  84 matrix_t *
  85 transpose_matrix (matrix_t * s)
  86 {
  87   matrix_t *d = NULL;
  88   if (s != NULL)
  89     d = make_matrix (s->rn, s->cn);
  90   if (d != NULL) {
  91     int i, j;
  92     for (i = 0; i < s->rn; i++)
  93       for (j = 0; j < s->cn; j++)
  94         *(d->e + j * d->cn + i) = *(s->e + i * s->cn + j);
  95     /* d->rn= s->cn; d->cn= s->rn; done in make_matrix */
  96   }
  97   return d;
  98 }
  99 
 100 void
 101 xchg_rows (matrix_t * m, int i, int j)
 102 {
 103   if (m != NULL && i >= 0 && i < m->rn && j >= 0 && j < m->rn) {
 104     int k;
 105     double tmp;
 106     for (k = 0; k < m->cn; k++) {
 107       tmp = *(m->e + i * m->cn + k);
 108       *(m->e + i * m->cn + k) = *(m->e + j * m->cn + k);
 109       *(m->e + j * m->cn + k) = tmp;
 110     }
 111   }
 112 }
 113 
 114 void
 115 xchg_cols (matrix_t * m, int i, int j)
 116 {
 117   if (m != NULL && i >= 0 && i < m->cn && j >= 0 && j < m->cn) {
 118     int k;
 119     double tmp;
 120     for (k = 0; k < m->rn; k++) {
 121       tmp = *(m->e + k * m->cn + i);
 122       *(m->e + k * m->cn + i) = *(m->e + k * m->cn + j);
 123       *(m->e + k * m->cn + j) = tmp;
 124     }
 125   }
 126 }
 127 
 128 matrix_t *
 129 mull_matrix (matrix_t * a, matrix_t * b)
 130 {
 131   if (a == NULL || b == NULL || a->cn != b->rn)
 132     return NULL;
 133   else {
 134     matrix_t *c = make_matrix (a->rn, b->cn);
 135     int i, j, k;
 136     if (c == NULL)
 137       return NULL;
 138     for (i = 0; i < c->rn; i++) {
 139       for (j = 0; j < c->cn; j++) {
 140         double s = 0;
 141         for (k = 0; k < a->cn; k++)
 142           s += *(a->e + i * a->cn + k) * *(b->e + k * b->cn + j);
 143         *(c->e + i * c->cn + j) = s;
 144       }
 145     }
 146     return c;
 147   }
 148 }
 149 
 150 matrix_t *
 151 ge_matrix (matrix_t * a)
 152 {
 153   matrix_t *c = copy_matrix (a);
 154   if (c != NULL) {
 155     int i, j, k;
 156     int cn = c->cn;
 157     int rn = c->rn;
 158     double *e = c->e;
 159     for (k = 0; k < rn - 1; k++) {      /* eliminujemy (zerujemy) kolumnę nr k */
 160       for (i = k + 1; i < rn; i++) {    /* pętla po kolejnych
 161                                            wierszach poniżej diagonalii k,k */
 162         double d = *(e + i * cn + k) / *(e + k * cn + k);
 163         for (j = k; j < cn; j++)
 164           *(e + i * cn + j) -= d * *(e + k * cn + j);
 165       }
 166     }
 167   }
 168   return c;
 169 }
 170 
 171 int
 172 bs_matrix (matrix_t * a)
 173 {
 174   if (a != NULL) {
 175     int r, c, k;
 176     int cn = a->cn;
 177     int rn = a->rn;
 178     double *e = a->e;
 179 
 180     for (k = rn; k < cn; k++) { /* pętla po prawych stronach */
 181       for (r = rn - 1; r >= 0; r--) {   /* petla po niewiadomych */
 182         double rhs = *(e + r * cn + k); /* wartość prawej strony */
 183         for (c = rn - 1; c > r; c--)    /* odejmujemy wartości już wyznaczonych niewiadomych *
 184                                            pomnożone przed odpowiednie współczynniki */
 185           rhs -= *(e + r * cn + c) * *(e + c * cn + k);
 186         *(e + r * cn + k) = rhs / *(e + r * cn + r);    /* nowa wartość to prawa strona / el. diagonalny */
 187       }
 188     }
 189     return 0;
 190   }
 191   else
 192     return 1;
 193 }

pivot.c

   1 #include "matrix.h"
   2 #include <stdlib.h>
   3 #include <math.h>
   4 
   5 matrix_t *
   6 pivot_ge_matrix (matrix_t * a, int *row_per)
   7 {
   8   matrix_t *c = copy_matrix (a);
   9   if (c != NULL) {
  10     int i, j, k;
  11     int cn = c->cn;
  12     int rn = c->rn;
  13     double *e = c->e;
  14     for (i = 0; i < rn; i++)
  15       row_per[i] = i;
  16     for (k = 0; k < rn - 1; k++) {      /* eliminujemy (zerujemy) kolumnę nr k */
  17       int piv = k;              /* wybór eleemntu dominującego - maks. z k-tej kol., poniżej diag */
  18       for (i = k + 1; i < rn; i++)
  19         if (fabs(*(e + i * cn + k)) > fabs(*(e + piv * cn + k)))
  20           piv = i;
  21       if (piv != k) {           /* jeśli diag. nie jest pivtem - wymień wiersze */
  22                                 int tmp;
  23         xchg_rows (c, piv, k);
  24                                 tmp= row_per[k];
  25         row_per[k] = row_per[piv];
  26         row_per[piv] = tmp;
  27       }
  28       for (i = k + 1; i < rn; i++) {    /* pętla po kolejnych
  29                                            wierszach poniżej diagonalii k,k */
  30         double d = *(e + i * cn + k) / *(e + k * cn + k);
  31         for (j = k; j < cn; j++)
  32           *(e + i * cn + j) -= d * *(e + k * cn + j);
  33       }
  34     }
  35   }
  36   return c;
  37 }
  38 
  39 matrix_t *
  40 symm_pivot_ge_matrix (matrix_t * a, int *row_per)
  41 {
  42   matrix_t *c = copy_matrix (a);
  43   if (c != NULL) {
  44     int i, j, k;
  45     int cn = c->cn;
  46     int rn = c->rn;
  47     double *e = c->e;
  48     for (i = 0; i < rn; i++)
  49       row_per[i] = i;
  50     for (k = 0; k < rn - 1; k++) {      /* eliminujemy (zerujemy) kolumnę nr k */
  51       int piv = k;              /* wybór elementu dominującego - maks. z k-tej kol., poniżej diag */
  52       for (i = k + 1; i < rn; i++)
  53         if (fabs(*(e + i * cn + k)) > fabs(*(e + piv * cn + k)))
  54           piv = i;
  55       if (piv != k) {           /* jeśli diag. nie jest pivtem - wymień wiersze */
  56                                 int tmp;
  57         xchg_rows (c, piv, k);
  58         xchg_cols (c, piv, k);
  59                                 tmp= row_per[k];
  60         row_per[k] = row_per[piv];
  61         row_per[piv] = tmp;
  62       }
  63       for (i = k + 1; i < rn; i++) {    /* pętla po kolejnych
  64                                            wierszach poniżej diagonalii k,k */
  65         double d = *(e + i * cn + k) / *(e + k * cn + k);
  66         for (j = k; j < cn; j++)
  67           *(e + i * cn + j) -= d * *(e + k * cn + j);
  68       }
  69     }
  70   }
  71   return c;
  72 }
  73 
  74 int *
  75 pivot_get_inv_per ( matrix_t *m, int *row_per)
  76 {
  77         /* odtwarzamy oryginalną kolejność niewiadowmych */
  78   int *iper = malloc( m->rn * sizeof *iper );
  79         int i;
  80 
  81         for( i= 0; i < m->rn; i++ )
  82                 iper[row_per[i]] = i;
  83         
  84   return iper;
  85 }

ge_solver.c

   1 #include "matrix.h"
   2 #include <stdio.h>
   3 
   4 int
   5 main (int argc, char **argv)
   6 {
   7   FILE *in;
   8   if (argc > 1 && (in = fopen (argv[1], "r")) != NULL) {
   9     matrix_t *m = read_matrix (in);
  10     if (m != NULL) {
  11       matrix_t *c = NULL;
  12       printf ("\nMacierz:\n");
  13       write_matrix (m, stdout);
  14       c = ge_matrix (m);
  15       if (c != NULL) {
  16         printf ("\nPo elim. Gaussa:\n");
  17         write_matrix (c, stdout);
  18         if (bs_matrix (c) == 0) {
  19           int i, j;
  20           printf ("\nPo podstawieniu wstecz:\n");
  21           write_matrix (c, stdout);
  22           printf ("Rozwiązania: \n");
  23           for (j = 0; j < c->cn - c->rn; j++) {
  24             printf ("Nr %d:\n", j + 1);
  25             for (i = 0; i < c->rn; i++)
  26               printf ("\t%g", *(c->e + i * c->cn + j + c->rn));
  27             printf ("\n");
  28           }
  29         }
  30       }
  31       else
  32         printf ("nie lula!\n");
  33     }
  34     return 0;
  35   }
  36   else
  37     return 1;
  38 }

pivot_ge_solver.c

   1 #include "matrix.h"
   2 #include <stdio.h>
   3 #include <string.h>
   4 #include <stdlib.h>
   5 
   6 int
   7 main (int argc, char **argv)
   8 {
   9   FILE *in;
  10   if (argc > 1 && (in = fopen (argv[1], "r")) != NULL) {
  11     matrix_t *m = read_matrix (in);
  12     int sym= 0;
  13     if (m != NULL) {
  14       matrix_t *c = NULL;
  15       int *row_per = malloc (m->rn * sizeof *row_per);
  16       printf ("\nMacierz:\n");
  17       write_matrix (m, stdout);
  18       if( argc > 2 && strcmp( argv[2], "-s" )==0 ) {
  19         c = symm_pivot_ge_matrix (m, row_per);
  20                                 sym= 1;
  21                         } else
  22         c = pivot_ge_matrix (m, row_per);
  23       if (c != NULL) {
  24         int i;
  25         printf ("\nPo elim. Gaussa:\n");
  26         write_matrix (c, stdout);
  27         printf ("Permutacja:");
  28         for (i = 0; i < c->rn; i++)
  29           printf (" %d", row_per[i]);
  30         printf ("\n");
  31         if (bs_matrix (c) == 0) {
  32           int j;
  33           int *iper = pivot_get_inv_per (c, row_per);
  34           printf ("Permutacja odwrotna:");
  35           for (i = 0; i < c->rn; i++)
  36             printf (" %d", iper[i]);
  37           printf ("\n");
  38           printf ("\nPo podstawieniu wstecz:\n");
  39           write_matrix (c, stdout);
  40           printf ("Rozwiązania: \n");
  41           for (j = 0; j < c->cn - c->rn; j++) {
  42             printf ("Nr %d:\n", j + 1);
  43             for (i = 0; i < c->rn; i++) {
  44               int oi = sym ? iper[i] : i;
  45               printf ("\t%g", *(c->e + oi * c->cn + j + c->rn));
  46             }
  47             printf ("\n");
  48           }
  49         }
  50       }
  51       else
  52         printf ("nie lula!\n");
  53     }
  54     return 0;
  55   }
  56   else
  57     return 1;
  58 }

Makefile - tak, mówiliśmy też troszeczkę o programie make

CFLAGS = -Wall --ansi -pedantic

solver: ge_solver.o matrix.o
        $(CC) $(CFLAGS) -o $@ $^

test-solver: solver a56
        ./solver a56

pivot: pivot_ge_solver.o matrix.o pivot.o
        $(CC) $(CFLAGS) -o $@ $^

test-pivot: pivot a1011
        ./pivot a1011
        ./pivot a1011 -s

mull: test_mull_matrix.o matrix.o
        $(CC) $(CFLAGS) -o $@ $^

test-mull: mull a5 v5
        ./mull a5 v5

generator: mk_seq.c
        $(CC) $(CFLAGS) -o seq mk_seq.c

.PHONY: clean
clean:
        -rm *.o *~

Na koniec kilka plików z danymi: a5

5 5
 4 -1 -1  0  0 
-1  4 -1 -1  0 
-1 -1  4 -1 -1
 0 -1 -1  4 -1
 0  0 -1 -1  4

v5

5 1
50
40
30
20
10

a1011

10 11
0.840188 0.394383 0.783099 0.79844 0.911647 0.197551 0.335223 0.76823 0.277775 0.55397 29.4476
0.477397 0.628871 0.364784 0.513401 0.95223 0.916195 0.635712 0.717297 0.141603 0.606969 32.6739
0.0163006 0.242887 0.137232 0.804177 0.156679 0.400944 0.12979 0.108809 0.998925 0.218257 20.2714
0.512932 0.839112 0.61264 0.296032 0.637552 0.524287 0.493583 0.972775 0.292517 0.771358 33.1302
0.526745 0.769914 0.400229 0.891529 0.283315 0.352458 0.807725 0.919026 0.0697553 0.949327 33.4921
0.525995 0.0860558 0.192214 0.663227 0.890233 0.348893 0.0641713 0.020023 0.457702 0.0630958 15.8318
0.23828 0.970634 0.902208 0.85092 0.266666 0.53976 0.375207 0.760249 0.512535 0.667724 32.8602
0.531606 0.0392803 0.437638 0.931835 0.93081 0.720952 0.284293 0.738534 0.639979 0.354049 31.8288
0.687861 0.165974 0.440105 0.880075 0.829201 0.330337 0.228968 0.893372 0.35036 0.68667 30.7581
0.956468 0.58864 0.657304 0.858676 0.43956 0.92397 0.398437 0.814767 0.684219 0.910972 39.8569

a56

5 6
 4 -1 -1  0  0  -1
-1  4 -1 -1  0   0
-1 -1  4 -1 -1   0
 0 -1 -1  4 -1   6
 0  0 -1 -1  4  13

pivot1 - układ równań, którego nie da się rozwiązać bez wyboru elementu głównego

3 4
1 -1  2   8 
0  0 -1 -11
0  2 -1  -3

Na sam koniec powiedzieliśmy też kilka słów o bardziej efektywnym przechowywaniu macierzy pod kątem wyboru elementu głównego i zamiany wierszy - ale o tym porozmawiamy jeszcze raz na kolejnym wykładzie.


2015-09-23 06:44